Repurposing FDA-approved phytomedicines, natural products, antivirals and cell protectives against SARS-CoV-2 (COVID-19) RNA-dependent RNA polymerase

Construction of FDA approved drugs dataset

A total of 1,664 FDA approved drugs were retrieved from commercial compounds databases at Selleck Inc., USA. All compounds were imported to Ligprep software, desalted and 3D optimized using OPLS2005 force field at physiological pH. The compounds dataset are provided in Supplemental Files 1.

SARS-CoV-2 RdRP optimization and structure preparation

The protein preparation module in Maestro software package (Schrodinger LLC, NY, USA) was used to optimize the protein structure for docking. The protein (PDB ID 7BV2) was protonated, optimized at cellular pH conditions, and energy minimized using OPLS2005 force field. All ions, except the active site magnesium, and non-relevant crystallographic materials were removed. During preparation, two structural files were prepared comprising RdRP in the presence or absence of magnesium ions. The docking grid was generated by the selection of active site residues. The docking box was centered around the active site residues ASP255 and ASP256 with a size of 20 Å. For comparative purpose, the sequence similarity between SARS CoV, SARS CoV-2 and MERS CoV RdRPs were retrieved and aligned (Supplemental Methods File). The reliability of docking procedures were assigned by redocking of the co-crystalized inhibitor (remdesivir). The obtained poses were comparable to the co-crystalized structure (Fig. 1).

Virtual screening

The glide docking module of Schrodinger maestro package was used for virtual screening. Two-step docking protocol was carried out. At first, the docking was set to run under the standard precision (SP docking) (Friesner et al., 2004; Halgren et al., 2004). The results were ranked by the docking score. Remdesivir was used as standard RdRP inhibitor. It is a broad-spectrum antiviral with potent inhibition of MERS CoV RdRP (Gordon et al., 2020). Compounds with a docking score of −7.00 or higher were retrieved and reanalyzed by extra-precision docking protocol (XP-docking). Magnesium ions are important in coordinating the polymerase actions of RdRP. Therefore, the two docking interactions were run comprising the presence and absence of Mg ions in their binding sites.

Molecular dynamics simulations

Groningen Machine for Chemical Simulations (GROMACS) 5.1.4. software was used in all molecular dynamics (MD) simulations (Abraham et al., 2015; Van Der Spoel et al., 2005). The ligand parameter, topology and restraint and the protein were handled by AMBERFF14SB force field and general AMBER force field (GAFF), respectively. The RdRP-ligand complexes were solvated in a cubic box with 1.0 nm from protein to box edge. The box was filled with a single point charge (SPC) water model. The solvated complexes were minimized for 5,000 steps. During water and ions coupling, the heavy atoms of protein and ligand were restrained. Two-step system equilibration was adopted comprising 50 ps NVT at 300 K followed by 1 ns NPT ensemble at 300 K. Production stages were done over simulation times of 20 and 100 ns with NPT ensemble. Parrinello–Rahman algorithm maintain constant pressure at 1 bar and V-rescale thermostat algorithmfor temperature coupling at 300 K. Long-range electrostatics under periodic conditions with a direct space cut-off 12 Å were restrained by Particle Mesh Ewald (PME) method. Time step was set to 2 fs. Coordinates and output values were collected every 10 ps.

Trajectory analysis

GROMACS MD simulation toolkits were used in the trajectory analysis. The g_rms and g_rmsf were used to calculate the root mean square deviation (RMSD) of the starting structure to the end of the simulation time and the per-residue root mean square fluctuation (RMSF) of the protein residues.

Binding energy calculations

The molecular mechanics–generalized Born surface area (MM/GBSA) was calculated by the g_mmpbsa tool.

Statistical analysis

Multiple correlation comparisons were performed between the obtained docking score and the properties and interactions of the compounds with SARS-CoV-2, including the molecular weight, the hydrogen bonding and the lipophilic interaction parameters. The analysis was carried out with the Pearson’s correlation coefficient, implemented in GraphPad Prism version 7 (GraphPad Software, La Jolla, CA, USA).

Virtual screening and docking

The virtual screening data for 1,664 FDA-approved drugs against SARS-CoV-2 are provided in Table S1. The presented parameters include the docking scores, ligand efficiencies, and lipophilic and hydrogen bonding interactions. The assessment of docking protocol was confirmed by docking of remdesivir monophosphate (RMP). The redocked pose was compared with the co-crystalized ligand (Fig. 1). It is was evident that the presence of magnesium ions increased the docking scores of RMP in both of the standard precision and extra-precision docking. In the absence of magnesium ions, the SP and XP docking scores of RMP were −8.9 and −9.7, respectively. These scores was decreased by about 15% in the absence of magnesium ions. In contrast, the docked ligands from the FDA approved drugs dataset showed higher scores in the absence of magnesium ions. This might be attributed to the favorable interaction of magnesium ions with the phosphorylated substrates.

Table 1 shows the 26 compounds with the highest docking scores, which were −7 or higher. Within the top hits, there were five general tonics, cell protectives, and antioxidants; one antiviral; one anti-parasitic; one antibiotic; one immune stimulant; three anticancer drugs; and other miscellaneous systemically or topically acting drugs.

A considerable number of natural products, natural product derivatives, and biological molecules were present among the top hits. Plant flavonoids (quercetin, troxerutin, taxifolin), plant glycosides (sennosides A and B, digoxin), a citrus-derived sweetener (neohesperidin), and biological molecules (oxiglutatione, nadid) constituted about 36% of the top binding drugs with SARS-CoV-2 RdRP.

In order to get deeper insights into the interaction of drugs with SARS-CoV-2 RdRP, the binding mode and ligand interaction profiles of the drugs and also of remdesivir were examined (Fig. 2). The active site residues comprised positively charged amino acids (LYS476, ARG507, LYS681 and GLU694), negatively charged amino acids (ASP328, ASP501, ASP506, ASP643, ASP644 and ASP748), aromatic residues (TRP461, TYR474, TRP500, TYR502, TRP683 and PHE726), polar amino acids (GLN327, ASN330, SER475, SER565, SER642, CYS696, SER697, GLN698 and PRO715) and non-polar amino acids (GLY329, ALA313, ALA331, ILE333, VAL471, ILE472, GLY473, MET484, LEU641 and ALA645). Naringen formed two hydrogen bonds with ARG555 and THR591 (Fig. 2A). In contrast, echinacoside showed multiple interaction routes, including six hydrogen bonds with ARG555, CYC622, ASP623, SER682, ASP760 and ARG836 (Fig. 2B). It was noticeable that the number of hydrogen bonds formed between the docked ligands with the active is affected by the presence of magnesium ions. While metal coordinates is evident with the docked poses, this was underscored by lower hydrogen bonds interactions (Fig. 2). This might contributed to the observed lower docking scores in structures containing magnesium.

Statistical analysis revealed a significant positive correlation between the docking score and the hydrogen bond, lipophilic interaction, and ligand efficiency scores (p < 0.001; Table 2). There was a negative low but significant correlation between the obtained docking scores and the molecular weights of the compounds. This explains the superior docking scores of the selected set of drugs over the standard inhibitor remdesivir, which involved the lower molecular weights and improved hydrogen bonding and lipophilic interactions.

For further filtration to identify the most potent hits, compounds with a score of −7.00 or higher were subjected to XP-docking followed by MM/GBSA calculations (Table 3). Six compounds produced high docking scores exceeding −11.00: echinacoside, salvianolic acid B, ginsenoside, neohesperidin, troxerutin, and docetaxel. These top six were then used for further MD simulations.

Molecular dynamics simulations for 20 ns

Combining docking with MD simulation is a powerful tool for drug discovery, as these techniques rank compounds on their binding affinity and can also evaluate target-receptor binding energy and drug-receptor dynamics. In this study, the top-ranked compounds from XP-docking were subjected to MD simulation followed by post dynamic analysis of binding energy. A two-step filtering process was adopted. First, the six compounds were run in MD simulation for 20 ns and their RMSD, RMSF, and MM/GBSA values evaluated. In the second stage, the top two compounds were tested in a more comprehensive 100 ns MD simulation. Structural changes in the backbone residues of RdRP relative to the initial structure (RMSD) were compared after 20 ns simulations (Fig. 3).

MM/GBSA binding energies at 20 ns

According to the MM/GBSA binding energies, strong binding was indicated for docetaxel (−69.19 kcal/mol) and neohesperidin (−65.754 kcal/mol). Ginsenoside showed a moderate binding energy of −47.52 kcal/mol (Table 4). The other three compounds yielded low binding energy values, all less than −30.708 kcal/mol. In addition, all of the tested compounds showed low structure RMSD during the 20 ns simulations, with values <0.24 nm (Table 4). Based on the 20 ns simulation data, docetaxel and neohesperidin were selected for further analysis in 100 ns MD simulations.

Molecular dynamics simulations for 100 ns

In order to obtain deeper insight into the strongest-binding drugs, the two with binding energy >−60 kcal/mol (docetaxel and neohesperidin) were subjected to 100 ns MD simulations followed by analysis of binding free energy, RMSD, RMSF, hydrogen bond length, and radius of gyration. The estimated MM/GBSA binding energies were −67.273 and −63.669 for docetaxel and neohesperidin, respectively. These values imply superior binding for both of the compounds.

The 100 ns simulations yielded average RMSD values of 0.25 and 0.22 nm for docetaxel and neohesperidin, respectively. These values indicate stability of RdRP in complexes with docetaxel or neohesperidin. Both RMSD and per-residue RMSF showed almost similar profiles for RdRP complexes with docetaxel and neohesperidin (Fig. 4).

Radius of gyration

Radius of gyration can be taken as a measure of the compactness of the system. Fig. 5 shows Rg variations observed during 100 ns MD simulations. For both compounds, the average Rg value was 3.18 nm. High Rg values indicate lower compactness or more unfolded protein, while low Rg indicates more stable structures. The similar Rg values indicate comparable stability of the two examined drugs when complexed with RdRP.

Decomposition of MM/GBSA binding energy at 100 ns

Post-dynamic energy decomposition analysis was used to evaluate the dominant interactions during recognition of the drugs by RdRP (Table 5). The results revealed that van der Waals and electrostatic interactions were the predominant forces underlying binding of the drugs to RdRP. Specifically, van der Waals interactions were the major forces for both docetaxel and neohesperidin, followed by electrostatic interactions.

The binding site of docetaxel and neohesperidin

The binding characteristics of docetaxel and neohesperidin with RdRP are provided in Fig. 6. The binding site is mostly composed of hydrophobic residues (ILE548, VAL557, ILE589, THR680, and ALA688), with the presence of few positively charged residues (LYS551, LYS545, and ARG555), negatively charged residues (ASP623 and ASP761), and neutral residues (SER549, SER682, ASN691, SER759, and SER814). The interaction of docetaxel and neohesperidin with RdRP was respectively supported by six and seven hydrogen bonds. Neohesperidin formed hydrogen bonding with ASP618, SER682, SER759, while docetaxel interacted with ILE548, SER549, LYS551 and SER814.

The number of hydrogen bonds

The number of hydrogen bonds between docetaxel and neohesperidin and RdRP was traced during MD simulation (Fig. 7; Fig. S3). The average number of hydrogen bonds were 4.3 and 6.2 for docetaxel and neohesperidin, respectively. Additionally, the distance between the two drugs and two selected residues in the active site was traced during simulation (Fig. 7B; Fig. S4). The selected residues were LYS521 for docetaxel and ASP588 for neohesperidin. After initial stabilization within the first few ns, the results indicated the lack of large fluctuations in the position of both drugs in the active site.