Phytoconstituents from ten natural herbs as potent inhibitors of main protease enzyme of SARS-COV-2: In silico study.

Background: Lack of treatment of novel Coronavirus disease led to the search of specific antivirals that are capable to inhibit the replication of the virus. The plant kingdom has demonstrated to be an important source of new molecules with antiviral potential. Purpose: The present study aims to utilize various computational tools to identify the most eligible drug candidate that have capabilities to halt the replication of SARS-COV-2 virus by inhibiting Main protease (Mpro) enzyme.
Methods: We have selected plants whose extracts have inhibitory potential against previously discovered coronaviruses. Their phytoconstituents were surveyed and a library of 100 molecules was prepared. Then, computational tools such as molecular docking, ADMET and molecular dynamic simulations were utilized to screen the compounds and evaluate them against Mpro enzyme.
Results: All the phytoconstituents showed good binding affinities towards Mpro enzyme. Among them laurolitsine possesses the highest binding affinity i.e. -294.1533 kcal/mol. On ADMET analysis of best three ligands were simulated for 1.2 ns, then the stable ligand among them was further simulated for 20 ns. Results revealed that no conformational changes were observed in the laurolitsine w.r.t. protein residues and low RMSD value suggested that the Laurolitsine-protein complex was stable for 20 ns.
Conclusion: Laurolitsine, an active constituent of roots of Lindera aggregata, was found to be having good ADMET profile and have capabilities to halt the activity of the enzyme. Therefore, this makes laurolitsine a good drug candidate for the treatment of COVID-19.

Introduction

Novel coronavirus disease (COVID) broke out in late 2019 in Wuhan city (China), from where, it spread to other parts of the world (Wang et al., 2020). On Mar 11th, 2020, the World Health Organization declared the COVID-19 as a pandemic and till Jan 7th, 2021, the virus had already infected approx. 85.1 million people and caused 1.8 million deaths across the globe. Countries such as USA, Brazil and India have been severely affected by this pandemic (Cucinotta and Vanelli, 2020; WHO, 2020) although at present other countries are suffering new outbreaks. Common symptoms of this disease include fever, cough, shortness of breath, loss of taste and pneumonia; whereas severe conditions include respiratory, gastrointestinal, hepatic and neurological complications among others (Pan et al., 2020; Singhal, 2020). After exposure, it generally takes 2–14 days to develop mild or severe symptoms. Direct contact with the infected person or respiratory droplets are considered as major modes of transmission (Huang et al., 2020).

COVID was caused by Novel Coronavirus 2019 (nCOV-19) or severe acute respiratory syndrome coronavirus-2 (SAR-COV-2) which belongs to beta coronavirus (COV) category of Coronaviridae family (Chen et al., 2020; Kanhed et al., 2020). This virus was thought to be sourced from bats (major reservoirs of α-COV and β-COV) but their transmission pattern has not been elucidated yet (Cui et al., 2019). nCOV-19 has a bi-lipid layer envelope that protects and stores genetic material i.e. RNA and N protein, Moreover, envelope surface has other proteins such as spike glycoprotein, membrane protein, envelope protein and hemagglutinin/esterase that help the virus to invade the host cell and perform other replication related tasks (Fig. 1 ).

Fig. 1

Diagrammatic representation of SARS-CoV-2.

Among them spike glycoprotein guides nCOV-19 to invade host cell by binding to the host's ACE-2 receptor (Fig. 2 ) with 1000 times more binding affinity than ACE-2 inhibitors discovered till now (He et al., 2020). After binding, the virus utilizes the receptor and enters the cytosol of the host cell where it un-coats and releases genetic material comprising of 30,000 nucleotides. Later, the genetic material gets translated by ribosomes to synthesize polyproteins. Maturation of these synthesized polyproteins is carried out by main protease enzyme (Mpro) which performs autolytic cleavage at not less than 11 conserved sites. These mature proteins merge again to generate a new viral envelope and accessory proteins that results in the formation of a new copy of virus (Mathewson et al., 2008).

Fig. 2

Schematic representation of life cycle of nCOV-19.

Therefore, suggesting two methods to halt the replication of virus. Firstly, by preventing the invasion of virus into the host cell by blocking ACE-2 receptor, and secondly, by inhibiting the activity of Mpro to pause the replication of virus. Furthermore, a study has reported that normal levels of ACE-2 receptors are required to combat the inflammatory response in the lungs. Angiotensin-II (an activator of Angiotensin-II type-1 receptor which on activation leads to acute lung injury) is degraded by ACE-2 receptor (Fig. 3 ). In case of nCOV-19 infection, ACE-2 receptors are occupied by the virus for its invasion into the host cell which indirectly increases the Angiotensin-II levels (Shenoy et al., 2011). This increase leads to activation of Angiotensin-II type-1 receptor causing vasoconstriction, proliferation and inflammation of the tissue leading to acute cell injury Kai and Kai (2020). Pulmonary arteries and veins have ACE-2 receptors whose blockage can cause constriction of arteries and veins, which can worsen the condition of the patient or may lead to lung failure also (Guignabert et al., 2018; Jia, 2016; Shenoy et al., 2011). On the other hand, the inhibition of Mpro is crucial for preventing the maturation of polyproteins that ultimately stops the spread of virus to other healthy cells. Moreover, the absence of closely related homologues of Mpro in humans also makes it a safe and potential target for the development of antivirals Hall and Ji (2020).

Fig. 3

schematic representation of injury caused by blocking of ACE-2 receptor (A) normal pathway (B) Altered pathway due to nCOV-19.

Currently, several Mpro inhibitors are proposed by researchers across the globe as potential drug targets for the treatment of coronavirus disease but they still lack optimum specificity and efficacy towards this protein. The conventional development of a new drug or vaccine requires a huge investment of money and time. From the last few years, computer-aided drug design (CADD) has played a prominent role in the development of numerous clinical drug candidates including Dorzolamide, Saquinavir and Aliskiren (Talele et al., 2010). Moreover, CADD provides more efficient and much cheaper and alternative path to design and develop bioactive molecules with high specificity and efficacy Surabhi and Singh (2018). In addition, various tools of CADD such as virtual screening, molecular docking, molecular dynamic simulations etc. are very effective in predicting the potential leads or drug candidates for various diseases. On the other hand, medicinal plants have successfully delivered a wide range of drug candidates to clinical practice. Thus, screening of phytoconstituents from antiviral medicinal plants using CADD is an efficient approach to identify potential leads against nCOV-19. Therefore, during literature survey, we came across ten natural herbs whose phytoconstituents or extracts were found to be beneficial in the treatment of previously discovered coronaviruses (Lin et al., 2014). The triterpene glycosides: saikosaponins (A, B, C and D), derived from the plants Bupleurum spp. (Apiaceae), Heteromorpha spp. (Apiaceae) and Scrophularia scorodonia L. (Scrophulariaceae) were found to prevent early-stage infection of HCoV-22E9 by inhibiting viral attachment and penetration (Cheng et al., 2006). Moreover, anti-SARS-COV effect was also shown by the extracts derived from Lycoris radiata (L'Hér.) Herb. (Amarylidaceae), Artemisia annua L., (Asteraceae), Pyrrosia lingua (Thunb.) Farw. (Polypodiaceae) and Lindera aggregata (Sims) Kosterm. (Lauraceae) (Li et al., 2005). Phenolic compounds from Isatis tinctoria L.(syn: Isatis indigotica Fortune ex Lindl.(Brassicaceae) and Torreya nucifera (L.) Siebold & Zucc. (Taxaceae) were found to inhibit nsP13 helicase and 3CL protease of SARS- COV (Lin et al., 2005; Ryu et al., 2010; Yu et al., 2012). Water extract from Houttuynia cordata Thunb. (Saururaceae) was also found to inhibit viral 3CL protease and RNA dependent RNA polymerase (Lau et al., 2008). As the extracts of these plants possess antiviral properties. Therefore, their phytoconstituents were chosen and screened against the Mpro using various in silico techniques in order to precisely identify the drug candidate for the treatment of COVID-19.

Materials and methods

A schematic workflow was followed to perform various in silico studies as shown in Fig. 4 .

Fig. 4

schematic representation of workflow.

Computational details

Molecular docking studies were carried out using DS Biovia's Discovery studio (DSBDS) software with hardware configuration of Intel(R) Core (TM) i7-3612QM CPU @ 2.10 GHz with 12 GB of RAM running on windows 10 × 64 bit operating system. Further, Molecular Dynamics Simulations were performed using DE Shaw's Desmond software with the hardware configuration of Intel(R) Pentium Dual-core E5400 CPU @ 2.70 GHz with Nvidia Geforce GT 710 and 4 GB of RAM running on Ubuntu 18.04 × 64 bit operating system.

Library preparation

As discussed above, ten natural herbs whose phytoconstituents or extracts were found to be beneficial in the treatment of previously discovered coronaviruses were chosen. Structures of their phytoconstituents were surveyed and downloaded from the PubChem database in structure data format (SDF) format (complete details of Phytoconstituents mentioned in supplementary excel sheet no. S1). These SDF files were then prepared using “prepare ligands” module and filtered using “Filter by Lipinski and Veber's Rules” module of DSBDS software. This process removed duplicate entries, computed isomers and tautomer, and generated 3D conformations and minimized them Dassault Systèmes BIOVIA (2020).

Protein preparation

3D structure of Mpro (PDB code: 5RE4) in complex with N-(4-methylpyridin-3-yl)acetamide inhibitor (NI) was downloaded from RCSB Protein Data Bank. The PDB file of Mpro was prepared using “prepare protein” module of DSBDS software. This process cleaned the PDB file, optimized the side-chains conformation for residues, removed water molecules, remodelled missing loop regions utilizing SEQRES information and protonated the structure at pH 7.4. Later the prepared PDB file was typed with CHARMM forcefield in order to score the docked poses and to calculate binding energies. In addition, the co-crystallized compound “NI” in 5RE4 PDB file which was used to detect the active binding site by generating binding sphere Dassault Systèmes BIOVIA (2020).

High throughput screening

Prior to docking the prepared ligands, the co-crystallized ligand NI was removed from the active binding and then High throughput screening of prepared ligands w.r.t Mpro was done using “LibDock” module of DSBDS. A total of 100 hotspots were setup with docking tolerance of 0.25 and a high-quality docking parameter. A Fast conformation method was chosen with a cut off value of 255 and minimization algorithm were set to “do not minimize” in order to rigidly dock the ligands on the active binding site of Mpro. Parallel processing was utilized with a batch size of 25 in order to speed up the docking process Dassault Systèmes BIOVIA (2020)

Scoring and calculation of binding energies of ligand-Mpro complex

After successful high throughput screening of ligands, an empirical fitting approach and knowledge-based statistical approach of “Score ligand poses” module of DSBDS was utilized. This process scored the confirmations and helped to speed up the process of estimation of binding energies of ligands Dassault Systèmes BIOVIA (2020). Later, these scored ligands were used by “Calculate binding energies” module of DSBDS to estimate the binding energy between the enzyme and the ligand

Prior to calculation of binding energy and conformation entropy “In situ ligand minimization” was turned on. This “In situ ligand minimization” utilizes smart minimizer to minimize each ligand up to 1000 steps with RMS gradient of 0.001 and removes any van der waals clashes of ligands. After minimization, binding energies were calculated which were used to prioritize the ligand poses Dassault Systèmes BIOVIA (2020).

ADMET analysis

All the ligands were screened using standard parameters of “ADMET descriptors” and “toxicity prediction (extensible)” modules of DSBDS for the calculation and estimation of various ADMET parameters (Dassault Systèmes BIOVIA, 2020; Han et al., 2019).

Molecular dynamics simulation

In order to mimic the biological system, a system of protein and ligand complex was generated by using “system builder” module of DE Shaw Desmond software. A SPC solvent model was selected with orthorhombic box shape with 10 × 10 × 10 dimensions. Later, Counter ions (Na+ or Cl−) for neutralization were added, system's concentration was set to 0.15 M and OPLS3e force field was applied to the system. Later, A short 1.2 ns MD simulation was performed between Mpro (5RE4) and best 3 ligands by using “Molecular Dynamics” module of DE Shaw Desmond. A 5-stage minimization of model system was performed prior to final MD simulation. These stages are: (a) stage 1 - simulates Brownian Dynamics NVT, T = 10 K, small timesteps, and restraints on solute heavy atoms, 100ps; (b) stage 2 - simulates NVT, T = 10 K, small timesteps, and restraints on solute heavy atoms, 12ps; (c) stage 3 - simulate NPT, T = 10 K, and restraints on solute heavy atoms, 12ps; (d) stage 4 - simulate, NPT and restraints on solute heavy atoms, 12ps; (d) stage 5 - simulate, NPT and no restraints, 24ps; (e) stage 6 – final simulation for 1.2 ns that covers about 200K steps. Result of these simulations were analyzed and the compound which showed minimum variations was then again evaluated by performing simulation for further 20 ns with same standard parameters (Bowers et al., 2006).

Results and discussion

The crystalline structure of Mpro (resolution 1.88 Å) was prepared and minimized using” prepare protein” module and its analysis showed that there were no missing residues within the protein. The Ramachandran plot showed no violation of amino acid's torsion angles and the energies of all the residues were within the low energy zones with glycine as an exception (Figure no. 5A). Hence, concluding that the protein has been minimized. Moreover, this 33.8 kDa protein consists of 304 amino acids and has 3 domains ranging from (a) 8–101 (red), (b) 102–184 (green) and (c) 201-303 (blue). The main catalytic site (grey colored) with 63 amino acids is present in a cleft between the first and second domain (Fig. 5 B). His-41, Cys-145 and His-164 residues are the crucial amino acids involved in the catalytic activity and these must be blocked to halt the activity of the Mpro (Jin et al., 2020). Therefore, these amino acids involved in catalytic activity of enzyme and native ligand NI were selected to generate a binding sphere at x: 8.660620, y: -1.671627, z: 24.136016 coordinates with a radius of 7.30 Å, this was done to give coordinates of active binding site to the software to perform docking.

Fig. 5

(A) Ramachandran plot for prepared Mpro (PDB code:5RE4), (B) Cartoon representation of Mpro with domain-1 (Red), domain-2 (Green), domain-3 (Blue) and active binding site represented by grey colored surface, (C) Showing (i) interacting residues with NI in original crystalized structure, (ii) interacting residues with NI after redocking.

Before performing high throughput screening of selected antiviral ligands, validation of docking protocol was done by re-docking co-crystalized ligand NI with Mpro. The docking results revealed that the interaction pattern of re-docked NI inhibitor was same as that in the crystallized state i.e., interacting with Glu-166, Cys-145 and Asn-142 amino acid residues with a RMSD value of 0.31659 and binding energy of -123.024 Kcal/mol (Fig. 5C).

After validation of the parameters, the prepared library which contains 30% of alkaloid, 21% terpenes, 14% of glycosides, 8 % each of terpenoids and saponins, 7% of flavonoids, 4% of phenolic and remaining 8% comprises of alcohols, phytosterols, flavone, alkane and other organic compounds was screened using “Lipinski's and Veber's rule of five”. This ruled out 30 molecules from the library as they were violating the parameters and then high throughput screening of remaining ligands were performed with the same validated parameters. This rigid docking of the ligands generated approx. 15,327 conformations which were then scored and sorted according to their dock score and binding energies values.

On analyzing the docking results, most of the molecules were found to be interacting with the amino acids involved in catalytic activity i.e. His-41, Cys-145 and His-164 (Fig. 6 A). The hydrophobic interaction of Cys-145 residue was found to be the most favorable, but also showed fewer hydrogen bonding interactions with the molecules. Whereas, other residues participating in hydrophobic interaction are Leu-27, His-41, Met-49, His-163, Met-165 and His-172. Other residues such as Thr-25, Thr-26, His-41, Cys-44, Ser-46, Phe-140, Leu-141, Asn-142, Gly-143, Ser-144, Cys-145, His-163, His-164, Met-165, Glu-166 and Gln-189 that showed hydrogen bonding interaction.

Fig. 6

(A) Histograms representing the count of different ligands interacting with various amino acid residues of Mpro‘s active binding site; (B) 2D interaction diagrams of top 10 molecules: (1) Laurolitsine, (2) Laurotetanine, (3) Ginkgetin, (4) Norisoboldine (5) Isoboldine, (6) Secoboldine, (7) Boldine, (8) Pseudolycorine, (9) Syringic acid, (10) Lauroscholtzine.

Later, the docked molecules were sorted according to their binding affinities towards Mpro enzyme and their binding energies of molecules were compared with co-crystalized ligand NI (values of remaining 70 molecules mentioned in supplementary excel sheet no. S3). Binding energies and libdockscore values of top 10 molecules are mentioned in Table 1 (Jin et al., 2020). The values signify that the molecules have capabilities to block the active site with good binding affinities. Among the top 10, laurolitsine showed the highest binding affinity i.e. -294.1533 kcal/mol which was approx. 2 times more than NI (-123.024 Kcal/mol) and considerably higher than N3 inhibitor (-262.071 Kcal/mol).

Table 1

Docking results and Lipinski's rule of five parameter values of top 10 molecules with their biological source.

S.no.	Ligand Name	Compound type	Plant source	Family	Docking results		Lipinski's rule of five parameters						References
					BE (kcal/mol)	Lib Dock Score	MW	Alog P	NRB	MR	HA	HD
1	Boldine	Aporphine alkaloid	Lindera aggregata (Sims) Kosterm	Lauraceae	-211.3621	74.59	288.318	-0.017	1	75.98	5	2	(Gan et al., 2009; Yang et al., 2020)
2	Ginkgetin	Biflavonoid	Houttuynia cordata Thunb.	Saururaceae	-269.3681	97.165	564.495	2.151	5	157.64	10	2	(Lau et al., 2008; Li et al., 2017a)
3	Isoboldine	Isoquinoline alkaloid	Lindera aggregata (Sims) Kosterm	Lauraceae	-214.4377	93.71	328.382	1.522	2	63.36	4	3	(Gan et al., 2009; Yang et al., 2020)
4	Laurolitsine	Aporphine alkaloid	Lindera aggregata (Sims) Kosterm	Lauraceae	-294.1533	65.433	313.348	0.613	2	78.36	4	2	(Gan et al., 2009; Yang et al., 2020)
5	Lauroscholtzine	Isoquinoline alkaloid	Lindera aggregata (Sims) Kosterm	Lauraceae	-202.2071	74.794	340.393	2.605	3	53.99	5	0	(Gan et al., 2009; Yang et al., 2020)
6	Laurotetanine	Isoquinoline alkaloid	Lindera aggregata (Sims) Kosterm	Lauraceae	-289.0029	76.903	327.374	0.839	3	67.36	4	1	(Gan et al., 2009; Yang et al., 2020)
7	Norisoboldine	Isoquinoline alkaloid	Lindera aggregata (Sims) Kosterm	Lauraceae	-242.8458	88.576	313.348	0.613	2	78.36	4	2	(Gan et al., 2009; Yang et al., 2020)
8	Pseudolycorine	Phenanthridine alkaloid	Lycoris radiata (L'Hér.) Herb.	Amaryllidaceae	-205.5765	60.259	196.157	-0.993	3	81.65	5	0	(Kihara et al., 1991; Park et al., 2021; Wang et al., 2009)
9	Secoboldine	Phenanthrene alkaloid	Lindera aggregata (Sims) Kosterm	Lauraceae	-212.6199	77.592	328.382	1.758	5	75.53	4	3	(Gan et al., 2009; Yang et al., 2020)
10	Syringic acid	Phenolic derivative	Isatis tinctoria L.(syn: Isatis indigotica Fortune ex Lindl)	Cruciferae	-205.0694	87.158	328.382	1.522	2	63.36	4	3	(Lin et al., 2005)

# BE: binding energy; MW: molecular weight; MR: molecular refractivity; HA: hydrogen acceptor, HD: Hydrogen bond donors; NRB: no. of rotatable bonds

Moreover, it was observed that among top 10 molecules, 8 molecules were alkaloid derivatives and one bioflavonoid and one phenolic compound which signifies that structure of alkaloids derivatives have the desired pharmacophoric features to perfectly bind and halt the activity of the enzyme. The 2D interaction diagrams of top 10 molecules that are laurolitsine (aporphine alkaloid), laurotetanine (isoquinoline alkaloid), ginkgetin (biflavonoid), norisoboldine (isoquinoline alkaloid), isoboldine (isoquinoline alkaloid), secoboldine (phenanthrene alkaloids), boldine (aporphine alkaloids), pseudolycorine (alkaloid), syringic acid (phenolic derivative), lauroscholtzine (isoquinoline alkaloid) are shown in Fig. 6B which suggest that the molecules were interacting with all major amino acids involved in the catalytic activity of the Mpro.

Further, the drug-likeliness of top ten molecules was evaluated through Lipinski's rule of five parameters and the values are summarised in Table 1 (values of all the molecules are mentioned in supplementary excel sheet no. S3). It was found that parameters of all molecules were within standard ranges except ginkgetin, which was found to be violating all the parameters. The molecular weight of molecules (except ginkgetin) was less than 500, log p was lower than 5.0 which suggest that molecules were not very lipophilic and their molar refractivity index was in the range of 98.4-111.85 which was within the limits of accepted value of 40–130 (Lipinski, 2004).

Further, the physicochemical and structural properties of all the docked molecules were utilized to check the mutagenicity and ADMET parameters. These parameters were evaluated by “ADMET descriptors” and “toxicity prediction (extensible)” modules of DSBDS (Han et al., 2019; Reddy et al., 2012). Ponnan et al., (2013) specified a set of rules (supplementary excel sheet no. S2) that were utilized to specify the threshold values of parameters. The values of predicted parameters of top 10 molecules were summarised in Table 2 , which showed that all compounds except ginkgetin had good intestinal absorption capacity and good aqueous solubility. None of the compound was found to inhibit cytochrome P450 2D6 enzyme. Whereas secoboldine and lauroscholtzine were found to have mutagenic nature, and all the compounds were found to be hepatotoxic. Moreover, the plot (Fig. 7 ) between Alog P98 and 2D polar surface area (PSA_2D) suggested that all the selected molecules possess good ADMET profile except Ginkgetin.

Table 2

Predicted ADMET profile of top 10 molecules.

Property/ molecule	Boldine	Ginkgetin	Isoboldine	Laurolitsine	Lauroscholtzine	Laurotetanine	Norisoboldine	Pseudolycorine	Secoboldine	Syringic acid
2D polar surface area (PSA_2D) with upper limit of 131.62 at 95 %	71.214	146.091	62.695	55.976	47.443	44.091	55.976	69.762	59.491	62.695
Absorption level (human intestinal absorption)	0	2	0	0	0	0	0	0	0	0
AlogP98	-0.017	2.151	1.522	0.613	2.605	0.839	0.613	-0.993	1.758	1.522
AMES toxicity	Non-Mutagen	Non-Mutagen	Non-Mutagen	Non-Mutagen	Mutagen	Non-Mutagen	Non-Mutagen	Non-Mutagen	Mutagen	Non-Mutagen
Aqueous solubility level	4	3	3	3	2	3	3	4	3	3
Blood brain barrier penetration level (BBB)	3	4	3	3	2	3	3	3	3	3
CYP2D6 inhibition	0	0	0	0	0	0	0	0	0	0
hepatotoxicity	1	1	1	1	1	1	1	1	1	1
Plasma protein binding level (PPB)	0	0	0	2	2	2	2	2	0	0

Fig. 7

plot between Alog P98 and PSA 2D.

After analyzing interaction pattern and ADMET profile of the molecules, top 3 molecules were chosen for Molecular Dynamics (MD) simulation to understand the time-dependent ligand-protein complex stability. Although ginkgetin violates all the ADMET parameters, it possesses good binding energies with Mpro which led us to choose this molecule among the other two i.e. laurolitsine and laurotetanine for the MD studies. Preliminary evaluation was done by a short simulation of 1.2 ns performed using DE Shaw Desmond which simulated the Newtonian dynamics of the generated model system. This simulation generated various trajectories of the particles, coordinates, velocities, and energies which were then analyzed statistically to get the desired results. Initial ligand-Mpro complex pose was considered as a reference for the calculation of Root Mean Square Deviation (RMSD) and other statistical parameters.

Results of MD were analyzed and on aligning ligand RMSD with respect to (w.r.t.) protein's backbone RMSD revealed that the variations of top 3 ligands were reported to be less than 3 Å (Fig. 8 i). Moreover, ligand RMSD of laurolitsine showed very low deviation as compared to RMSD's of laurotetanine and ginkgetin. For laurolitsine, RMSD varied near to the value of 0.20 for the first 200 ps, but later varied between a range of 0.40 to 0.60 till 1200 ps. A slight increase was observed but was within the 3 Å margin, hence suggesting that the protein-ligand complex was stable for 1.2 ns without showing any major orientational changes.

Fig. 8

Representing results for 1.2 ns MD simulation of top 3 molecules (i) RMSD plots of top 3 potent molecules aligned with the backbone of Mpro, (ii) RMSF plot protein w.r.t ligands, (iii) Type of interaction of ligand with Mpro, (iv) Timeline of ligand residues contact varying w.r.t. time, (v) Represents fluctuation in rGyr, intra HB, MoISA, SASA and PSA plots of ligands

Furthermore, the local changes within the protein residues were understood by calculating protein's root mean square fluctuations (RMSF) where green colored vertical lines in the graph shows contact of ligand's atoms with amino acid residues of protein (Fig. 8ii). The results revealed that protein's backbone showed no significant fluctuations w.r.t. these compounds. RMSF of backbone w.r.t laurolitsine was between 0.4 Å and 1.6 Å, but the fluctuation was minimal where the amino acid residues were in contact with laurolitsine atoms. Moreover, RMSF of laurotetanine and ginkgetin shared the same pattern of fluctuation. Thus, it can be concluded that there were no major fluctuations in the backbone of the protein w.r.t. ligands.

Later, interaction patterns for these ligands were studied w.r.t. time (Fig. 8iv) and it was observed that the interactions circulated within the residues Thr-26, Asn-142, Gly-143, Ser-144, Cys-145, Glu-166 and Gln-189 throughout the simulation. At a time, a maximum of 9 residues were observed to be in contact with the ligand. In the case of laurolitsine and laurotetanine, the interaction with Glu-166 was observed to be consistent throughout the timeline but was inconsistent in the case of ginkgetin. Moving towards the interaction fraction of ligands with protein residues (Fig. 8iii), most interactions in case of laurolitsine were through hydrogen bonding and water bridges; but showed few hydrophobic and ionic interactions. Whereas, in case of laurotetanine and ginkgetin, water bridges were more prominently observed than other interactions. This concludes that laurolitsine has more tendency to interact with the residues of the Mpro directly as compared to laurotetanine and ginkgetin.

Apart from above observations, solvent accessible surface area analysis, Radius of gyration (rGyr), Molecular Surface Area (MoISA) and Polar Surface Area (PSA) further supported the stability of protein-ligand complex (Fig. 8v). But, ginkgetin showed quite significant fluctuations in intramolecular hydrogen bonding which may be due to its flexibility.

The three compounds were found to be stable but laurolitsine showed least variations therefore the stability of Mpro-laurolitsine complex was further analyzed at 20 ns. The same previously used parameters were employed for MD simulation and initially docked pose was considered as a reference. On aligning the RMSDs of laurolitsine with Mpro protein, the RMSD variation pattern was observed to be the same as was in initial 1.2 ns MD simulations, after 10 ns the variation lowered to 0.4 (approx.) which suggested that the system has further moved toward equilibrium which further support the stability of protein-ligand complex (Fig. 9 A). On the other hand, RMSF of protein was more fluctuating as compared to 1.2 ns simulation which was smoother (Fig. 9B). But the fluctuations observed were within the range of 0.5 and 2.0 Å that is an acceptable range. Loss of a few contact points with the ligand was observed. In conclusion, these RMSD and RMSF patterns confirmed stability of the laurolitsine and Mpro complex.

Fig. 9

Representing results of 20 ns MD simulations of Laurolitsine (A) Aligned RMSD of Laurolitsine and Mpro Complex, (B) Fluctuation of protein residues during the simulation RMSF-P, (C) Representing protein ligand contact points via amino acids residues, (D) Histograms show the type of protein ligand interaction, (E) 2D interaction diagram of Laurolitsine with protein residues and (F) Represents plots for rGyr, intra HB, MoISA, SASA and PSA values of Laurolitsine.

The interaction timelines of laurolitsine and Mpro residues w.r.t. time revealed that the interactions of laurolitsine with Glu-166 (contacts between 3 or 4) and Asn-142 (contacts between 1 or 2) residues were stable and observed to be consistent throughout the simulation (Fig. 9C). Other residues interacting with laurolitsine were Thr-25, Thr-26, His-41, Ser-46, Met-49, Asn-119, Phe-140, Leu-141, Gly-143, Ser-144, Cys-145, His-163, His-164, Met-165, Pro-168, Arg-188, Gln-188, Gln-189, Thr-190 and Gln-192. Histogram plot of interaction fractions showed that laurolitsine was mainly interacting via water bridges and H-bonds but have fewer hydrophobic and ionic interactions with the Mpro (Fig. 9D and 9F). Other statistical parameters, such as SASA analysis, rGyr, MoISA and PSA were within the acceptable ranges that added their support for stability of Mpro-ligand complex .

Conclusion

All the phytoconstituents interacted well with Mpro protein and possesses good binding affinities but alkaloid-class derivatives were found to be dominating than the others. Among them, laurolitsine was the most promising compound with highest binding energy and its complex with Mpro was found to be stable for 20 ns. In addition to this, laurotetanine and ginkgetin were also among the top 3 drug candidates as they also have good binding affinity towards Mpro protein with good ADMET profile but their complex with Mpro showed much more variations than the laurolitsine- Mpro complex. Apart from this laurolitsine also possesses a good ADMET profile and it is an active constituent derived from roots of Lindera aggregata whose extract has already proven for its anti-SARS-COV effect. All these properties make laurolitsine a promising candidate for the treatment of COVID-19.

CRediT authorship contribution statement

Nitish Kumar: Conceptualization, Methodology, Resources, Validation, Writing – original draft, Writing – review & editing. Atamjit Singh: Data curation, Supervision, Formal analysis, Investigation, Writing – review & editing. Harmandeep Kaur Gulati: Data curation, Supervision, Formal analysis, Investigation, Writing – review & editing. Kavita Bhagat: Data curation, Supervision, Formal analysis, Investigation, Writing – review & editing. Komalpreet Kaur: Data curation, Supervision, Formal analysis, Investigation, Writing – review & editing. Jaspreet Kaur: Data curation, Supervision, Formal analysis, Investigation, Writing – review & editing. Shilpa Dudhal: Data curation, Supervision, Formal analysis, Investigation, Writing – review & editing. Amit Duggal: Data curation, Supervision, Formal analysis, Investigation, Writing – review & editing. Puja Gulati: Data curation, Supervision, Formal analysis, Investigation, Writing – review & editing. Harbinder Singh: Data curation, Supervision, Formal analysis, Investigation, Writing – review & editing. Jatinder Vir Singh: Data curation, Supervision, Formal analysis, Investigation, Writing – review & editing. Preet Mohinder Singh Bedi: Data curation, Supervision, Formal analysis, Investigation, Writing – review & editing.

Declaration of Competing Interest

None.