Withasomniferol C, a new potential SARS-CoV-2 main protease inhibitor from the Withania somnifera plant proposed by in silico approaches

Collection of phytocompounds

The collection and identification of withanolides from the Withania somnifera plant was manually done through literature curation. The NCBI PubMed database search engine (https://www.ncbi.nlm.nih.gov/pubmed/) was used to collect peer-reviewed research articles. All research articles with the search term “withania somnifera” were collected for the analysis. Each article was reviewed separately, and the phytocompounds reported in the WS plant with valid experimental evidence were selected for the analysis. The structure information of the phytocompounds was retrieved from the PubChem database in the SMILES format.

Molecular fingerprinting

A structure based molecular screening of the collected phytocompounds was performed against compounds with a determined SARS-CoV-2 inhibitory activity. A structure based molecular screening was performed using 2D molecular fingerprints (FP) in RDKit (2020.03.1). In the present study, the Morgan circular FP (i.e., the extended connectivity) (Rogers & Hahn, 2010) was used to generate the FP of the phytocompounds (the FP radius was set to four). The molecular similarity exploration was performed using the Tanimoto coefficient (Tc) (Bajusz et al., 2015). The Tc similarity score ranges from zero to one, with zero representing the minimum and one representing the maximum similarity.

The generated FPs were compared to the calculated FP of 8702 molecules collected from the ChEMBL database (Davies et al., 2015) (release version ChEMBL 27) with the assay ID CHEMBL4495582, and the target ID CHEMBL4523582. The assay ID CHEMBL4495582 corresponds to the SARS-CoV-2 inhibitors and target ID CHEMBL4523582 corresponds to the replicase polyprotein 1ab of the SARS-CoV-2.

Molecular docking

Detailed investigation of phytocompounds against SARS-CoV-2 3CL^pro was performed by molecular docking studies. Before docking, 3D conformations of phytocompounds were generated and optimized from their SMILES notations using RDKit and MMFF94 force field. The 3D structures were then converted into the pdbqt file format using AutoDockTools 4 script prepare_ligand4.py (Morris et al., 2009). Similarly, the receptor was prepared as described below using UCSF Chimera 1.14 (Pettersen et al., 2004). In the present study, we used the SARS-CoV-2 3CL^pro conformation which was generated in our previous study based on the k-means clustering of the unbound 3CL^pro during a course of a 900 ns MD simulations (PDB ID: 6LU7) (Novak et al., 2021a). The k-means clustering was based on all non-hydrogen backbone atoms, and it gave two different conformations of the SARS-CoV-2 3CL^pro, the primary A conformation and the secondary B conformation, which were present 86.7% and 13.3% of the simulation time, respectively. The structural differences of these conformations compared to the crystallographic conformation of the SARS-CoV-2 3CL^pro (6LU7) are thoroughly discussed in (Novak et al., 2021a), while for the present work, the A conformation was chosen, as it is the more prevalent one. The A conformation of the SARS-CoV-2 3CL^pro was prepared by adding the Gasteiger charges to all atoms, followed by merging of the nonpolar hydrogens. Such structure was then saved in the pdbqt file format. Docking was performed using AutoDock Vina (Trott & Olson, 2010a) locally on a personal computer with 8 Intel^® CoreTM i7-6700K CPU @ 4.00 GHz, 32 GB RAM, and 64-bit Windows 10 Pro operating system. The center of the grid box was set on the Cys145 CA atom (dimensions: x = 13.3, y = 58.2, z = 45.4), with a box size of 20  ×25  × 25 Å. The grid spacing was set to 1 Å and number of modes and exhaustiveness were both set to 100.

The phytocompounds’ physicochemical descriptors, ADME properties, and druglike nature were evaluated using the SwissADME server (Daina, Michielin & Zoete, 2017).

Selection of docked poses using the AlteQ method

Validation of the conformations generated by the docking tools can be assessed using different methods (B-H & Brenk, 2009), among which the RMSD (root-mean-square deviation) based methods are the most widely used. A docked ligand pose with an RMSD score less than 2 Å, compared to the reference ligand (the crystallographically determined ligand pose), is generally considered as a good pose (Trott & Olson, 2010b). However, selection of docked poses based solely on RMSD has several flaws and can lead to misclassification of both the correct and incorrect poses (Cole et al., 2005; Kroemer et al., 2004). In the current work, we used the SARS-CoV-2 3CL^pro conformation from a molecular dynamics simulation, and hence measuring the correctness of docking poses based on the RMSD with the reference model is not an appropriate method. Besides, the covalently bound peptidomimetic ligand of SARS-CoV-2 3CL^pro (6LU7) was used for the MD simulations and measuring the RMSD between covalently bound peptidomimetic ligand with the non-covalently docked ligand is also not appropriate. Therefore, in the present work distance-based measures (i.e., ligand–receptor contacts) were used to check for the correctness of the docked poses. Many distance-based measures have been developed (Rueda et al., 2010; Hawkins et al., 2008), and in these methods, the cut-off length for the ligand–receptor contacts has to be defined. The complementarity principle coupled with the AlteQ method is a newly developed method for determining ligand–receptor contacts in which no cut-offs have to be introduced (Potemkin & Grishina, 2008); it calculates ligand–receptor contacts based on the electron density overlaps between the ligand and the receptor atoms using the Slater’s type atomic contributions (Potemkin, Grishina & Potemkin, 2017; Grishina & Potemkin, 2019). Recently, we used this method to calculate electron density overlaps in EGFR (Kandagalla et al., 2021) and CDK (Rimac, Grishina & Potemkin, 2020) receptor–ligand complexes, where we showed that all interactions determined by overlaps of electron clouds follow the complementarity principle, expressed by the following equation: (1)${\displaystyle ln\left({\rho }_{ligand}\times {\rho }_{enzyme}\right)=b+a\times RLRE.}$

where ρ_ligand represents ligand’s contribution to electron density in the m^th point of the molecular space, ρ_enzyme represents enzyme’s contribution to electron density in the same point, and RLRE is defined in a following manner (Eq. 2): (2)${\displaystyle RLRE=dis{t}_{enzyme}\times dis{t}_{ligand}}$ where dist_ligand (dist_enzyme) represents the distance between the mth point and the ligand’s (enzyme’s) atom having the highest contribution to ρ_ligand (ρ_enzyme) at that point.

The above complementarity model (Eq. 1) was used for validation of the docking poses. In this regard, the complementarity model was developed from the experimentally solved crystal structure of the SARS-CoV-2 3CL^pro. A total of 29 experimentally solved SARS-CoV-2 3CL^pro crystal structures with non-covalently bound inhibitors were selected from the RCSB Protein Data Bank (https://www.rcsb.org/) (Berman et al., 2002) for the analysis (PDB ID’s 5R7Y, 5R7Z, 5R80, 5R81, 5R82, 5R83, 5R84, 5RE4, 5RE9, 5REB, 5REH, 5REZ, 5RF1, 5RF2, 5RF3, 5RF6, 5RF7, 5RFE, 5RG1, 5RGH, 5RGI, 5RGK, 5RGU, 5RGV, 5RGW, 5RGX, 5RGY, 5RGZ, 5RH0, 5RH1, 5RH2, 5RH3, 5RH8, 5RHd, 6M2N, 6W79, 7JU7, 7KX5, and 7L5D). On the other hand, the five conformations of each phytocompound showing the lowest binding energy with 3CL^pro were collected from docking studies for further analysis. All the collected conformations were prepared using the UCSF Chimera 1.14 (University of California, USA) (Pettersen et al., 2004). The chain A conformations with the ligands were retained for the electron density analysis. The 3D maps of the electron density (ρ) were calculated for all conformations using the in-house developed quantum free-orbital AlteQ method (Potemkin & Grishina, 2008). The AlteQ method represents molecular electron density as a sum of Slater’s type atomic increments, and it can be expressed as (Eq. 3): (3)${\displaystyle \rho \left(x_m,y_m,z_m\right)=\sum _{A=1}^N{\rho }_A}$ (4)${\displaystyle {\rho }_A=\sum _{i=1}^{n_A}{a}_{A_{isp}}exp\left(-b_{A_{isp}}{R}_A\right)+\sum _{i=3}^{n_A-1}{a}_{A_{id}}exp\left(-b_{A_{id}}{R}_A\right)+\sum _{i=4}^{n_A-2}{a}_{A_{if}}exp\left(-b_{A_{if}}{R}_A\right)}$ where N is the number of atoms in the molecule, ρA is the A atomic increment in molecular electron density, a_Aisp, b_Aisp, a_Aid, b_Aid, a_Aif, and b_Aif are the AlteQ atomic parameters describing the i-th sp-orbital, d-orbital, and the f-orbital of the A atom respectively, nA is the period number of the A atom, RA is the distance between the A atomic center and point m. The units of the AlteQ coefficients are [b] = 1/Å, [a] = e/Å3, and consequently [ρ] = e/Å3.

Thus, the electron density of the outer shell, which plays the most important role in the formation of covalent bonds and intermolecular contacts, can be calculated at each point of the molecular space as follows: (5)${\displaystyle {\rho }_{A\left(outer\right)}=a_{A_{nsp}}exp\left(-b_{A_{nsp}}{R}_A\right)+a_{A_{\left(n-1\right)d}}exp\left(-b_{A_{\left(n-1\right)d}}{R}_A\right)+a_{A_{\left(n-2\right)f}}exp\left(-b_{A_{\left(n-2\right)f}}{R}_A\right).}$

The AlteQ method separates the molecular density into atomic contributions and allows one to separately consider contributions of the enzyme and the ligand atoms (Eq. 5). The generated electron density maps were used to construct a linear regression model to establish a correlation between the distance and the electron density overlap between the ligand and receptors atoms using Eqs. (1) and (2). The generated 3D electron density maps were statistically processed using the scikit-learn (Pedregosa et al., 2011) and the SciPy library (Virtanen et al., 2020) in Python (v 3.7.6), and the plots were generated using the Matplotlib library (Hunter, 2007). The AlteQ software is freely available online via the ChemoSophia webpage (http://www.chemosophia.com/).

Molecular dynamics (MD) simulations

MD simulations for all ligands (WS1, WS4, WS11, and two conformations of the WS7 ligand) were run in complex with the SARS-CoV-2 3CL^pro in the previously acquired conformation (PDB: 6LU7). The best docked ligand positions, which were determined using the AlteQ method, were used as the starting points for the MD simulations. The AMBER ff14SB force field (Maier et al., 2015) was used to model the protein and the GAFF force field (Wang et al., 2004) was used to model the ligand. Other simulation parameters such as periodic boundary conditions, NVT conditions, pressure, temperature were the same as described in our previous article (Novak et al., 2021b; Pathak et al., 2021). In brief, protein-ligand complexes were solvated in a truncated octahedral box of TIP3P water molecules spanning a 12 Åthick buffer, and Na⁺ and Cl⁻ ions were added according to (Machado & Pantano, 2020) to achieve a neutral environment with a salt concentration of 0.15 M. Such structures were then submitted for geometry optimization in the AMBER16 program (Case et al., 2016) employing periodic boundary conditions in all directions. For the first 1,500 cycles, the complex was restrained and only water molecules were optimized, after which another 2,500 cycles of optimization followed where both water molecules and the complex were unrestrained. Optimized systems were gradually heated from 0 to 310 K and equilibrated during 30 ps using NVT conditions, followed by productive and unconstrained MD simulations of 300 ns employing a time step of 2 fs at constant pressure (1 atm) and temperature (310 K), the latter held constant using Langevin thermostat with a collision frequency of 1 ps⁻¹. Bonds involving hydrogen atoms were constrained using the SHAKE algorithm (Ryckaert, Ciccotti & Berendsen, 1977), while the long-range electrostatic interactions were calculated employing the Particle Mesh Ewald method (Darden, York & Pedersen, 1993). The non-bonded interactions were truncated at 11.0 Å. Analysis of the trajectories was performed using the cpptraj module of AmberTools16 (Roe & Cheatham, 2013).

Binding free energy calculations and decomposition

The binding energies, ΔG_BIND, of the simulated complexes were calculated using the MM-GBSA (Molecular Mechanics–Generalized Born Surface Area) and the MM-PBSA (Molecular Mechanics–Poisson–Boltzmann surface area) protocols (Genheden & Ryde, 2015; Hou et al., 2011), available as a part of AmberTools16 (Ferenczy, 2015). ΔG_BIND is calculated from snapshots of MD trajectories (Ferenczy, 2015) with an estimated standard error of 1–3 kcal/mol (Genheden & Ryde, 2015). ΔG_BIND is calculated in the following manner: ${\displaystyle \Delta G_{BIND}=<G_{complex}>-<G_{protein}>-<G_{ligand}>}$ where the symbol <> represents the average value over 100 snapshots collected from the last 30 ns part of the corresponding MD trajectories (every 150th frame was taken for the calculation). The calculated MM-GBSA and MM-PBSA binding free energies were decomposed into specific residue contribution on a per-residue basis according to established procedures. This protocol calculates the contributions to ΔG_BIND arising from each amino acid side chains and identifies the nature of the energy change in terms of interaction and solvation energies (Gohlke, Kiel & Case, 2003; Rastelli et al., 2010). The entropy term was not calculated.

Collection of Withania somnifera compounds from literature and molecular fingerprinting

A search of the NCBI PubMed database resulted in 1,401 research articles which included the term “withania somnifera”, with the dates until February 28, 2021. All the articles were reviewed and withanolides were collected. Altogether, a total of 80 compounds were collected. However, several studies have already reported WS withanolides against the SARS-CoV-2 3CL^pro andhence these compounds were excluded from the analysis. A total of 11 specific types of withanolides, such as withaphysalins, withasomniferols, and withafastuosins, were considered for a detailed binding interaction since there were no previous studies focused on their binding interactions with the SARS-CoV-2 3CL^pro. Details about withanolides selected in the present investigation and their respective article information are provided in File S1. Among the selected withanolides, withasomniferol A, withasomniferol B, and withasomniferol C are newly isolated WS phytoconstituents (Anjaneyulu & Rao, 1997).

Similarity of the collected phytocompounds with reported SARS-CoV-2 inhibitors was checked using a molecular similarity analysis. Tc scores of the phytocompounds were calculated relative to the reported SARS-CoV-2 inhibitors, and they ranged from zero to 0.4 (File S2).

Molecular docking

From all tested phytocompounds, compounds with ID’s WS1, WS4, WS7, and WS11 showed the lowest energies of −8.0, −8.2, −7.6, and −7.8 kcal/mol, respectively, when binding to the SARS-CoV-2 3CL^pro (File S3). Their molecular interactions with the SARS-CoV-2 3CL^pro active pocket residues are shown in Figs. 1 and 2, and 2D structures of phytocompounds are shown in Fig. 3. Their binding energies against the 3CL^pro and their ADME properties are shown in Table 1. The analysis of the ADME properties of potential drug candidates is essential in the early stage of drug discovery to reduce failure rates in the clinical phase of drug discovery. According to Lipinski’s rule of 5 (RO5), to be orally active, drug-like compounds should have molecular weight below 500, number of hydrogen bond donors below 5, number of hydrogen bond acceptors below 10, and the log P value should not exceed 5. As shown in Table 1, all the selected compounds obey these rules except for the WS1 compound, which has a molecular weight above 500. All four compounds were found to be non-inhibitors of the cytochrome P450 enzymes (CYP3A4, CYP2D6, CYP2C9, CYP2C19, and CYP1A2). The cytochrome P450 enzymes play an essential role in metabolism of various molecules. Their inhibition can cause serious drug-drug interactions that can cause unanticipated adverse effects (Lynch & Neff, 2007). Since all our lead compounds were predicted to be non-inhibitors of the cytochrome P450 enzymes, this suggests that they possess a low potential for drug-drug interactions. Additionally, none of the suggested lead compounds are Pan-assay interference compounds (PAINS), which indicates that they will not give false positive results in high-throughput screens. All PAINS interact nonspecifically with numerous biological targets, as opposed to with specific targets. Furthermore, stability of these compounds with SARS-CoV-2 3CL^pro was verified through MD simulations. Additionally, before the MD analysis, poses of the compounds generated by docking were evaluated through the complementarity centered AlteQ method.

Selection of the docked poses with the complementarity-centered AlteQ method

Firstly, a linear regression model was constructed to establish a correlation between the distance between the ligand and receptors atoms and their electron density overlap (Eq. 1) for the experimental SARS-CoV-2 3CL^pro complexes and this information was used for evaluating the correctness of the docked conformations. The correlation coefficients (R²) in the 5R7Y, 6M2N, 5RGH, 5RF7, 5RG1, 5RH8, and 5REB complexes were found to be below 0.50 and hence were not included in model generation. The average R² of all the other experimental conformations was 0.732 (File S4). In a recent article, we described two ways of using the complementary based AlteQ method for validating docking scores (Rimac, Grishina & Potemkin, 2021). In the present work, we evaluated the correctness of the docked poses from the intercept-versus-slope graph. The intercept and the slope coefficients for all experimental complexes collected from Eq. 1 are correlated and fall onto the same regression line (R² = 0.97). The following model (Eq. 6) was obtained (Fig. 4). (6)${\displaystyle a=-2.90-0.33\times b\left(R^2=0.97\right)}$ where, a is the intercept value and b is the slope value.

The distance between the receptor and the ligand decreases with the slope, as can be observed in Fig. 4. In the 5RH1 complex, the ligand and the receptor atoms are very close, and the slope value was found to be −1.83. In contrast, the slope value of the 5R81 complex was found to be very high (−0.88), which indicates long distances between the ligand and the receptor atoms. The regression line obtained from the experimental SARS-CoV-2 3CL^pro complexes was used a reference line to measure appropriateness of the generated docked conformations.

The top five docking poses of WS1, WS4, WS7, and WS11 were evaluated, as they show the lowest binding energies when binding to the 3CL^pro. Their slope and intercept values were collected from Eq. 1 and plotted against the experimental conformations (Fig. 4). As expected, different electron density overlap patterns were observed in the docked conformations; they are slightly away from the regression line of the experimental conformations due to different conformations and their orientation within the active pocket, and they correspond to different interaction patterns. The distance of the docked conformations from the regression line was measured and the complete information is provided in File S5. In the top five conformations, the distance of the conformation 1 (based on the docking results) of WS1(0.062), WS4 (0.045), WS11(0.105) ligands were found to be much lower compared to the other four conformations, which is in accordance with the experimental conformations. The only exception was conformation 1 of the WS7 ligand, whose distance was found to be much greater (0.111). The electron density overlap patterns of a few conformations of the ligands WS4, WS7, and WS11 were more distant from the rest of the points (Fig. 4, dotted circle), and these conformations show a unique pattern in the electron density overlap, which was not observed in the experimental conformations.

Molecular dynamics simulations and free energy calculations

To check for the stability of the conformations, MD simulations were performed. For ligands WS1, WS4, and WS11 the best docked conformations were chosen since their electron density overlap patterns were very similar to the electron density overlaps of the experimental conformations. For the WS7 ligand, two different conformations were considered, namely the first (−7.6 kcal/mol, WS7_v1) and the fifth conformation (−7.3 kcal/mol, WS7_v2), with both conformations having different electron density overlap patterns with different binding energies. The distance of the conformation 1 of the WS7 ligand (0.111) was found to be much higher when compared to the fifth conformation (0.062) from the regression line, but the fifth conformation showed a unique pattern in the electron density overlap, and its points (i.e., slope and intercept values from Eq. 3) are located far away from the rest of the points (Fig. 4, dotted circle).

From the five complexes tested, three complexes were found to be stable for the entire duration of the simulation (WS4, WS7_v1, and WS11), while two complexes showed a brief dissociation (∼8.5 ns in the case of WS7_v2 and ≈ 26.0 ns in the case of WS1), which can also be seen from Fig. 5. The root-mean-square fluctuation (RMSF), radius of gyration and intermolecular hydrogen bond plots are provided in the File S6.

From Table 2, it can also be seen that WS1 and WS7_v2 bind the weakest (calculated using the MM-GBSA approach). From the three stable compounds (WS4, WS7_v1, and WS11), WS7_v1 and WS11 have approximately the same ΔG_BIND, which is significantly lower than that of WS4. Additional calculations of ΔG_BIND were performed using the MM-PBSA protocol. These results are completely in accordance with the MM-GBSA results and show the same trend; WS7_v1 and WS11 have virtually the same binding constant (∼−26.00 kcal/mol), with the other three ligands having at least 12 kcal/mol higher ΔG_BIND (detailed information can be found in File S7). The top ten contributing amino acids for the binding of the two best ligands (i.e., WS7_v1 and WS11) are shown in Table 3 and their positions relative to the catalytic Cys145 residue are shown in Fig. 6. Relative to the figure, the WS7_v1 ligand is located slightly more to the left, and the WS11 ligand slightly more to the right. This results in a fact that these ligands share nine out of ten most contributing residues (Table 3), with exception being Thr45, which interacts only with the WS7_v1 far left hydroxy group, and Met165, which interacts only with the WS11 cyclopentanyl group (Fig. 7, Table 3, bolded). While the hydrogen bond interaction with Thr45 is moderately important in case of WS7_v1 binding, the van der Waals interaction with Met165 seems to be the most important one for WS11_v1 binding. Additionally, the WS7_v1 ligand forms two very strong interactions with partial ΔG_BIND lower than −2.0 kcal/mol, with one of them being the Cys145 residue, a crucial residue for the 3CL^pro function (Cole et al., 2005), and the other one being Thr25 which forms a very strong hydrogen bond with the keto oxygen atom of WS7. On the other hand, WS11 forms much stronger interactions with the other crucial residue, (His41), as well as with Cys44, with both interactions being hydrophobic in nature, while in the case of WS7_v1 these interactions have a van der Waals character. This indicates that WS7 and WS11 interactions have a preference for different parts of the binding pocket and that they could both be good candidates for 3CL^pro inhibitors. This also implies that there is still room for developing or finding a ligand which could interact more strongly with both parts of the binding pocket at the same time. The complete binding free energy results of all the ligands are provided in the File S7.