Study outline
Initially, 305 alkaloids constituents, coumarins, flavonoids, lignans, phenols, saponins, steroids, stilbenes, tannins and terpenoids obtained from PubChem were used as ligands for the current study. The absorption, distribution, metabolism and excretion (ADME) profiles of the above-mentioned study compounds were evaluated and screened out 195 ligands (Table S1) for docking study, based on Lipinski rule of five, where no violations (Molecular mass less than 500 Dalton, High lipophilicity (expressed as LogP less than 5), less than 5 hydrogen bond donors and less than 10 hydrogen bond acceptors. These properties which were followed as strategic process in drug discovery and development has been implemented in the current study [93]. Fig. 2 shows the workflow of the present study.
Target receptors selection
As stated in the introduction, a full understanding the life cycle mechanism of SARS-CoV-2 is necessary in order to identify druggable targets for the discovery and development of potent coronavirus therapeutics. In the present study, we have used the potential of natural compounds derived from medicinal plants targeting the SARS-CoV-2 can be divided into two categories, with the first group targeting virus-host interactions (ACE2 and TMPRSS2) or inhibiting viral enzymes (Mpro, PLpro and RdRp). These targets play a major role in the viral entry, replication / transcription, the many pharmacological target groups that can be addressed by the traditional de novo drug discovery strategy are reflected in genome and protein synthesis. Accordingly, release and assembly inhibitors are less explored steps in drug discovery against COVID-19 [94]. Despite, most of the drugs for 2019-nCoV are registered for control rather than prophylaxis; ACE2, integrin, and S protein targets are importance for drug repurposing or development programs aimed at preventing viral entrance and fusion. Considering the basic goal of all viruses is to transport and replicate their genetic code in competent host cells, inhibiting the Mpro, PLpro, angiotensin-converting enzyme 2 (ACE2), TMPRSS2, and RdRp are suggested to be promising targets for anti-2019-nCoV drug discovery of small molecules.
In addition, targeting human proteins is also a strategy to avoid viral scape caused by mutation. Another option is a multi-target strategy that combines the potential antiviral agents that acting on several targets to enhance efficacy and prevent virus resistance [95]. In light of the COVID-19 pandemic's urgency, the repurposing drugs is a better opportunity to develop a timely and effective cure. In fact, some of the pharmaceuticals currently being tested in clinical trials were originally approved for other purposes [96]. However, past and current coronavirus epidemics necessitate our preparedness not only for the immediate situation, but also for the possibility of novel coronaviruses re-emerging in the future. In this regard, finding equivalent medications that work as pancoronavirus antivirals or using a multi-target approach is crucial to minimize any loss of effectiveness due to viral mutation escape.
Molecular Docking Study
Total around 195 compounds were docked against five receptors including, main protease (Mpro) (PDB ID: 6LU7), Human Angiotensin Converting Enzyme (PDB ID: 1O86), primer RNA (Rdrp) (PDB ID: 7BV2), Trans membrane Protease Serine 2 (TMPRSS2) (PDB ID: 7MEQ) and Papain-like protease (PLpro) (PDB ID: 6WX4). The docking scores of all the study compounds were between -11.2 and -4.3 kcal/mol (Table S1). These compounds were compared with clinical drugs such as ritonavir, lopinavir, hydroxychloroquine, chloroquine, favipiravir, remdesivir, camostat, nafamostat and disulfiram (Table S2).
Docking results for main protease showed that 21 compounds showed stronger binding energy in the range of -9.1 and -7.9 kcal/mol as compared to control drugs, ritonavir (-7.1 kcal/mol) and lopinavir (-7.8 kcal/mol). As shown in Table 1, main protease amino acid residues MET6, ILE152, ASP153, ARG298, SER301 and GLY302 exhibited hydrogen bond with 3-HDH-withanolide F. Additionally, the alkyl interaction with PRO9 was also observed. These interactions were resulted in the binding energy of − 9.1 kcal/mol. The possible binding modes of 3-HDH-withanolide F at main protease active sites have been shown in Fig. 3.
Table 1
Docking scores and interaction modes of top hit selected antiviral compounds with COVID-19 Mpro receptor
No.
|
Compound
|
Score
Kcal/mol
|
Interactions
|
H
|
Hydrophobic
|
9
|
3-HDH-withanolide F
|
-9.1
|
MET6, ILE152, ASP153, ARG298, SER301, GLY302
|
Alkyl: PRO9
|
12
|
Limonin
|
-8.7
|
ARG131, LYS137, LEU287
|
π-Alkyl: LEU286, LEU287
|
2
|
Withanolide
|
-8.6
|
LYS137, ASN238, LEU286, LEU287
|
Alkyl: LEU286
|
16
|
(18R)-Withaphysalin F
|
-8.6
|
ARG40, TYR54, GLU55, PHE181
|
|
8
|
Obacunone
|
-8.6
|
THR199, GLY275, MET276
|
π-Alkyl: LEU286
|
89
|
Mimulone
|
-8.5
|
THR199, ASN238
|
π-Sigma: LEU272
Alkyl: VAL171, ALA194
|
1
|
Withanolide J
|
-8.4
|
TYR37, LYS88, TYR101
|
Alkyl: VAL35, LYS90
|
79
|
Inophyllum C
|
-8.4
|
GLN110, THR292, PRO293
|
π-Sigma: ILE249
π-π Stacked: PHE294
π-Alkyl: PRO293
|
5
|
Withanolide N
|
-8.3
|
MET276, ASN277, GLY278
|
π-Alkyl: TYR239
|
78
|
Inophyllum P
|
-8.3
|
ARG131, THR199, ASN238
|
|
On the other hand, ligands were docked at the active site of ACE target, the results demonstrated that 182 ligands could inhibit target receptors with a potential energy from -11.2 to -7.0 kcal/mol which were most potent than hydroxychloroquine (-6.8 kcal/mol) and chloroquine (-6.9 kcal/mol). It was found that withanolide exhibited highest docking score of -11.2 kcal/mol (Table 2). As shown in Fig. 4, the molecular interactions holding by withanolide in the active site resulted in the formation of five H-bond with amino acid residues namely GLN281, ASN374, GLU376, HIS383 and ASP415. Furthemore, four hydriphobic interactions were observed including two alkyl and π-alkyl interaction with PRO163, VAL380, TRP279 and HIS383 aminoacid residues.
Table 2
Docking scores and interaction modes of top hit selected antiviral compounds with COVID-19 ACE receptor
No.
|
Compound
|
Score
Kcal/mol
|
Interactions
|
H
|
Hydrophobic
|
Other
|
2
|
Withanolide
|
-11.2
|
GLN281, ASN374, GLU376, HIS383, ASP415
|
Alkyl: PRO163, VAL380,
π-Alkyl: TRP279, HIS383
|
|
16
|
(18R)-Withaphysalin F
|
-11.2
|
GLN281, THR282, ALA354, HIS513
|
Alkyl: ALA354, VAL380,
VAL518
π-Alkyl: HIS353, HIS383, HIS387, PHE512, TYR523
|
|
13
|
Withasomnilide
|
-11.1
|
ASN70, ARG124, ALA356, GLU143
|
Alkyl: LEU139, LEU140, VAL518
π-Alkyl: HIS353, TRP357, HIS387, PHE512, HIS513
|
|
89
|
Mimulone
|
-11.0
|
ALA356, GLU384, HIS387, HIS513
|
π-π T-shaped: TRP357
π-Alkyl: HIS353, HIS383, PHE457, HIS513, TYR523, PHE527
|
π-Sigma: PHE457
|
10
|
Somniferanolide
|
-11.0
|
GLY404, HIS410, GLU411, ARG522
|
Alkyl: MET223, PRO407,
π-Alkyl: TRP357, HIS387
|
π-Sigma: TRP357, HIS410
|
12
|
Limonin
|
-10.9
|
HIS353, SER355, TYR394, HIS410,
|
π-π T-shaped: HIS353
|
π-Cation: HIS353
|
1
|
Withanolide J
|
-10.8
|
ARG522, GLU411
|
Alkyl: VAL351, VAL518
π-Alkyl: HIS353, HIS387, PHE512
|
π-Sigma: HIS387
|
9
|
3-HDH-withanolide F
|
-10.7
|
ASN66, GLU411
|
Alkyl: LEU81, LEU140, VAL518
π-Alkyl: TYR69
|
|
14
|
7-Deacetylgedunin
|
-10.7
|
ALA354, ASP453, LYS454
|
Alkyl: VAL380
π-Alkyl: VAL379, PHE457
|
π-Sigma: GLU376
π-Anion: GLU376
|
6
|
Withasomniferanolide
|
-10.6
|
ASN70, SER355, ALA356, ARG522
|
π-Alkyl: TRP357, HIS387, HIS410
|
π-Sigma: TRP357, HIS387
|
Meanwhile, the minimum binding energies of docked ligands to the RdRp enzyme indicated that 50 of them (-7.2 to -8.2 kcal/mol) could strongly inhibit target protein as compared to the standard drugs, favipiravir and remdesivir with binding energy of -5.8 and -7.1 kcal/mol, respectively. The ligands were ranked according to their protein-ligand binding energies, aurantinidin was found most active compound with binding energy of -8.2 kcal/mol (Table 3). The interactions of aurantinidin and COVID-19 RdRp established five hydrogen bonds (ARG553, TRP617, ASP618, TYR619 and ASP760) to the amino acid residues. Moreover, another interaction was also observed with ASP618 and ASP760 residues via π-anion interaction type (Fig. 5).
Table 3
Docking scores and interaction modes of top hit selected antiviral compounds with COVID-19 RdRp receptor
No.
|
Compound
|
Score
Kcal/mol
|
Interactions
|
H
|
hydrophobic
|
Other
|
28
|
Aurantinidin
|
-8.2
|
ARG553, TRP617, ASP618, TYR619, ASP760
|
|
π-Anion: ASP618, ASP760
|
87
|
Tomentin A
|
-8.1
|
ASP618, ASP623, SER759, ASP760
|
π-Alkyl: CYS622, ALA688
|
π-Anion: ASP760
|
10
|
Somniferanolide
|
-8.0
|
ARG555, THR687, ASN691, ASP760
|
|
|
31
|
Europinidin
|
-7.9
|
TYR456, ARG555, THR556, SER759
|
π-Alkyl: ALA688
|
π-Cation: ARG555
π-Anion: ASP623
|
91
|
Scutellarein
|
-7.8
|
ASN691
|
π-Sigma: SER682
|
π-Cation: ARG555
π-Anion: ASP623
|
29
|
Cyanidin
|
-7.8
|
ASP452, TYR456, ARG555, THR556
|
|
π-Anion: ASP623
|
31
|
Herbacetin
|
-7.7
|
ASN691, ASP760
|
π-Alkyl: ALA558
π-Sigma: SER682
|
π-Cation: ARG555
π-Anion: ASP623
|
45
|
Pedalitin
|
-7.7
|
ASN691, SER681, SER682
|
π-Sigma: SER687
|
π-Cation: ARG555
π-Anion: ASP623
|
60
|
Silymarin
|
-7.7
|
THR556, THR680, SER682, ASP760
|
|
π-Cation: ARG555
π-Anion: ASP623, ASP760
|
83
|
Herbacetin
|
-7.7
|
ARG555, SER682, SER759
|
|
π-Cation: ARG555
π-Anion: ASP623
|
Trans Membrane Protease Serine 2 (TMPRSS2) receptor have been inhibited by 7 tested compounds, the binding score of these compounds revealed that they are most active as compared to standard drugs camostat (-7.6 kcal/mol) and nafamostat (-7.7 kcal/mol) as shown in Table 4. Silibinin with a binding affinity score of -8.4 kcal/mol showed signiifcant inhibitory activity against target receptor. The interactions between silibinin and residues in TMPRSS2 were shown in Fig. 6. This compound binds with VAL280, HIS296, THR393 and GLY439 residues of receptor through H-bond. Four amino acid residues VAL278, VAL280, CYS281 and CYS297 were formed π-alkyl interaction type, whereas two alkyl interactions were also established between VAL280 and LEU302 residues and Silibinin.
Table 4
Docking scores and interaction modes of active compounds with COVID-19 TMPRSS2 receptor
No.
|
Compound
|
Score
Kcal/mol
|
Interactions
|
H
|
hydrophobic
|
Other
|
59
|
Silibinin
|
-8.4
|
VAL280, HIS296, THR393, GLY439
|
π-Alkyl: VAL278, VAL280, CYS281, CYS297
Alkyl: VAL280, LEU302
|
|
60
|
Silymarin
|
-8.4
|
VAL280, HIS279, GLN438, GLY464
|
π-Sigma: VAL280
Amide-π Stacked: CYS437
π-Alkyl: CYS465
|
|
16
|
(18R)-Withaphysalin F
|
-8.4
|
HIS296, SER441
|
|
|
87
|
Tomentin A
|
-8.0
|
CYS297, SER436
|
π-Alkyl: CYS465
|
π-Cation: HIS296
|
89
|
Mimulone
|
-8.0
|
VAL280, GLN438, SER441
|
|
π-Cation: HIS296
π-Sulfur: CYS465
|
12
|
Limonin
|
-7.9
|
HIS296, GLY462, GLY464
|
Amide-π Stacked: TRP461
|
|
2
|
Withanolide
|
-7.9
|
HIS279, THR393
|
Alkyl: VAL280
|
|
The results of the binding score energy of the same selected compounds for Papain-like protease (PLpro) protein are presented in Table S1 and 195 compounds were scored higher than disulfiram (-4.1 kcal/mol). Therefore, somniferanolide appeared as the most active inhibitor for target protein. The binding site interaction between Somniferanolide and 6WX4 happens at TYR207 and VAL202 via π-alkyl and alkyl, respectively. Furthermore, four H-bonds were also observed between this compound and PLpro residues including ARG166, SER170, GLN174 and GLU203 (Fig. 7).
Table 5
Docking scores and interaction modes of top hit selected antiviral compounds with COVID-19 PLpro receptor
No.
|
Compound
|
Score
Kcal/mol
|
Interactions
|
H
|
hydrophobic
|
Other
|
10
|
Somniferanolide
|
-8.6
|
ARG166, SER170, GLN174, GLU203
|
π-Alkyl: TYR207
Alkyl: VAL202
|
|
16
|
(18R)-Withaphysalin F
|
-8.3
|
GLU161, TYR268
|
π-Sigma: TYR264
|
|
78
|
Inophyllum P
|
-8.2
|
TRP106, ALA288
|
π- π Stacked: TRP106
π-Alkyl: LYS105, ALA288
|
π-Cation: LYS105
|
1
|
Withanolide J
|
-8.2
|
ARG166, TYR207, LYS232
|
Alkyl: MET206
|
|
79
|
Inophyllum C
|
-8.1
|
ALA288, LEU289
|
π- π Stacked: TRP106
π-Alkyl: LYS105, ALA288, LEU289
|
π-Cation: LYS105
|
89
|
Mimulone
|
-8.1
|
GLU161, ASP164
|
π-Sigma: TYR264
π-Alkyl: TYR264
Alkyl: PRO247, PRO248
|
π-Cation: LYS157
|
14
|
Soulattrolide
|
-8.1
|
LYS157
|
π- π Stacked: TYR264
π-Alkyl: PRO248
|
|
77
|
Inophyllum B
|
-8.0
|
ARG166, SER170
|
π- π Stacked: TYR171
π-Alkyl: VAL202
|
|
16
|
Absinthin
|
-8.0
|
GLU167, TYR264
|
|
|
6
|
Withasomniferanolide
|
-7.9
|
LYS157, TYR268
|
Alkyl: LYS157
|
|
Worth mention, among docked compounds, withanolide, (18R)-withaphysalin F, limonin and mimulone showed the highest binding affinity as compared to the selected drugs in this study. Therefore, the stability of withanolide complexes with each Mpro and ACE receptors of COVID-19 were further analyzed using MD simulations for 100 ns duration.
In-silico ADMET prediction
In silico ADMET analysis of the target phytochemicals is presented in Table 6. The obtained results which were compared with the results as reported in literature [97], as well as referring to Lipinski's rule of five (LRO5) [98]. The first parameter observed was Mol Wt, where all compounds met the LRO5 criteria, which were ≤500 g/mol. Then compounds silibinin and silymarin showed TPSA values of more than 140 A², exceeding the criteria in LRO5 [99]. The number of hydrogen donor and acceptor groups of all compounds meets the LRO5 criteria, with no more than 5 and 10 groups, respectively. One compound exceeded the LRO5 criteria for the Log P value: compound mimulone with Log P of 5.75, although this value was still within the optimal range of Chander et al [100, 101]. All compounds showed ideal Log S values, although all of them were not in the optimal range of CaCo2 values. Only withasomnilide compounds were in the range of Log BB values from Chander et al [100, 101], although they all met the criteria for the number of Rot below 15. Finally, there are four compounds (3-HDH-withanolide F, somniferanolide, withasomnilide and (18r)-withaphysalin F) which were predicted to have high toxicity with an LD50 value of less than 50 mg/kg and were included in class II of the globally harmonized system of classification of labeling of chemicals (fatal if swallowed) [102], although all compounds did not showed potential as mutagens by the Ames test.
Further analysis of the target toxicity of each compound with ProTox-II was presented in Fig. 8. Some of the targets observed were potential carcinogenicity (carcino), immunotoxicity (immuno) and cytotoxicity (cyto); disrupting nuclear receptor signaling pathways on aryl hydrocarbon receptor (nr_ahr), androgen receptor ligand-binding domain (nr_ar_lbd), estrogen receptor alpha (nr_er) and estrogen receptor ligand-binding domain (nr_er_lbd); as well as interference with stress response pathways to mitochondrial membrane potential (sr_mmp) and tumor suppressor p53 (sr_p53). The results obtained were varied, with the most potential toxicity targets shown by aurantinidin with five targets in the probability range between 0.61 to 0.92. On the other hand, two compounds had the least potential toxicity targets (one target): mimulone and somniferanolide. However, the probability of somniferanolide (0.94) was much higher than mimulone (0.67). In addition, the pLD50 of somniferanolide (34 mg/kg) also much lesser than mimulone (2000 mg/kg), so the potential for toxicity of mimulone was lower than somniferanolide and other compounds. Meanwhile, the most reported target of toxicity was the immunotoxicity shown by nine compounds (except mimulone and aurantinidin). Overall, ADMET results exhibited that all tested compounds possessed the drug-likeness properties. However, some compounds such as limonin, withanolide, tomentin A and aurantinidin showed more ideal ADME properties, while the mimulone exhibited lower toxicity potential.
Table 6
The ADMET results of selected phytochemicals
Compounds
|
Mol Wta
|
TPSAb
|
HBD
|
HBA
|
log P
|
log Sc
|
Caco2d
|
log BB
|
Rot
|
Acute tox.e
|
pLD50f
|
Ames
|
Silymarin
|
482.441
|
155.14
|
5
|
10
|
2.3627
|
-3.204
|
0.435
|
-1.207
|
4
|
IV
|
2000
|
No
|
Limonin
|
470.518
|
104.57
|
0
|
8
|
3.1374
|
-4.379
|
0.952
|
-0.841
|
1
|
III
|
244
|
No
|
Withanolide
|
470.606
|
96.36
|
2
|
6
|
3.4954
|
-5.127
|
0.831
|
-0.315
|
2
|
III
|
300
|
No
|
Withasomnilide
|
470.606
|
96.36
|
2
|
6
|
3.4954
|
-5.209
|
0.834
|
0.022
|
2
|
II
|
7
|
No
|
(18R)-withaphysalin F
|
484.589
|
105.59
|
2
|
7
|
2.8318
|
-4.973
|
0.812
|
-0.422
|
1
|
II
|
7
|
No
|
Mimulone
|
408.494
|
86.99
|
3
|
5
|
5.7451
|
-3.817
|
0.651
|
-1.034
|
6
|
IV
|
2000
|
No
|
Silibinin
|
482.441
|
155.14
|
5
|
10
|
2.3627
|
-3.204
|
0.435
|
-1.207
|
4
|
IV
|
2000
|
No
|
Tomentin A
|
442.508
|
116.45
|
4
|
7
|
4.5348
|
-3.502
|
-0.366
|
-1.388
|
5
|
IV
|
2000
|
No
|
3-HDH-withanolide F
|
488.621
|
124.29
|
4
|
7
|
2.738
|
-5.209
|
0.668
|
-0.9
|
2
|
II
|
34
|
No
|
Somniferanolide
|
468.59
|
96.36
|
2
|
6
|
3.4155
|
-5.108
|
0.897
|
-0.347
|
2
|
II
|
34
|
No
|
Aurantinidin
|
287.247
|
114.29
|
5
|
6
|
2.9089
|
-2.968
|
-0.881
|
-1.612
|
1
|
V
|
3919
|
No
|
a Molecular weight in g/mol
b Topological polar surface area in Ų
c Aqueous solubility in log mol/L
d Predicted apparent Caco-2 cell permeability in 10-6 cm/s
e Acute toxicity according to the globally harmonized system (GHS) of classification and labeling of Chemicals
f Predicted LD50 in mg/kg
Molecular dynamics simulation
The best-docked pose with the highest binding affinity was utilized as the starting structure for the 100 ns Molecular dynamics simulation run. The best pose output is employed to build up this method in a high-throughput manner for studying the binding mechanism of the complex (ligand-protein) under clearly defined water environments. The different structures represented in Fig. 9 (a,b) give a visual impression of the sequence of events and the dynamics of the process, during the 1, 10, 20, 50, and 100 ns MD simulation runs. To examine structural stability, MD data is processed by calculating the RMSD. Complexes of 1O86 and 6LU7 with withanolide formed stable conformation after ~60 ns and ~85 ns, respectively, with an appropriate RMSD value of 2.65 and 2.91 Å, respectively, as seen in the RMSD plot (Fig. 10). The most acceptable RMSD value range is < 3.0 Å, as the lower RMSD value indicates superior stability of the system [102]. This finding shows that the withanolide develop a stable protein-ligand combination. The drop in the RMSD value of protein-ligand complexes reflects a conformational alteration in the protein secondary structure due to ligand binding. To examine the average flexibility and fluctuation of individual amino acids, the RMSF of protein and complexes was plotted from a 100 ns MD trajectory (Fig 11). Several times during the ligand-bound state, the fluctuations of the residues in the protein could be found in the RMSF plot. The average RMSF values of 1.481 and 1.936 Å was observed for complexes 1O86-withanolide and 6LU7-withanolide, respectively. The findings show that ligand-protein interaction brings protein chains closer and reduces the gap between them, represented in Fig 12. For evaluating and comparing protein structures, RR distance maps represent the average distance and standard deviation for all amino acid pairings between two conformations. The RR distance maps are represented in Fig. 13, which plots patterns of spatial interaction [103, 104]. The white diagonal on the map shows the zero distance between two residues, while the red and blue elements represent residue pairings with the greatest distance variances in the two conformations. The lowest radius of gyration (Rg) value of ~23.47 and ~21.04 Å was observed for 1O86-withanolide and 6LU7-withanolide complexes, respectively. Along the simulation time, Rg for the protein-ligand combination decreased, indicating that the structure became more compact (Fig. 14).
Grid-search on 17x19x18 and 14x18x16 grids, rcut = 0.35 for 1O86 and 6LU7, respectively, revealed the interactions of ligand with protein via hydrogen bonds, which were plotted against time (Fig. 15). On calculating hydrogen bonds between ligand (36 atoms) and proteins (6007 and 3038 atoms for receptors 1O86 and 6LU7, respectively), 835 donors and 1663 acceptors were found for 1O86-withanolide complex, while 435 donors and 853 acceptors were found for the 6LU7-withanolide complex. The average number of hydrogen bonds per timeframe was observed to be 0.662 and 0.775 out of 694302 and 185528 possible for 1O86-withanolide and 6LU7-withanolide complexes, respectively. It can conclude that the interactions between ligand-protein substantially enhanced the number of H-bonds. Fig. 16 shows the solvent accessible surface area (SASA) values changed due to the binding of the ligand to the protein. The reduced value of SASA of protein upon binding to ligand indicates the alteration of conformation in the protein structure and reduction in pocket size with increased hydrophobicity around it.