In silico study of medicinal plants with cyclodextrin inclusion complex as the potential inhibitors against SARS-CoV-2 main protease (Mpro) and spike (S) receptor.

The current outbreak of novel coronavirus disease (COVID-19) causes an alarming number of deaths in 221 countries around the world. Nowadays, there is no specific and effective drug regimen for curing COVID-19. Since the COVID-19 pandemic, several medicinal plants with promising results in the previous SARS-CoV could be used to treat SARS-CoV-2 infected patients. This work assesses proven medicinal plants as potential inhibitors against SARS-CoV-2 main protease (Mpro) and spike (S) receptors by employing in silico methods. Molecular docking studies and 3D structure-based pharmacophore modeling were performed to identify the molecular interactions of potential active molecules with the Mpro and (S) receptor of SARS-CoV-2. The drug-likeness and ADME properties were also predicted to support the drug-like nature of the selected active molecules. The results indicated that the most favorable ligand was Terrestriamide with (ΔG: ─8.70 kcal/mol; Ki: 0.417 μM) and (ΔG: ─7.02 kcal/mol; Ki: 7.21 μM) for Mpro and (S) receptor, respectively. Terrestriamide is also supported with a high drug-likeness value and appropriate ADME profile. Furthermore, to improve drug delivery, the cyclodextrin inclusion complex was calculated based on semi-empirical quantum mechanical methods. Terrestriamide/γ-cyclodextrin is the most favorable pathway of inclusion complex formation and could be used to treat COVID-19.

Introduction

An outbreak of novel coronavirus disease (n class="Disease">COVID-19) has been started since December 2019 in Wuhan, China. It has been declared as a global pandemic causing the death of 2,460,792 in 221 countries by February 22, 2021 [[1], [2], [3]]. The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is positive-sense, enveloped, and single-stranded RNA with a clinical manifestation closely resembling viral pneumonia [4]. SARS-CoV-2 is belonging to the Coronaviridae family and similar to SARS-CoV and MERS-CoV [5,6]. Coronaviruses represent the largest known genome of RNA viruses with more than six open reading frames (ORFs) [7]. Most recent studies have confirmed that the SARS-CoV-2 main protease (Mpro) and spike (S) receptors are promising drug targets [8,9]. The major ORF 1 ab is responsible to encode pp1a and pp1ab overlapping proteins. These proteins were then cleaved into 16 nonstructural proteins by the main protease (Mpro) [10,11]. This condition indicates that Mpro has an essential role in viral replication and transcription. The other genome encodes structural proteins including envelope protein (E), nucleocapsid phosphoprotein (N), and spike glycoprotein (S) [12]. Recently, SARS‐CoV‐2 has been confirmed utilizing the SARS‐CoV receptor spike (S) protein that composes a subunit containing a receptor-binding domain (RBD) [13]. This RBD engages with the host cell receptor angiotensin-converting enzyme (ACE2) [14]. Thus, inhibiting the function of Mpro and (S) protein would be useful to treat COVID19.

Several vaccines and drugs are being developed to prevent and cure the disease caused by SARS-CoV-2. However, there is no effective treatment that has yet been generated [15]. Most of the initial efforts to combat n class="Disease">COVID-19 are drug repurposing. Numerous FDA-approved drugs like chloroquine phosphate, hydroxychloroquine, umifenovir, remdesivir, ribavirin, and lopinavir have been used to treat infected patients [[16], [17], [18], [19]]. However, the clinical manifestation of these drugs against COVID-19 is still not fully understood [20]. Recently, secondary metabolites from medicinal plants have been used to combat COVID-19 with promising inhibitory effects against those of previous types of coronaviruses [21]. Several of them pose promising SARS-CoV protease inhibitory effects such as Tribulus terrestris, Psoralea corylifolia, and Angelica keiskei [[22], [23], [24]]. The potential inhibitory effect of active compounds derived from these medicinal plants could be predicted by employing in silico approach prior to an experimental effort.

In this work, we study computationally the molecular binding affinity and interactions between the chosen active molecules from medicinal plants with Mpro and (S) receptors as the drug target. To identify the potential inhibitors, literature searches were performed in various journals that focused on proven medicinal plants for treating n class="Species">coronaviruses. We employed the molecular docking simulation 3D structure-based pharmacophore to assess the molecular interaction between the ligand and the target receptor. The drug-likeness and ADME properties are predicted as further analysis to generate the best possible ligand. The best ligand was then further analyzed for protein structure flexibility compared to standard ligand-receptor complexes. Furthermore, the inclusion complex calculation based on semiempirical quantum mechanical methods is also performed to enhance drug delivery properties using a widely used drug carrier system, namely, cyclodextrins (CYD).

Computational methods

Data collection and construction of ligand structures

The X-ray crystal structures of two employed SARS-CoV-2 proteins, namely, n class="Disease">COVID-19 main protease in complex with an inhibitor N3 (PDB ID: 7BQY with 1.70 Å of resolution) [25] and SARS-CoV-2 spike receptor-binding domain bound with ACE2 (PDB ID: 6M0J with the resolution of 2.45 Å) [26] were retrieved from RCSB protein data bank (PDB). The medicinal plant's data was generated using a literature search focusing on proven inhibitory properties of secondary metabolites for combating coronaviruses. The literature search was conducted using PubMed, Scopus, WHO website, and Google Scholar. Scientific journals and case reports edited by WHO were included, while thesis and dissertations were not considered. The medicinal plant data were compiled in Table 1 . The molecular structures of the ligands were obtained from PubChem National Library of Medicine, National Center of Biotechnology Information (https://pubchem.ncbi.nlm.nih.gov/). The three-dimensional (3D) structures of ligands were optimized by Chem3D (PerkinElmer Inc.). The energy of ligands was also minimized using MM2 energy minimization in ChemDraw Professional 20.0 (PerkinElmer Inc.). The best ligand and cyclodextrins (CYDs) were then fully optimized by the semiempirical quantum mechanical PM6 and PM7 methods using Molecular Orbital PACkage (MOPAC) 2016 and Avogadro 1.2 software. The SARS-CoV-2 known spike receptor inhibitor (Arbidol) was retrieved from the literature [18].

Table 1

Selected medicinal plants with their active molecules.

No	Common name	Scientific name	Active molecules	Ref
1	Puncture vine	Tribulus terrestris L.	N-trans-Feruloyloctopamine, p-Coumaroyltyramine, n-Caffeoyltyramine, Terrestrimine, Terrestriamide	[22]
2	Babchi	Psoralea corylifolia	Bavachinin, Neobavaisoflavone, Isobavachalcone, 40-O-Methylbavachalcone, Psoralidin, Corylifol A	[23]
3	Ashitaba	Angelica keiskei	Isobavachalcone, 4-hydroxyderricin, Xanthoangelol, Xanthoangelol B-G	[24]
4	Red sage	Salvia miltiorrhiza	Tanshinone IIA, Tanshinone IIB, Methyl Tanshinonate, Cryptotanshinone, Tanshinone I, Dihydrotanshinone I, Rosmariquinone	[47]
5	Japanese alder	Alnus japonica	Platyphyllenone, Hirsutenone, Hirsutanonol, Oregonin, Rubranol, Rubranoside A, Rubranoside B, Amentoflavone, Herbacetin, Gallocatechin gallate, Pectolinarin, Rhoifolin, Emodin	[48]
6	Brown algae	Ecklonia cava	Eckol, Dioxinodehydroeckol, 2-phloroeckol, 7-phloroeckol, Fucodiphloroethol G, Dieckol, Phlorofucofuroeckol A	[49]
7	Japanese nutmeg-yew	Torreya nucifera L.	Bilobetin, Ginkgetin, Sciadopitysin, Apigenin, Luteolin, Quercetin, Amentoflavone	[50]
8	Chinese elder	Sambucus javanica subsp. chinensis (Lindl.) Fukuoka	Caffeic acid, Chlorogenic acid, Gallic acid	[51]

Molecular docking simulation and validation

The receptors and ligands were prepared for molecular docking simulation using AutoDockTools 1.5.6 [27]. The active site of the n class="Gene">Mpro receptor was defined by the redocking procedure of inhibitor N3 and yielded a Root Mean Square Deviation (RMSD) of 1.842, which is in the acceptable range. The active site of the (S) receptor was determined by utilizing CASTp 3.0 [28]. The receptor grid was determined by choosing the amino acid residues of the active sites (residue number 338, 339, 342, 343, 367, 368, 371, 373, 374) based on the CASTp 3.0 results. Both receptors and ligands were protonated. The Kollman charges were added to the receptor and the Gasteiger charges were added to the ligands [29]. The grid parameter file was based on the active site amino acid residue information from CASTp 3.0 that composed of 52 × 52 × 54 points with a space of 0.375 Å, then centered to the active site of the receptor (x = ─32.00; y = 11.00; z = 28.00). This grid parameter of the active site of spike proteins has been validated by a previous study [30]. Molecular docking simulation was performed by using AutoDock 4.2 (The Scripps Research Institute). The docking parameter file was based on the Lamarckian Genetic Algorithm (LGA) with 150 population sizes, 100 runs, 5,000,000 energy of evaluation, 0.8 rates of crossover, and 0.02 rate of gene mutation. The RMSD tolerance of 1.0 Å was employed to cluster the conformation results of the docking simulation [31]. The generated complexes of receptor and ligand were visualized using BIOVIA Discovery Studio Visualizer 2020 [32] and PyMOL 2.4 [33]. LigandScout Advanced 4.4 (Inte: Ligand GmbH, Vienna, Austria) was employed to define the features of ligand interaction for each pose within the active site of the receptor [34]. The docking simulation for the inclusion complex procedure between ligand and CYDs was also carried out by using AutoDock Tools 1.5.6. The guest and host structures were protonated. The grid boxes comprise 40 × 40 × 40 points with 0.375 Å space and were centered on the center site of host molecules. The box dimensions were (10.674 Å × 4.481 Å × 3.231 Å); (15.282 Å × 1.22 Å × 1.105 Å) and (−1.577 Å × 2.765 Å × 5.737 Å) for α−CYD, β−CYD and γ-CYD, respectively [35].

3D structure-based pharmacophore modeling

Pharmacophore is defined as an ensemble of steric that is important to ensure the optimal molecular interactions with a specific target and to trigger or inhibit its biological responses [36]. The 3D structure-based pharmacophore modeling was performed to assess the pharmacophore profile of ligands in the active pocket of the receptor. The validation of the interaction feature model was obtained based on a previous study [31]. The n class="Gene">ligands were screened by the validated 3D structure-based pharmacophore modeling method that was performed using LigandScout 4.4 Advanced algorithms [34].

ADME parameter prediction

The ligands with potency for inhibiting the receptors involved in viral replication were investigated for their bioavailability potency for oral administration and n class="Disease">toxicity properties using a free web tool to identify pharmacokinetics, namely SwissADME (http://www.swissadme.ch/) [37].

Toxicity and drug-likeness analysis

The ligaene">nds were further evaluated for various types of n class="Disease">toxicity properties such as tumorigenic, mutagenic, irritant, reproductive effectiveness, and drug similarity utilizing OSIRIS DataWarrior V5.2.1 [38].

Protein structure flexibility

To evaluate the effect of ligand interactions to the individual amino acid residues of n class="Species">SARS-CoV-2, Mpro, and (S) receptor, RMSF was calculated using CABS-flex 2.0 with 100 cycles [39]. The RMSF profile of the best ligand-receptor was also compared to the apo-protein (ligand-free protein), and standard ligand-receptor complexes.

Complexion energy calculation

The selected docked conformations of the best ligand/α, β, and γ−n class="Chemical">CYD inclusion complexes were then optimized geometrically by using the semi-empirical quantum mechanical PM6 and PM7 methods. The most stable conformation of the ligand-CYD inclusion complex was chosen based on the energy of complexation (ΔE) which is defined as the difference between the heat of the complex formation and the heat of formation of the free molecules involved which is represented by the formula [40].where E , E , and E represent the heat of formation of the inclusion complex, the isolated ligand molecule, and isolated cyclodextrin molecule, respectively.

Results

The entire ligands were docked into the most likely binding site of the receptor based on the defined coordinates. It is an important parameter to devote each ligand to having a binding site of greater affinity. Molecular docking simulation assessed the ligands until the convergence to the minimum energy was reached. The molecular docking studies on receptor-ligand interactions were focused on the amino acid residues of the active site. The most favorable conformation of each ligand was evaluated by calculating the affinity scoring function (free binding energy) (ΔG), inhibition constant (Ki), and hydrogen bonds as presented in Table 2 . The value of ΔG of ligand was ranged from −11.50 to −4.17 kcal/mol and −7.84 to −3.37 kcal/mol for ligand-Mpro and ligand-(S) receptor complexes, respectively. The value of Ki was ranged from 0.004 to 877.78 μM and 1.79 to 3400 μM for ligand-Mpro and ligand-(S) protein complexes, respectively. The most favorable ligand was terrestriamide with ΔG: ─ 7.02 kcal/mol; Ki: 7.21 μM and ΔG: ─ 8.70 kcal/mol; Ki: 0.417 μM for Mpro and (S) receptor. This result was obtained not only based on the lowest ΔG but also based on the ADME parameters, toxicity, and drug-likeness analysis. Furthermore, terrestriamide interacted with the key residues at the active site of the Mpro receptor, namely, Glu166, Gln189, and Thr190 that similar to inhibitor N3 as a standard inhibitor as can be seen in Fig. 1 . Terrestriamide also had hydrogen interactions with active residues of (S) receptors such as Phe342 and Asn343. The interactions of terrestriamide and arbidol (standard inhibitor) are depicted in Fig. 2 . These results indicated that these drugs could be the most favorable ligands in the ligand-binding domain (LBD) of the receptors.

Table 2

The docking simulation results of active molecules in the ligand-binding domain (LBD) of the target proteins.

Molecule name	Chemical structure	ΔG		Ki		Hydrogen bonds
		(kcal/mol)		(μM)
		6M0J	7BQY	6M0J	7BQY	6M0J	7BQY
Arbidol	Image 1	─6.16	N/A	30.74	N/A	Ser373	N/A
Inhibitor N3	Image 2	N/A	─6.00	N/A	40.26	N/A	Phe 140, Gly143, Cys145, His164, Glu166, Gln189, Thr190
N-trans-Feruloyloctopamine	Image 3	─6.03	─7.60	38.00	2.70	Cys336, Gly339, Asn343	Ser144, Glu166, Gln192
p-Coumaroyltyramine	Image 4	─6.50	─8.06	17.17	1.23	Ser371, Ser373	Tyr54, Gln192
n-Caffeoyltyramine	Image 5	─6.28	─7.51	25.04	3.13	Ala344, Ser373, Arg509	Ser144, His163, Gln192
Terrestrimine	Image 6	─5.95	─8.12	43.38	1.11	Asn343	Glu166, Gln189, Gln192
Terrestriamide	Image 7	─7.02	─8.70	7.21	0.417	Phe342,Asn343, Trp436	Tyr54, Glu166, Asp187, Gln189, Thr190
Bavachinin	Image 8	─7.14	─9.68	5.85	0.083	Cys336, Asn343	Gly143, Glu166, Asp187
Neobavaisoflavone	Image 9	─7.27	─9.11	4.67	0.211	Ser373	Thr190, Gln192
Isobavachalcone	Image 10	─6.69	─8.70	12.50	0.422	Cys336, Asn343	Glu166, Gln192
40′-O-Methylbavachalcone	Image 11	─7.00	─7.99	7.44	1.40	Val367	Tyr54, Glu166
Psoralidin	Image 12	─7.84	─9.71	1.79	0.077	Ser371, Ser373	His164, Thr190
Corylifol A	Image 13	─7.79	─9.40	1.96	0.128	─	Thr190, Gln192
4′-hydroxyderricin	Image 14	─6.31	─8.68	23.60	0.431	Asn343, Ser371	Glu166, Gln192
Xanthoangelol	Image 15	─6.77	─8.34	10.86	0.772	Trp436	Gln189, Thr190, Gln192
Xanthoangelol B	Image 16	─6.31	─7.89	23.54	1.66	Asp364, Ser371, Ser373	Gln192
Xanthoangelol C	Image 17	─6.45	─8.46	18.86	0.626	Cys336, Ser371	Gln189, Gln192
Xanthoangelol D	Image 18	─5.85	─7.66	51.92	2.41	Phe342, Val367, Asn440	Glu166, Gln192
Xanthoangelol E	Image 19	─6.12	─7.75	32.39	2.08	Cys336, Gly339	Asn142, Glu166, Gln192
Xanthoangelol F	Image 20	─7.12	─8.01	6.08	1.35	Asn440	Glu166, Gln192
Xanthoangelol G	Image 21	─6.20	─8.14	28.67	1.09	Asn440	His41, Glu166, Gln192
Tanshinone IIA	Image 22	─7.33	─8.73	4.22	0.397	Ser371	Cys145, Glu166
Tanshinone IIB	Image 23	─6.72	─8.49	11.87	0.60	Cys336, Gly339	Glu166
Methyl Tanshinonate	Image 24	─6.95	─8.89	8.04	0.303	Ser373	His163
Cryptotanshinone	Image 25	─7.47	─8.92	3.33	0.291	Asn343, Ser371	Glu166
Tanshinone I	Image 26	─6.95	─8.00	8.08	1.37	Ser373	Glu166
Dihydrotanshinone I	Image 27	─7.02	─8.31	7.15	0.815	Asn343, Ser371	Glu166
Rosmariquinone	Image 28	─7.38	─8.92	3.92	0.289	Ser371	Glu166
Platyphyllenone	Image 29	─6.84	─8.24	9.71	0.906	Asp364	Tyr54, Glu166, Gln192
Hirsutenone	Image 30	─6.54	─7.71	15.95	2.23	Phe342, Asp364	Ser144, Gln192
Hirsutanonol	Image 31	─5.86	─6.66	51.07	13.04	Cys336, Phe338, Gly339, Phe342, Asn343	Glu166, Gln189
Oregonin	Image 32	─4.72	─6.16	347.88	30.63	Asn343, Ala344, Arg509	Leu141, Asn142, Thr190
Rubranol	Image 33	─5.56	─6.35	84.27	22.22	Cys336, Gly339, Asn343	His163, Glu166
Rubranoside A	Image 34	─4.72	─7.26	348.39	4.79	Asp364, Ser371, Ser373	Ser144, His163, Gln192
Rubranoside B	Image 35	─4.98	─6.54	223.74	16.06	Asn343	His163, Gln192
Amentoflavone	Image 36	─7.56	─10.54	2.85	0.019	Asp364, Val367	Gly143, Ser144, Thr190
Herbacetin	Image 37	─5.98	─8.23	41.41	0.925	Asp364	Tyr54, Arg188, Gln189, Gln192
Gallocatechin gallate	Image 38	─5.80	─8.67	55.60	0.44	Cys336, Phe338, Asn343, Ser371	Tyr54, His163
Pectolinarin	Image 39	─4.86	─7.72	274.19	2.19	Asn343, Val367	Glu166, Thr190, Gln192
Rhoifolin	Image 40	─5.95	─7.39	43.37	3.83	Asn343, Ser373	His41, Gly143, Gln189, Gln192
Emodin	Image 41	─6.26	─7.24	25.61	4.90	Ser371, Ser373	Gln189, Gln192
Eckol	Image 42	─5.78	─7.90	57.92	1.62	Phe342, Asn343, Arg509	Ser144, His163, Arg188
Dioxinodehydroeckol	Image 43	─6.14	─7.85	31.79	1.75	Asp364	Asn142, His163, His164, Asp187
2′-phloroeckol	Image 44	─6.38	─8.48	20.96	0.603	Cys336, Ser371	Leu141, Ser144, His163, Glu166
7′-phloroeckol	Image 45	─5.94	─7.47	44.51	3.37	Asn343, Ser371, Arg509	Asn142, Ser144, Glu166, Gln192
Fucodiphloroethol G	Image 46	─4.48	─6.65	516.29	13.45	Cys336, Gly339, Glu 340, Asp364, Val367	His163, His164, Glu166
Dieckol	Image 47	─6.30	─6.95	23.98	8.05	Asn343, Ser373	Asn142, Gly143, Glu166, Gln189
Phlorofucofuroeckol A	Image 48	─7.06	─8.15	6.63	1.06	Asn343, Asp364	His41, Thr190
Bilobetin	Image 49	─6.74	─11.05	11.49	0.008	Ser373, Arg509	Gly143, His163, Thr190
Ginkgetin	Image 50	─6.73	─11.50	11.63	0.004	Asn343, Ser371	Gly143, Ser144, His163, His164, Thr190
Sciadopitysin	Image 51	─6.71	─11.24	12.12	0.006	Asn343, Arg509	Gly143, His163, Thr190
Apigenin	Image 52	─6.30	─8.40	23.96	0.698	Phe342, Ser371, Ser373	His164, Gln192
Luteolin	Image 53	─6.40	─8.20	20.24	0.979	Asp364, Val367	His164, Glu166, Gln192
Quercetin	Image 54	─6.28	─8.20	25.09	0.979	Cys336, Val367	His164, Glu166, Asp187, Thr190
Caffeic acid	Image 55	─4.71	─4.91	352.47	253.79	Ser371, Ser373	Glu166, Gln192
Chlorogenic acid	Image 56	─5.59	─6.64	79.34	13.62	Cys336, Asn343, Asp364, Ser373	Gly143, His163, Thr190
Gallic acid	Image 57	─3.37	─4.17	3400	877.78	Cys336, Phe338, Gly339, Asp364	Glu166, Thr190, Gln192

Fig. 1

Ligand interaction diagram of (a) inhibitor N3 and (b) terrestriamide in Ligand Binding Pocket (LBD) of SARS-CoV-2 Mpro (PDB ID: 7BQY).

Fig. 2

Ligand interaction diagram of (a) arbidol and (b) terrestriamide in Ligand Binding Pocket (LBD) of SARS-CoV-2 (S) receptor (PDB ID: 6M0J).

Ligaene">nd interaction diagram of (a) inhibitor n class="Gene">N3 and (b) terrestriamide in Ligand Binding Pocket (LBD) of SARS-CoV-2 Mpro (PDB ID: 7BQY).

Ligaene">nd interaction diagram of (a) n class="Chemical">arbidol and (b) terrestriamide in Ligand Binding Pocket (LBD) of SARS-CoV-2 (S) receptor (PDB ID: 6M0J).

The structure-based best-docked ligand conformation of n class="Chemical">terrestriamide was evaluated to investigate the specific molecular interaction in the active pocket of the receptors. The molecular properties including hydrophobic, hydrogen bond donor and hydrogen acceptor are presented as yellow spheres, green arrows, and red arrows, respectively, as shown in Fig. 3 and Fig. 4 for Mpro and (S) receptor, respectively. Pharmacophore modeling was used to investigate the crucial parts of each ligand that affect the molecular behavior of the receptors. This modeling was according to the best-docked ligand conformations. Terrestriamide has shown hydrophobic interactions with the active site of the receptor. This phenomenon happened due to the LBD of the receptor is predominantly formed as the hydrophobic cavity composed of amino acid residues. Terrestriamide as the most promising ligand has shown the hydrophobic interaction between the benzene ring with residues of Met 165 and Leu 167 in the Mpro complex. The hydrogen bond acceptor (HBA) occurred between the carboxyl group of residues Glu166 and Gln192. On other hand, hydrogen bond donor (HBD) happened between the hydroxyl group of Tyr54, Asp187, and Thr190 (Fig. 3). Furthermore, terrestriamide has shown hydrophobic interactions between benzene rings of Leu 441, Phe342, Leu 368, and Phe 374. The prop-2-enamide group indicated HBA with Trp436 and HBD with Asn343 (Fig. 4). This result denoted that the pharmacophore of this ligand was predominantly derived from the hydrophobic and hydrogen interactions in the cavity of the receptor and showed a very high agreement with the docking simulation result, especially for the key residue interaction.

Fig. 3

(a) 3D and (b) 2D structure-based pharmacophore modeling of the best-docked pose of terrestriamide in SARS-CoV-2 Mpro (PDB ID: 7BQY). Hydrophobic, hydrogen bond donor, and hydrogen bond acceptor interactions are represented as yellow spheres, green, and red arrows (spheres), respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 4

3D and (b) 2D structure-based pharmacophore modeling of the best-docked pose of terrestriamide in SARS-CoV-2 (S) receptor (PDB ID: 6M0J). Hydrophobic, hydrogen bond donor, and hydrogen bond acceptor interactions are represented as yellow spheres, green, and red arrows (spheres), respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

(a) 3D and (b) 2D structure-based pharmacophore modeling of the best-docked pose of terrestriamide iene">n n class="Species">SARS-CoV-2 Mpro (PDB ID: 7BQY). Hydrophobic, hydrogen bond donor, and hydrogen bond acceptor interactions are represented as yellow spheres, green, and red arrows (spheres), respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

3D and (b) 2D structure-based pharmacophore modeling of the best-docked pose of terrestriamide iene">n n class="Species">SARS-CoV-2 (S) receptor (PDB ID: 6M0J). Hydrophobic, hydrogen bond donor, and hydrogen bond acceptor interactions are represented as yellow spheres, green, and red arrows (spheres), respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

The prediction of ADME parameters of ligands was performed for ligands that have lower values of ΔG than inhibitor and/or arbidol as standard inhibitors as can be seen in Table 3 . This prediction is based on Lipinski's rule of Five regarding the active entity administered orally with four physicochemical parameters (molar mass ≤ 500, logP ≤ 5, hydrogen bond donor ≤ 5, and hydrogen bond acceptors ≤ 10) correlated with 90% of the drug administered orally that has reached phase II of a clinical trial [41]. The Topological Polar Surface Area (TPSA), Csp 3 fraction, and cytochrome (CYP) inhibitor were also evaluated. ADME prediction was performed to evaluate the chemical properties and bioavailability potential of ligands for oral administration. Among the ligands, amentoflavone, dieckol, phlorofucofuroeckol A, bilobetin, ginkgetin, and sciadopitysin showed a molecular mass above the recommended values (500 g/mol). The other ligands have appropriate properties for oral administration based on Lipinski's rule. For TPSA, the degree of polarity was also predicted. The value of TPSA was ranged from 34.14 to 287.14 Å2. The higher the value of TPSA, then the higher possibility to hinder the metabolic biotransformation of the structure and subsequent absorption. Terrestriamide showed a TPSA value of 95.86 Å2. Thus, this ligand has good biotransformation during oral absorption due to the ADME guideline requirement (TPSA < 140 Å2, good intestinal absorption) [42]. Most of the ligands also interact with CYP1A2 and CYP3A4, presenting a high gastrointestinal absorption [43].

Table 3

ADME parameters prediction results of selected active molecules. (MM: molecular weight; HBD: hydrogen bond donor; HBA: hydrogen bond acceptor, TPSA: topological polar surface area).

Molecule name	MM (g/mol)	x log P3	HBD	HBA	TPSA (Å2)	CYP inhibitor
p-Coumaroyltyramine	283.32	2.72	3	3	69.56	CYP2D6, CYP3A4
n-Caffeoyltyramine	299.32	2.36	4	4	89.79	CYP2C9, CYP2D6, CYP3A4
Terrestriamide	327.33	2.32	3	5	95.86	CYP2C9, CYP3A4
Bavachinin	338.40	4.45	1	4	55.76	CYP2C19, CYP2C9, CYP2D6, CYP3A4
Neobavaisoflavone	322.35	4.40	2	4	70.67	CYP1A2, CYP2C9, CYP3A4
Isobavachalcone	324.37	5.10	3	4	77.76	CYP1A2, CYP2C9, CYP3A4
40′-O-Methylbavachalcone	352.42	5.76	1	4	55.76	CYP1A2, CYP2C19, CYP2C9, CYP3A4
Psoralidin	336.34	4.69	2	5	83.81	CYP1A2, CYP2C19, CYP2C9
Corylifol A	390.47	6.25	2	4	70.67	CYP2C19, CYP3A4
4′-hydroxyderricin	338.40	5.43	2	4	66.76	CYP1A2, CYP2C19, CYP2C9, CYP3A4
Xanthoangelol	392.49	6.96	3	4	77.76	CYP1A2, CYP2C9, CYP3A4
Xanthoangelol B	408.49	5.97	4	5	97.99	CYP1A2, CYP2C9, CYP3A4
Xanthoangelol C	366.41	4.43	3	5	94.83	CYP1A2, CYP2C9, CYP3A4
Xanthoangelol F	406.51	7.29	2	4	66.76	CYP1A2, CYP2C9, CYP3A4
Xanthoangelol G	422.51	6.30	3	5	86.99	CYP1A2, CYP2C9, CYP3A4
Tanshinone IIA	294.34	4.33	0	3	47.28	CYP1A2, CYP2C19, CYP2C9, CYP2D6, CYP3A4
Tanshinone IIB	310.34	2.93	1	4	67.51	CYP1A2, CYP2C19, CYP2C9, CYP2D6, CYP3A4
Methyl Tanshinonate	338.35	3.20	0	5	73.58	CYP1A2, CYP2C19, CYP2C9, CYP3A4
Cryptotanshinone	296.36	3.80	0	3	43.37	CYP1A2, CYP2C19, CYP2C9, CYP3A4
Tanshinone I	276.29	3.69	0	3	47.28	CYP1A2, CYP2C19, CYP3A4
Dihydrotanshinone I	278.30	3.16	0	3	43.37	CYP1A2, CYP2C19, CYP2C9, CYP2D6, CYP3A4
Rosmariquinone	282.38	4.88	0	2	34.14	CYP2C9
Platyphyllenone	296.36	3.80	2	3	57.53	CYP1A2, CYP2C19, CYP2D6, CYP3A4
Hirsutenone	328.36	3.09	4	5	97.99	CYP1A2, CYP2C9, CYP2D6, CYP3A4
Amentoflavone	538.46	5.04	6	10	181.80	─
2′-phloroeckol	496.38	3.27	8	12	198.76	CYP2C9
Dieckol	742.55	4.87	11	18	287.14	CYP2C9
Phlorofucofuroeckol A	602.46	4.66	9	14	232.13	CYP2C9
Bilobetin	552.48	5.36	5	10	170.80	CYP2C9
Ginkgetin	566.51	5.69	4	10	159.80	CYP2C9
Sciadopitysin	580.54	6.02	3	10	148.80	─
Apigenin	270.24	3.02	3	5	90.90	CYP1A2, CYP2D6, CYP3A4
Luteolin	286.24	2.53	4	6	111.13	CYP1A2, CYP2D6, CYP3A4
Quercetin	302.24	1.54	5	7	131.36	CYP1A2, CYP2D6, CYP3A4

ADME paran class="Chemical">meters prediction results of selected active molecules. (MM: molecular weight; HBD: hydrogen bond donor; HBA: hydrogen bond acceptor, TPSA: topological polar surface area).

The results of the toxicity and drug-likeness property analysis are shown iene">n Table 4 for selected n class="Gene">ligands. Among all screened ligands, terrestriamide possessed the highest drug likeliness property, followed by bilobetin, ginkgetin, sciadopitysin, apigenin, and luteolin. Based on the training dataset employed by OSIRIS DataWarrior V5.2.1, the ligands with a higher or positive value of drug-likeliness are considered as good drug candidates. Terrestriamide as the best docking conformation, non-mutagenic, and non-tumorigenic ligand was further subjected to protein flexibility and inclusion complex studies.

Table 4

Toxicity and drug-likeness analysis of selected active molecules.

Molecule Name	Drug likeness	Mutagenic	Tumorigenic	ReproductiveEffective	Irritant
p-Coumaroyltyramine	0.26453	None	None	None	None
n-Caffeoyltyramine	0.26453	None	None	None	None
Terrestriamide	1.05980	None	None	None	None
Bavachinin	─0.16737	None	None	None	None
Neobavaisoflavone	─0.41739	None	None	None	None
Isobavachalcone	─0.47336	None	High	None	High
40′-O-Methylbavachalcone	─0.26787	None	None	Low	High
Psoralidin	─0.53359	None	None	High	None
Corylifol A	─1.91870	None	None	None	None
4′-hydroxyderricin	─0.26787	None	None	None	None
Xanthoangelol	─1.86470	None	High	None	High
Xanthoangelol B	─13.10500	None	High	None	High
Xanthoangelol C	─2.77470	High	High	None	High
Xanthoangelol F	─1.65310	None	None	None	None
Xanthoangelol G	─12.94700	None	None	None	None
Tanshinone IIA	─7.78620	None	None	High	None
Tanshinone IIB	─12.76000	None	None	High	None
Methyl Tanshinonate	─9.93700	None	None	High	None
Cryptotanshinone	─7.22240	None	None	High	None
Tanshinone I	─3.70550	None	High	High	High
Dihydrotanshinone I	─3.11350	None	High	High	High
Rosmariquinone	─7.77380	None	None	High	None
Platyphyllenone	─4.92300	None	None	None	None
Hirsutenone	─4.92300	None	None	None	None
Amentoflavone	0.28194	None	None	None	None
2′-phloroeckol	─2.20540	Low	None	None	None
Dieckol	─2.20540	Low	None	None	None
Phlorofucofuroeckol A	─2.06380	Low	None	None	None
Bilobetin	0.40331	None	None	None	None
Ginkgetin	0.40331	None	None	None	None
Sciadopitysin	0.40331	None	None	None	None
Apigenin	0.28194	High	None	None	None
Luteolin	0.28194	None	None	None	None
Quercetin	─0.08283	High	High	None	None

Toxicity aene">nd drug-likeene">ness aene">nalysis of selected active molecules.

To assess the flexibility profile of individual amino acid residues of SARS-CoV-2 n class="Gene">Mpro and (S) receptor, root means square fluctuation (RMSF) was determined using CABS-flex 2.0 based on the protein dynamics. The dynamics structures of a protein represent its biological functions. The CABS models utilize the asymmetric Metropolis and Monte Carlo dynamics scheme, which meets the requirements of microscopic reversibility [39]. The Mpro (7BQY)-inhibitor and Mpro (7BQY)-terrestriamide complexes formed the disruption of hydrogen bonds with Glu166, Gln189, and Thr190 residues thus showed a higher fluctuation compared to apo-protein of Mpro (7BQY) as can be seen in Fig. 5 . The increase of fluctuation is also shown at key residues of (S) receptor (6M0J) caused by arbidol and terrestriamide as the inhibitors. Residue numbers 368–373 have a fluctuation movement while the apo-protein structure was more stable as depicted in Fig. 6 . Thus, terrestriamide has the potential to be an inhibitor due to the similarity with the standard inhibitors (inhibitor N3 and arbidol) with a higher value of flexibility.

Fig. 5

RMSF profile of apo-protein and Mpro-ligand (inhibitor N3 and terrestriamide) complexes.

Fig. 6

RMSF profile of apo-protein and (S) receptor-ligand (arbidol and terrestriamide) complexes.

RMSF profile of apo-protein and Mpro-n class="Gene">ligand (inhibitor N3 and terrestriamide) complexes.

RMSF profile of apo-protein and (S) receptor-ligaene">nd (n class="Chemical">arbidol and terrestriamide) complexes.

To optimize the drug delivery system of terrestriamide, we employed the iene">nclusion complex with n class="Chemical">cyclodextrins (CYDs). Cyclodextrins (CYDs) are cyclic oligosaccharides containing a hydrophilic shell and a hydrophobic core with a hydroxyl group coated to the outside and a glucose residue connected to the inside structure (Fig. 7 ). CYD improves stability, solubility, and bioavailability according to low immunogenicity and nontoxic properties [44,45]. The inclusion complexes between terrestriamide with α-CYD, β-CYD, and γ-CYD obtained from PM6 and PM7 calculations were stabilized in a water environment. The free energy values of terrestriamide/CYDs from the PM6 and PM7 methods were slightly similar that ranged from ─5.35 to ─4.33 kcal/mol and ─5.32 to ─4.30 kcal/mol for PM6 and PM7, respectively (Table 5 and Table 6 ). The inclusion complex calculation yielded terrestriamide/γ-CYD as the most favorable inclusion complex due to the lowest value of complexation energy (ΔE) with −197.68 and −269.37 kcal/mol for PM6 and PM7 calculation methods, respectively (see Table 7).

Fig. 7

Top prospective of the basic structure of cyclodextrin.

Table 5

The molecular docking results at 298.15 K. The geometrical structures of terrestriamide and CYDs were performed by the semiempirical quantum mechanical method of PM6.

No	Guest/Host	RMSD	Cluster	ΔG (kcal/mol)
				Lowest	Average
1	Terrestriamide/α−CYD	1.33	9	−4.94	−4.33
2	Terrestriamide/β−CYD	1.74	11	−4.86	−4.48
3	Terrestriamide/γ−CYD	1.33	26	−5.68	−5.35

Table 6

The molecular docking results at 298.15 K. The geometrical structures of terrestriamide and CYDs were performed by the semiempirical quantum mechanical method of PM7.

No	Guest/Host	RMSD	Cluster	ΔG (kcal/mol)
				Lowest	Average
1	Terrestriamide/α−CYD	1.57	17	−5.10	−4.47
2	Terrestriamide/β−CYD	1.15	16	−5.02	−4.30
3	Terrestriamide/γ−CYD	0.36	24	−5.79	−5.32

Table 7

Heat of formation energy (E) and the complexation energy (ΔE) of the minimized inclusion complexes based on semiempirical quantum mechanical PM6 and PM7 methods.

No	Molecule	PM6		PM7
		E (kcal/mol)	ΔE (kcal/mol)	E (kcal/mol)	ΔE (kcal/mol)
	Isolated molecule
1	Terrestriamide	−138.67		−141.33
2	α−CYD	−1340.77		−1342.04
3	β−CYD	−1550.12		−1555.74
4	γ−CYD	−1624.38		−1628.26
	Inclusion complex
56	Terrestriamide/α−CYD	−1676.52	−197.08	−1679.48	−196.11
	Terrestriamide/β−CYD	−1786.14	−97.35	−1792.35	−95.28
7	Terrestriamide/γ−CYD	−1960.73	−197.68	−2038.96	−269.37

The molecular docking results at 298.15 K. The geometrical structures of n class="Chemical">terrestriamide and CYDs were performed by the semiempirical quantum mechanical method of PM6.

The molecular docking results at 298.15 K. The geometrical structures of n class="Chemical">terrestriamide and CYDs were performed by the semiempirical quantum mechanical method of PM7.

The guest molecule (terrestriamide) was located sn class="Gene">lightly centered to the γ−CYD as the host molecule in all complex conformations, as shown in Fig. 8 a and Fig. 8b based on the existence of methyl groups at the primary hydroxyl group of all glucose units (C2 and C3 positions) of CYD. After the insertion process of terrestriamide into the cavity of γ−CYD, the methoxyl groups at the C2 and C3 position of γ−CD move deeper to the cavity according to the existence of carboxyl groups of terrestriamide. Thus, the terrestriamide moved deeper inside the cavity. The generated terrestriamide/γ−CYD obtained from PM7 calculation is a preferable complex with a lower ΔE compared to the inclusion complex generated from the PM6 calculation. The carboxyl and hydroxyl groups fall into the cavity from PM7 calculation and pushed the terrestriamide deeper inside the cavity due to the presence of the hydrophobic interaction. Terrestriamide/γ−CYD inclusion complex formed hydrogen bonds. The hydrogen bonds occurred between the hydroxyl and carboxyl groups of terrestriamide and the ether-like anomeric oxygen atom of γ−CYD. The hydrogen bonds maintain the terrestriamide tightly as the hydrophobic atoms of γ−CYD enclosing. The hydrogen bond interactions and their distance of terrestriamide/γ−CYD inclusion complexes in an aqueous environment are depicted in Table 8 .

Fig. 8

Hydrogen bond distance in 1:1 terrestriamide/γ−CD inclusion complex. (a) terrestriamide/γ−CD generated from semi-empirical quantum mechanical PM6 method, (b) terrestriamide/γ−CD generated from semi-empirical quantum mechanical PM7 method.

Table 8

The distance of hydrogen bonds between terrestriamide as guest and γ−CYD as the host obtained from PM6 and PM7 inclusion complexes.

No	Method	Hydrogen Bond	Distance (Å)
1	PM6	O(γ−CYD) … H(OH−TER)	2.07
		O(O−TER) … H(γ−CYD)	2.15
		O(O = TER) … H(γ−CYD)	2.92
		O(O−TER) … H(γ−CYD)	3.01
2	PM7	O(γ−CYD) … H(OH−TER)	1.93
		O(O = TER) … H(γ−CYD)	2.98
		O(O = TER) … H(γ−CYD)	3.28

Hydrogen bond distance iene">n 1:1 n class="Chemical">terrestriamide/γ−CD inclusion complex. (a) terrestriamide/γ−CD generated from semi-empirical quantum mechanical PM6 method, (b) terrestriamide/γ−CD generated from semi-empirical quantum mechanical PM7 method.

The distance of hydrogen boene">nds betweeene">n n class="Chemical">terrestriamide as guest and γ−CYD as the host obtained from PM6 and PM7 inclusion complexes.

Discussion

This study aimed to assess the medicinal plants in combating SARS-CoV maiene">n protease (n class="Gene">Mpro) and spike (S) receptors. The medicinal plants were chosen based on their proven inhibitory activity against previous SARS-CoV due to the whole genome of SARS-CoV-2 shares about 80% identity with SARS-CoV [46]. Tribulus terrestris L. as the best selected medicinal plant showed promising inhibitory effects for SARS-CoV Mpro with IC50 value ranged from 15.8 to 70.1 μM. The inhibitory activity was associated with the presence of polar substituent including alcohol and ketone on the methylene groups (C7′ and C8′) [22]. Based on these findings, we assessed the inhibitory activity of main bioactive molecules in Tribulus terrestris L. in the SARS-CoV-2 Mpro and (S) protein using in silico methods. The molecular docking simulation yielded terrestriamide ((E)-3-(4-hydroxy-3-methoxyphenyl)-N-[2-(4-hydroxyphenyl)-2-oxoethyl]prop-2-enamide) as the best bioactive constituent due to the lowest value of free binding energy and inhibitory constant. Figs. 1 and 2 showed that terrestriamide generated five and three hydrogen bonds with Mpro and (S) protein, respectively. The hydrogen bonds were occurred predominantly within hydroxyl and ketone groups. These hydrogen bonds stabilize the ligand-protein interactions. The pharmacophore modeling also supports the results by presenting the HBD and HBA within the hydroxyl and ketone groups in terrestriamide, while the benzene groups pose the hydrophobic interactions. These results indicate that the hydroxyl and ketone groups act as the essential parts and could generate the biological activity in the target receptors and meet agreement with previous studies in SARS-CoV. To identify the antagonist effect of terrestriamide in the receptors, flexibility profile of individual amino acid residues of SARS-CoV-2 Mpro and (S) receptor was determined. Terrestriamide poses a higher value of flexibility compared to apo-protein and also compared to inhibitor N3 and arbidol as the standard inhibitor agent for Mpro and (S) protein, respectively. This result indicates that terrestriamide has the better effect in disrupting the stability of receptors. Besides the molecular interactions, terrestriamide also supported by the appropriate properties for oral administration based on the Lipinski's rule of five and ADMET predictions. To optimize the delivery properties, the best bioactive molecule was complexed into the cyclodextrin (CYD) as the carrier. The inclusion complex was favorable for terrestriamide with γ-CYD based on the highest negative value of free binding energy (ΔG) and complexation energy (ΔE). Similar to the docking simulation and pharmacophore modeling, hydroxyl groups were also plays crucial role in binding with the CYD and keep the terrestriamide within the cavity of CYD. Thus, this bioactive molecule could be used for further analysis such as in vitro and in vivo tests.

Conclusions

In this study, secondary metabolites from mediciene">nal plants that have been reported to iene">nhibit other n class="Species">coronaviruses could act as a potential inhibitor to SARS-CoV-2 main protease (Mpro) and spike (S) receptors. From the screening, terrestriamide appears to be the best ligand having high potency against the target receptors via the lowest free binding energy and meets the ADME and drug-likeness requirements. The inhibitory properties of terrestriamide are also supported by the protein flexibility profile compared to the standard inhibitors. The drug delivery system of terrestriamide could be facilitated by the inclusion complex of γ−cyclodextrin (γ−CYD) due to the lowest value of complexation energy (ΔE). Based on this finding, we can conclude that terrestriamide has the potential to inhibit SARS-CoV-2 Mpro and S receptors. Hence, this molecule is worth being proposed for further in vitro and in vivo studies against Mpro and S receptors of SARS-CoV-2 to validate the results.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.