Identification of bioactive molecules from Triphala (Ayurvedic herbal formulation) as potential inhibitors of SARS-CoV-2 main protease (Mpro) through computational investigations

2.1 Protein preparation

The crystal structure of the Mpro of SARS-CoV-2 bound with non-covalent inhibitor X77 (N-(4-tert-butylphenyl)-N-[(1R)-2-(cyclohexylamino)-2-oxo-1-(pyridin-3-yl) ethyl]-1H-imidazole-4-carboxamide) at a resolution of 2.1 Å (PDB id: 6W63) (Rudrapal et al., 2021b) was retrieved from RCSB Protein data bank (https://www.rcsb.org/). The maestro's protein preparation wizard (version 2021-2 Schrodinger suite) was used to prepare proteins by removing all water within 5 Å distance from the ligand, filling amino acid side chains, adding missing hydrogen atoms, creating disulfide bonds, and optimizing hydrogen bonds. For restrained minimization, OPLS4 force field was used until the average RMSD (root-mean-square deviation) of the non-hydrogen atom converged to 0.30 Å, allowing heavy atoms to move sufficiently to release strained bonds, angles, and clashes (Roos et al., 2019; Kumar et al., 2021).

2.2 Ligand preparation

The Triphala formulation's bioactive molecules were chosen from the literature (Prasad and Srivastava, 2020), and their chemical structures were obtained from a popular online database, PubChem (https://pubchem.ncbi.nlm.nih.gov/). LigPrep was used to generate the ligands' three-dimensional (3D) coordinates (Chen and Foloppe, 2010). Ionization state (pH 7.4), tautomer, and chirality were predicted by Epik module. The OPLS4 force field was employed for energy minimization until an energetically stable confirmation was obtained (Roos et al., 2019).

2.3 Molecular docking study

The selected set of phytoconstituents was subjected for molecular docking against Mpro with the Extra Precision (XP) protocol of the glide docking module of Maestro (Friesner et al., 2006). The glide-receptor grid generation module was used with the default settings to identify the binding site, and the centroid of the bound ligand was taken into account for grid generation (Halgren et al., 2004). The single best pose for each molecule was saved and the compounds were then classified based on their docking scores, which represent free binding energies (Othman et al., 2021).

2.4 Molecular dynamics simulations

Following that, the top four compounds with the best binding free energies, namely terflavin A, corilagin, chebulagic acid, and chebulinic acid, were selected to further ensure the flexibility and stability at the binding site of M^pro. A 100 ns MD simulation was performed using GROnin-gen MAchine for Chemical Simulation (GROMACS) 2019.2 software (Abraham et al., 2015). The topologies of the protein and ligand were determined using the ‘pdb2gmx’ script and the GlycoBioChem PRODRG2 server (http://davapc1.bioch.dundee.ac.uk/cgi-bin/prodrg), respectively. Energy minimization was accomplished using the GROMOS96 54a7 force field. The complex was solvated using the SCP water model, electrically neutralized by adding sodium ions and chloride ions using ‘gmxgenion’ script, equilibrated under an 0.3 ns isothermal-isochoric (NVT) and 0.3 ns isothermal-isobaric (NPT) ensembles and subjected to MD simulations (100 ns) (Ghosh et al., 2021). The pressure and temperature of the system were maintained at 1 atm and 300 K, respectively. The leap frog algorithm was used to integrate the motion equation in triplicate at a time step of 2 fs. Finally, the protein–ligand complexes were put to generate trajectories of RMSD, RMSF, Rg and H bond by GROMACS scripts. The trajectories were analyzed and visualized by VMD 1.9.3 (Humphrey et al., 1996).

2.5 MM-PBSA calculations

The binding free energies of protein–ligand complexes were quantified by the Molecular Mechanics Poisson Boltzmann Surface Area (MM-PBSA). This technique allows us to study the biomolecular interactions between the target and bioactive molecules. The g_mmpbsa script was utilized for MM-PBSA binding free energy calculations (Kumari et al., 2014). Van der Waals, polar solvation, and electrostatic contributions to energy were calculated. Solvent accessible surface area (SASA) was used to calculate non-polar energy contribution.

2.6 Drug-likeness, ADMET and toxicity prediction

The selected bioactive molecules were further checked for drug-likeness properties using Lipinski’s rule of five (RO5) (Hussain et al., 2021; Rudrapal and Sowmya, 2019). The ADMET properties including absorption, distribution, metabolism, excretion and toxicity of the selected were determined by the admetSAR server (Cheng et al., 2012; Ghost et al., 2021). Oral in silico toxicity prediction was performed using web-based server ProTox-II (Banerjee et al., 2018).

3.1 Molecular docking

Triphala (Ayurvedic herbal formulation) is made up of dried fruits from three different herbal plants in equal parts (Jantrapirom et al., 2021; Tarasiuk et al., 2018). Following a thorough review of the literature on the three plants T. chebula, T. bellirica, and P. emblica, 19 bioactive molecules were chosen for molecular docking against the M^pro of SARS-CoV-2. The binding pocket around the inhibitor X77 of the target protein M^pro was used as a reference for determining the grid box. Docking analysis generated many configurations, which were scored to determine the most favorable binding modes. The bioactive molecules docked against M^pro of SARS-CoV-2 exhibited excellent binding, as evidenced by high docking scores that exceeded X77. The validation of the docking method was carried out by redocking of co-crystal ligand, X77 in the original ligand binding site of Mpro. The redocking was successful. The measured RMSD value of redocking was 0.4493 Å between the experimental docking pose of X77 and the natural binding mode of the native ligand. The superimposed conformation of the experimental ligand with the original co-crytal ligand, X77 is displayed in Fig. 2. The docking results and the binding free energies of 19 compounds are given in Table 1.

Terflavin A achieved the highest docking score, with a value of −14,521 kcal/mol. Corilagin, chebulagic acid, chebulinic acid and X77 had XP Gscore values of −12.845 kcal/mol, −11,359 kcal/mol, −11,217 kcal/mol, and −5,966 kcal/mol, respectively. Chemical structures of the four best-docked compounds and the native ligand X77 are depicted in Fig. 1. Chebulagic acid has a lower XP Gscore value than chebulinic acid, but the binding energy value of chebulinic acid is better (−64.05 kcal/mol). Furthermore, despite having a higher XP Gscore, the dG bind value of the X77 ligand is −85,045 kcal/mol, lower than corilagin, chebulagic acid, and chebulinic acid. The interactions with M^pro of the top four docked compounds terflavin A, chebulagic acid, chebulinic acid, corilagin, including the native ligand X77 are visualized in Fig. 3.

Close

Fig. 1

Chemical structures of best-docked compounds and native ligand X77.

Export to PPT

Close

Fig. 2

Redocking of X77 in the original binding site of Mpro. Superimposed conformation of the experimental ligand with the native ligand, X77. Gray indicates X77 docked and pink indicates X77cocrystal.

Export to PPT

Close

Fig. 3

Interaction between the ligands and main protease (M^pro) (PDB: 6W63). (a) Terflavin, (b) Chebulagic acid, (c) Chebulinic acid, (d) Corilagin, (e) X77 or N-(4-tert-butylphenyl)-N-[(1R)-2-(cyclohexylamino)-2-oxo-1-(pyridin-3-yl)ethyl]-1H-imidazole-4-carboxamide.

Export to PPT

Terflavin A, chebulagic acid, chebulinic acid, corilagin are types of hydrolysable tannins, which are typically found on Terminalia species leaves (L.). Terflavin A and corilagin are ellagitannin, while chebulagic acid and chebulinic acid are benzopyran tannin. Tannins can inhibit SARS-CoV-2 through binding with catalytic dyad residues of M^pro (Chang et al., 2019; Das et al., 2020). Tannic acid even inhibits both M^pro and TMPRSS2 protease, effectively suppressing SARS-CoV-2 (Wang et al., 2020). The catalytic residues of 3CL^pro (or M^pro in the case of SARS-CoV-2) are Cys145 and His41, which are surrounded by other residues that confer substrate specificity (Rudrapal et al., 2021b; Singh et al., 2020).

Terflavin A interacts with M^pro on its active site by H-bonds at residues such as Asn142 (2.10 Å), Gly143 (2.13 Å), Ser144 (2.50 Å), Cys145 (2.25 Å), Glu166 (1.84 Å), Arg188 (2.13 Å), Gln189 (2.13 Å), and Thr190 (2.50 Å) and then H-carbon bonds at Met165 (1.66 Å) and Pro168 (1.74 Å) and pi-sulphur interactions at Cys44 (2.14 Å) and Met49 (2.43 Å) residues. Terflavin A exhibits H-bond with Cys145 (1.96 Å), but has no interaction with His41. Chebulagic acid forms H-bonds with residues Asn142 (1.67), Gly143 (2.18 Å), Glu166 (2.01 Å), Pro168 (2.97 Å), Arg188 (2.16 Å), and Gln189 (1.67 Å), H-carbon bonds in Met49 (1.90 Å) and Asp187 (1.61 Å) and pi-alkyl interactions in Met165 (1.72 Å). While the binding of free energy value is high, chebulagic acid does not interact with Cys145 or His41. Chebulinic acid forms H-bonds on the active site, namely residues Thr26 (1.89 Å), Phe140 (2.07 Å), Asn142 (2.18 Å), Cys145 (2.16 Å), Glu166 (2.26 Å), Gln189 (2.27 Å) and Thr190 (1.52 Å) as well as pi-alkyl, pi-pi stacked and pi-sigma interactions on Pro168 (2.04 Å), His41 (2.19 Å), and Gln189 (2.67 Å), respectively. Chebulinic acid, like Terflavin A, has only one H-bond with Cys145 (2.08 Å) and no interaction with His41. Corilagin forms H-bonds at residues Thr25 (1.87 Å), His41 (2.11 Å), Asn142 (1.03 Å), Cys145 (1.66 Å), and Glu166 (1.78 Å) then H-carbon bonds at Met165 (2.06 Å) and pi-alkyl and pi-sulphur interactions at residues Met49 (2.13 Å) and Cys145 (2.12 Å), respectively. As evidenced, H-bond interactions exist between corilagin and Cys145 (1.72 Å) and His41 (1.08 Å). As a native ligand, the X77 forms H-bonds at residues Asn142 (2.16 Å), Gly143 (1.89 Å), His163 (1.67 Å), and Glu166 (2.05 Å), then H-carbon bonds at Thr26 (2.10 Å) and Met165 (2.90) and pi-pi stacked and pi-sulphur interactions at residues His41 (2.17 Å), Met49 (2.04 Å), and Cys145 (1.14 Å), respectively. These findings support previous findings that the selected compounds bind to the catalytic residues of M^pro (Silva et al., 2021). The substrate-binding cleft contains a catalytic dyad composed of His41 and Cys145, which is located between domain I and II at the N-terminus, and a hydrogen-bonded water molecule bound to His41 acts as a third component of the catalytic triad (Tan et al., 2005).

3.2 Molecular dynamics simulation

MD simulations and other related techniques are on the verge to act as routine techniques for the discovery of drugs (De Vivo et al., 2016). In this study, the stability and dynamic behavior of the four ligand–protein (bound) complexes were verified by MD simulations for a 100 ns time. The RMSD, root mean square fluctuation (RMSF), radius of gyration (Rg) and solvent accessible surface area (SASA) graphs were plotted of the apoprotein, co-crystal ligand complex and four ligand–protein complexes (Junejo et al., 2021). The RMSD analysis provides information insight into the structural deviation and conformational stability and convergence during the time of the simulation. Thereby, in an equilibrated, stable system, the RMSD towards the end of MD simulations should be within 1–3 Å for a globular protein. The time-dependent RMSD graph was plotted for the backbone of the apoprotein, co-crystal ligand–protein complex and ligand (terflavin A, corilagin, chebulagic acid, chebulinic acid) protein (M^pro) complex. Fig. 4 (A) shows the RMSD values of apoprotein and co-crystal ligand–protein complex fluctuate in the beginning from 0.2 nm and at the end at around ∼0.3 nm, confirming the stability of the protein during 100 ns of simulation. At first, the RMSD values of the apoprotein and the protein–ligand complex are different, but after 50 nm, the RMSD values become stable and remain stable until 100 ns. Fig. 4 (B) shows the RMSD values of protein (M^pro) and ligand (terflavin A, corilagin, chebulagic acid, chebulinic acid) complex, indicated with colour code including M^pro-TerA (green), M^pro-Cor (pink), M^pro-Cgic (cyan), and M^pro-Cnic (indigo). The RMSD of M^pro-TerA fluctuates at the beginning until 70 ns (0.2 nm to ∼ 0.6 nm) and at 100 ns, it fluctuates at around ∼0.3 nm. The RMSD values of both the Protein-Corilagin (M^pro-Cor) and Protein-Chebulagic acid (M^pro-Cgic) complex are stable from the beginning (0.2 nm), and stay stable until 100 nm (0.25 nm), confirming the stability of both protein–ligand complexes during 100 ns of simulation. Whereas, Protein-Chebulinic acid complex fluctuats (between 0.2 nm and 0.3 nm) up to 70 ns, but after 70 nm, it becomes stable and stays stable (0.2 nm) until 100 nm. The overall RMSD results reveal that the binding of the ligands (corilagin, chebulagic acid, chebulinic acid) in the active site of M^pro is stable.

Close

Fig. 4

(a) Root-mean-square deviation (RMSD) analysis of apoprotein, and co-crystal ligand–protein complex; (b) Root-mean-square deviation (RMSD) analysis of ligand–protein complexes (M^pro-TerA, M^pro-Cor, M^pro-Cgic, and M^pro-Cnic); (c) Root-mean-square fluctuation (RMSF) analysis of apoprotein, and co-crystal ligand–protein complex; (d) Root-mean-square fluctuation (RMSF) analysis of ligand–protein complexes (M^pro-TerA, M^pro-Cor, M^pro-Cgic, and M^pro-Cnic); (e) SASA analysis and (f) Radius of gyration (Rg) of the apoprotein, co-crystal ligand–protein complex and protein–ligand complexes (M^pro-TerA, M^pro-Cor, M^pro -Cgic, and M^pro-Cnic).

Export to PPT

The RMSF study gives data regarding the amino acid residues fluctuation from the reference structure. Fig. 4 (C) shows that the co-crystal ligand fluctuation is almost similar to apoprotein, whereas Fig. 4 (D) shows the fluctuation of the ligand's binding to active site residues Thr24, Thr25, His41, Cys44, Met49, Tyr54, Phe140, Asn142, Gly143, Cys145, His163, His164, Met165, Glu166, Leu167, Pro168, Asp187, Arg188, Gln189, and Thr190 in each of M^pro TerA, M^pro-Cor, M^pro-Cgic, and M^pro-Cnic indicates that there is nearly similar RMSF fluctuation compared to apoprotein and co-crystal ligand–protein residues. Thus, analysis of the RMSF suggested a stable protein–ligand complex since there are no significant changes at the active site's reference location after successful binding by the apoprotein.

Further, the ‘sasa' function was used to calculate the SASA to examine the accessibility of solvent molecules on the protein surface. The SASA value of apoprotein, co-crystal ligand–protein complex, and protein–ligand complexes M^pro-TerA, M^pro-Cor, M^pro-Cgic, and Mpro-Cnic were stable during 100 ns MD. Thus, it indicates that ligand–protein complexes have comparable stability to apoprotein and co-crystal ligand–protein complexes [Fig. 4 (E)].

During MD simulations, a study on the radius of gyration (Rg) indicates the protein compactness. Rg of the C α atoms of the Mpro SARS-CoV-2 protein was computed to ascertain conformational changes in it. The Rg graph was plotted for protein using 100 ns trajectories as represented in Fig. 4 (F). The Rg plot indicates that the Rg value apoprotein (mean: 2.20 nm) and protein–ligand complex (mean: 1.23 nm) are nearly identical. Finally, Rg results suggest that all secondary structure of protein complexes is closely packed during 100 ns MD simulation.

Previous literature has highlighted the potential impact of medicinal plants and their phytoconstituents against SARS-CoV-2 (Rudrapal et al., 2021b). Studies have reported that thioflavonoids and epicatechin-3-O-gallate interacted with (binding sites) Mpro and PLpro, and were considered as SARS-CoV-2 inhibitors (Gogoi et al., 2021; Ismail et al., 2021). Thus, bioactive molecules bind to or interact with the active residues which are present at the binding cavity of Mpro and PLpro.

3.3 MM-PBSA binding free energy (BFE)

Python script MmPbSaStat.py provided in the g_mmpbsa was utilized to calculate average BFE and its standard deviation error. The energy released during the interaction of protein–ligand is represented in the form of BFE. It is a cumulative sum of Van der Waals, polar salvation, electrostatic, and SASA energy. Precisely, the MM-PBSA method estimates the binding free energy of a ligand–protein complex as the difference between the free energy of the complex and the free energies of the unbound components (receptor and ligand). More negative the BFE, the more effectively the ligand binds to the protein (Table 2). The results indicate that all the forms of energy contributed fairly to the interactions between M^pro and the selected molecules. M^pro-TerA complex showed the least binding BFE (−233.782 kJ/mol) among the selected bioactive molecules followed by M^pro-Cgic (−221.916 kJ/mol). The binding free energy (BFE) of the X77 (−160.625 kJ/mol) is less negative compared to the selected bioactive molecules. A comparison of the BFE of all the complexes reveals that, except for corilagin, all of the other selected bioactive molecules outperform X77 in inhibiting the protein M^pro of SARS-CoV-2.

3.4 Drug-likeness, ADMET and toxicity properties

Lipinski's rule of five (Ro5) is a set of rules used to assess the drug-likeness of a chemical compound that has a pharmacological activity that could lead to it being an orally active drug for human use (Lipinski et al., 2001; Kalita et al., 2020). The Ro5 specifies that the molecular weight must be less than 500 Da, the number of H-bond donors must be less than 5, the number of H-bond acceptor must be less than 10, the xlogP must be less than 5, and the molar refractivity must be between 40 and 130 (Fatimawali et al., 2021; Kashyap et al., 2016; Sefren Geiner et al., 2020). Regarding this statement, Table 3 indicates that only compound X77 appears to meet the Ro5 criteria. However, strict adherence to Ro5 is not required for naturally occurring compounds or phytocomponents. Lipinski specifically states that the natural products and substrates of biological transporters are not covered by Ro5 (Lipinski et al., 2001; Rudrapal et al., 2017; Rudrapal et al., 2018). It is, therefore important to note here that the desired pharmacological efficacy of Ayurvedic/herbal formulations is possibility due to the multi-component nature (different components) of such formulations.

On its way to the target site, the drug candidate undergoes pharmacokinetic events such as absorption, distribution, metabolism, and excretion (ADME). ADME parameters, in conjunction with toxicity, are critical for the discovery and development of potential drug candidates (Guan et al., 2019). Not only should a high-quality drug candidate be efficacious against the therapeutic target, but it should also exhibit appropriate ADMET properties at the dosage for therapeutic purposes. Based on the results of the admetSAR and ProTox analyses (Table 4), all of the compounds were from any signs of toxicity such as hepatotoxic, carcinogenic, or mutagenic. The acute oral toxicity analysis reveals that all compounds are classified as class III (500 mg/kg LD₅₀ 5000 mg/kg), which means that they are not toxic, except for chebulagic acid, which is classified as category II (50 mg/kg LD₅₀ 500 mg/kg), indicating a moderate level of toxicity (Herkenne et al., 2008; Li et al., 2014). Typically, animal experiments are used to determine a compound's toxicity, carcinogenicity, and mutagenicity. However, this takes a long time and necessitates the sacrifice of animals. Thus, prediction of toxicity within in silico technique is a faster and less expensive alternative, as it is based on existing data on a compound's toxicity. As a result, it is used to develop models capable of predicting the toxicity, carcinogenicity, and mutagenicity of novel compounds.

Bioavailability is a pharmacological term that refers to the level and amount to which a bioactive compound can be absorbed and made available (Kim et al., 2014). A compound's low oral bioavailability is a frequent issue during the drug design process. As a result, a drug may fail in clinical trials despite its high efficacy in vitro and in vivo experiments (Fukunishi et al., 2014). A compound is predicted to have low oral bioavailability if its bioavailability score is lower than 0.5. If the score is greater than 0.5, the compound is expected to have a high oral bioavailability (Rath et al., 2021). According to the findings, all compounds, except for X77, have a high oral bioavailability.

Toxic doses are frequently expressed as LD₅₀ (mg/kg body weight). LD50 (median lethal dose) is a particular dose wherein 50% of the participants die after being exposed to a compound (Erhirhie et al., 2018). A compound is considered highly toxic if the LD₅₀ value is higher than 300 mg/kg. If the LD50 is between 300 and 1,000 mg/kg, a compound is regarded as moderately toxic. If, on the other hand, the LD₅₀ value of a compound is between 1,000 and 5,000 mg/kg, it is considered to be less or slightly toxic. Thus, chebulagic acid, chebulinic acid, and X77 are classified as moderately toxic, whereas corilagin and terflavin A are classified as less toxic.