Identification of phytochemical inhibitors against main protease of COVID-19 using molecular modeling approaches.

Severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2) is a novel corona virus that causes corona virus disease 2019 (COVID-19). The COVID-19 rapidly spread across the nations with high mortality rate even as very little is known to contain the virus at present. In the current study, we report novel natural metabolites namely, ursolic acid, carvacrol and oleanolic acid as the potential inhibitors against main protease (Mpro) of COVID-19 by using integrated molecular modeling approaches. From a combination of molecular docking and molecular dynamic (MD) simulations, we found three ligands bound to protease during 50 ns of MD simulations. Furthermore, the molecular mechanic/generalized/Born/Poisson-Boltzmann surface area (MM/G/P/BSA) free energy calculations showed that these chemical molecules have stable and favourable energies causing strong binding with binding site of Mpro protein. All these three molecules, namely, ursolic acid, carvacrol and oleanolic acid, have passed the ADME (Absorption, Distribution, Metabolism, and Excretion) property as well as Lipinski's rule of five. The study provides a basic foundation and suggests that the three phytochemicals, viz. ursolic acid, carvacrol and oleanolic acid could serve as potential inhibitors in regulating the Mpro protein's function and controlling viral replication. Communicated by Ramaswamy H. Sarma.

Introduction

Coronavirus disease (COVID-19) is a respiratory infectious disease caused by a novel virus strain, SARS-CoV-2 (Boopathi et al., 2020; Hemida & Ba Abduallah, 2020; Salata et al., 2019; Sarma et al., 2020; Seah & Agarwal, 2020; Su et al., 2016). In the past two decades, two other coronavriruses have caused global outbreaks, namely SARS-CoV (2002–2003) and Middle East respiratory syndrome coronavirus (2012–present) (de Wit et al., 2016; Gupta et al., 2020; Wang et al., 2020; Wu et al., 2020; Yuan et al., 2020). The Coronavirus Study Group (CSG) taxonomists working under the aegis of International Committee on Taxonomy of Viruses (ICTV) coined the nomenclature of SARS-COV-2 based on its 82% identity to the SARS coronavirus (SARS-CoV) genome (Coronaviridae Study Group of the International Committee on Taxonomy of Viruses, 2020; Hasan et al., 2020). Whole genome functional analysis revealed that both viruses phylogenetically belong to clade b of the genus Betacoronavirus (Chan et al., 2020; Muralidharan et al., 2020; WHO, 2020). The very first case of the novel COVID-19 was originally reported in Wuhan, Hubei Province, China, and has quickly spread over the 212 countries of the world (Elmezayen et al., 2020; Enayatkhani et al., 2020; Mackenzie & Smith, 2020; Xu et al., 2020). This contagious disease has led to over 4, 098, 018 confirmed cases and 283,271fatalities as on May 12, 2020; https://covid19.who.int/). The number of cases across the globe is increasing abruptly, and so far, there is no standard drug has proved to be effective for COVID-19 disease that have high mortality rate in immunocompromised patients (Aanouz et al., 2020; Bhatraju et al., 2020; Weiss & Murdoch, 2020). On 30 January, this respiratory infectious disease has been declared as a Public Health Emergency of International Concern by the World Health Organization (WHO) (Mackenzie & Smith, 2020; https://www.who.int/).

The SARS-CoV-2 genome is about 30 kb in size (29, 903 nucleotides) and encodes as many as 14 open reading frames (ORFs) (Chen et al., 2020; Elfiky & Azzam, 2020; Gorden et al., 2020; Woo et al., 2005). At the 5′ end of the viral genome, Orf1a/Orf1ab encodes several proteins, which are auto-proteolytically processed into 16 non-structural proteins (Nsp1-16) and form the replicase/transcriptase complex (RTC). Whereas, 3′ endencode structural viral proteins; spike (S), membrane (M), envelope (E) and nucleocapsid (N), and nine putative accessory factors. On the other hand, protease enzyme called the Mpro or also 3CLpro, play a crucial role in the life cycle of COVID-19 replication and maturation (Joshi et al., 2020; Khan et al., 2020; Pant et al., 2020; Wu et al., 2020). It has been demonstrated that no part of COVID-19 is more exposed than its main protease, as belongs to non-structural class proteins of the viral genome (Enmozhi et al., 2020; Zhang et al., 2020). Zhang et al., 2020 determined the three-dimensional crystal structure, at 1.75 Å resolution, of the Mpro of SARS-CoV-2 using x-ray crystallography.

Like other coronaviruses, the 3D structure of 3CLpro of SARS-COV-2 possesses three functional domains (Jin et al., 2020; Khan et al., 2020). The length of domain I, II, and III ranges from 8-101, 102-184, 201-306 amino acid residues, respectively. Among them, domain I and II are essentially beta-barrels and similar to the chymotrypsin. While, structural composition of domain III mainly consists of alpha-helices (Anand et al., 2005; Jin et al., 2020; Lu et al., 2006; Yang et al., 2006), protease has been characterized as one of the potential drug targets among coronaviruses (Anand et al., 2003). In this sense, many recent studies suggested that the selection of various FDA-approved antiviral compounds may yield promising results against COVID-19 infection (Chang et al., 2016; Chang et al., 2020; Contini, 2020; Gonzalez-Paz et al., 2020; Khan et al., 2020; Elfiky, 2020a, 2020b; Islam et al., 2020; Sinha et al., 2020; Wahedi et al., 2020; Umesh et al., 2020; Das et al., 2020; Abdelli et al., 2020). Chang et al., 2020 reported that chloroquine, an older antimalarial drug, has ability to inhibit the viral 3CL-protease activity. As crystal structure of protease provides a foundation for design of improved α-ketoamide inhibitors (Zhang et al., 2020), the capability of chloroquine to inhibit the protease activity its uses has been recommended in different countries including China, USA and India for the treatment of COVID-19 (Devaux et al., 2020; Gonzalez-Paz et al., 2020; Wang et al., 2020). However, many studies questioned the safety and reported the severe adverse effects of chloroquine (Kaisari & Borruat, 2020; Wang et al., 2020).

Ancient Indian scriptures including Rig-Veda, Atherveda, and Charka Sanhita demonstrated abundant benefits of plants for the treatment of various human aliments (Kumar et al., 2018). Plants are a remarkable natural source of high value alkaloids, flavonoids, phenols, chalcones, coumarines, lignans, polyketides, alkanes, alkenes, alkynes, simple aromatics, peptides, terpenes, and steroids. In the current era of drug discovery, enormous medicinal properties of plants allows the researchers to exclusively use them for the discovery of drug-like natural molecules (Jee et al., 2018; Kumar et al., 2018; Panchangam et al., 2016).

Ursolic acid (3-β-3-hydroxy-urs-12-ene-28-oic-acid) and oleanolic acid (3β-hydroxyolean-12-en-28-oic acid) are pentacyclic triterpenoid compounds with a widespread occurrence throughout the plant kingdom (Pollier & Goossens, 2012; Woźniak et al., 2015). Both molecules enrich various therapeutic properties such as antibacterial, antiviral, anticancer, antioxidant and tantimycotic activity (Jesus et al., 2015). Previous in vitro studies reported that these molecules exhibit antiviral activity against rotavirus, HIV, the influenza virus, hepatitis B and C viruses (Jesus et al., 2015; Khwaza et al., 2018; Tohmé et al., 2019).

Carvacrol (2-Methyl-5-(propan-2-yl)phenol) is a monoterpenoid phenol, possesses a wide range of strong antimicrobial and antiviral activity (Gilling et al., 2014; Kamalabadi et al., 2018; Marinelli et al., 2019). It is a major constituent of essential oil of plants of Labiatate family including oregano and thyme, and has been emerged as active molecule for therapeutic purpose (Hyldgaard et al., 2012). Using in-vitro methods, antiviral effects of carvacrol has been tested and validated on herpes simplex virus type 1, retrovirus and human respiratory syncytial virus (Kamalabadi et al., 2018).

Besides the uses of various FDA-approved antiviral compounds as mentioned above, there are many in-silico studies have been performed to screen the novel phytochemical molecules as a potential inhibitors of main protease of SARS-CoV-2 or develop new drugs against COVID-19 (Adem et al., 2020; Chandel et al., 2020; Gentile et al., 2020; Gonzalez-Paz et al., 2020; Khaerunnisa et al., 2020; Khan et al., 2020; Qamar et al., 2020; Sharma & Kaur, 2020; Sun et al., 2020). But there is not a single report available for the function of phytochemicals namely, ursolic acid, carvacrol and oleanolic acid derived from the antiviral herbs for the treatment of COVID-19. With the advent of molecular modeling approaches, targeted drug design may be possible and functional proteins of SARS-CoV-2 could be targeted with natural compounds to develop an effective treatment for COVID-19.

In the present study, we have targeted the protease of SARS-CoV-2 virus using available molecular modelling based methods and studied the interactions with selected natural compounds (ursolic acid, carvacrol and oleanolic acid) by molecular docking and molecular dynamics simulations followed by molecular mechanic/generalized Born/Poisson-Boltzmann surface area (MM/G/P/BSA) validation.

Material and methods

A flow chart of pipeline used in present study is summarized in Figure 1.

Figure 1.

Flowchart of pipeline used in present study to identify the phytochemical based inhibitors of Mpro.

Preparation of protease

The crystal structure of COVID-19 virus main protease (Mpro) in complex with Z45617795 (PDB ID: 5R7Y) was solved by PanDDA analysis group (https://www.rcsb.org/structure/5R7Y). This solved protein structure (PDB ID: 5R7Y) was extracted from the RCSB-Protein Data Bank (Berman et al., 2000; Burley et al., 2019). Crystal structure of Mpro have dimer in form with fragment of N-(2-phenylethyl) methanesulfonamide (Z45617795) compound provides a model for identifying potent inhibitors to target COVID-19 virus Mpro through in-silico study. The protein structure was prepared by removing water atoms, hetero atoms and adding polar hydrogen atom and kollman charges on it.

Ligand selection

Chemical structure of carvacrol (CID_10364), oleanolic acid (CID_10494) and ursolic acid (CID_64945) were extracted from PubChem database (Kim et al., 2019).

Molecular docking

In order to find out the potential drug targets, ligand based molecular docking between phytochemical compounds i.e. (carvacrol, oleanolic acid and ursolic acid) and Mpro protein were performed using AutoDock v4.2 (Morris et al., 2009). For docking experiments, the amino acid residues including Thr24, Thr26, Asn119, Phe140, Gly143, Cys145, His163, His164, Glu166, Gln189, and Thr190 were used as the active sites. Default mode parameters were selected in AutoDock during docking analysis.

Drug likeness

Absorption, distribution, metabolism, and excertion (ADME) calculations are important aspects of drug designing. All three ligand molecules were analysed based on the Lipinski’s rule of five using (Lipinski et al., 2001), Veber’s rule (Veber et al., 2002), Egan’s rule (Veber et al., 2002) and polar surface area (TPSA), and number of rotatable bonds (Daina et al., 2017), were calculated using SWISSADME tool (Daina et al., 2017).

Molecular dynamics (MD) simulations

The structural and dynamics transition at atomistic level in the Mpro of COVID-19 upon binding of small molecules were investigated by using MD simulations. MD simulations were performed using GROMOS96 43a1 force field embedded in GROMACS 5.1.1 suite on LINUX based platform (Berendsen et al., 1995). For the MD simulation, we followed the protocol described by Gajula, 2008; Gajula et al., 2016; Kumar et al., 2018. PRODRG server (Schüttelkopf & van Aalten, 2004) was used to generate the topology files of small molecules. The protein complexes were solvated in the dodecahedron box with simple point charge (SPC) waters, and 4 Na+ were added for the overall electrostatic neutrality of the system. Energy minimization of the system was performed by using the steepest descent algorithm for 50,000 iteration steps and cut-off up to 1000 kJmol−1 to reduce the steric clashes. The minimization of the system was performed at two phases for 50,000 steps, where the first phase equilibrated in two different phases for 50,000 steps. The first phase of equilibration was performed with a constant number of particles, volume, and temperature (NVT), each step 2 fs, the fallowed second phase of equilibration was performed with a constant number of particles, pressure, and temperature (NPT), the ensemble at 300 K. LINCS algorithm was utilized for covalent bond constraints in the equilibration steps. For the calculation of Lennard-Jones and Coulomb interactions, a 1.4 nm radius cut-off was used. Long-range electrostatics was calculated by using the Particle Mesh Ewald (PME) method with a Fourier grid spacing of 1.6 Å. The temperature inside the box was regulated by using V-rescale, a modified Berendsen temperature coupling method. Parrinello-Rahman pressure coupling method was utilized in NPT equilibration. The final production step of molecular dynamics simulation was carried out for 50 ns, each step of 2 fs. Trajectories were saved, and results were analyzed using XMGRACE. Root mean square deviation (RMSD) variation in protein backbone was calculated by using g RMS tool, which utilizes the least-square fitting method. Overall root mean square fluctuation (RMSF) in the atomic positions of protein C backbone was calculated by using the g rmsf tool. A rough measure of the compactness factor of protein during the course of the simulation was estimated by using the g gyrate tool of GROMACS. Gmxsasa was used for computation of the total solvent accessible surface area (SASA). Hydrogen bonds were calculated with 3.5 Å distance cut-off by using g h bond, and the distribution of intermolecular hydrogen bond lengths throughout the simulation was also analyzed.

MM/G/P/BSA binding free energy calculations

Molecular Mechanic/Poisson-Boltzmann Surface Area (MMPBSA) method was used to obtain the binding free energy of the interaction between ligand-protein complexes (Aldeghi et al., 2017). This employs ensembles derived from MD simulation. The g mmpbsa application of GROMACS module was used for the calculation of different components of the binding free energy of the Mpro and ligand complexes. Here, the binding energy is an average of three energetic terms, i.e. potential energy in the vacuum, polar-solvation energy, and non-polar solvation energy, respectively. The snapshots at every 100ps between 40 and 50 ns were collected, and MMPBSA was performed to predict the binding energy.

In the MMPBSA calculation, the binding free energy between a receptor and a ligand was calculated using following equations:

Computational facility details

The MD simulations and corresponding energy calculations were carried out on HP Gen7 server with 48 Core AMD processors and 32GB of RAM.

Results and discussion

The molecular docking approach to identify the drug targets has become one of the most popular methods for ligand-based computer-aided drug discovery (LB-CADD). In current era, with this approach, big data of drug libraries can be analyzed and annotated quickly and immense amount of energy, time, and costs related to CADD can be saved (Wadood et al., 2013; Yu & MacKerell, 2017). Currently, there are no effective treatments available to cure COVID-19 virus, and hence, identification of potential drug targets is urgently needed.

We used molecular modeling approach with molecular docking and MD simulation to identify potential phytochemcials active against the Mpro protein of COVID-19. These screened natural compounds may pave way for the development of drugs against COVID-19. On the basis of AutoDock binding affinity, phytochemicals namely, carvacrol, oleanolic acid and ursolic acid have shown satisfactory interactions with active site residues. These compounds have been found to have binding energy of −4.0 kcal/mol, −6.0 kcal/mol, −5.9 kcal/mol, respectively. The docked ligand molecules with protease are shown in Figure 2. Table 1 represented the 2D and 3D structures of docked ligand molecules along with AutoDock score. Hydrogen bonds play an important role in determining its specificity and affinity within protein-ligand complexes (Sakkiah et al., 2012).

Figure 2.

Schematic representation of molecular docking between Mpro and phytochemicals; (a) 3D structure of coronavirus derived from the RCSB-Protein Data Bank (PDB); (b) crystal structure of main protease of COVID-19 obtained from the PDB with PDB ID:5R7Y (c)interaction between Mpro and oleanolic acid with −6.0 kcal/mol docking energy; (d) interaction between Mpro and carvacrol with docking energy −4.0 kcal/mol; (e) interaction between Mpro and ursolic acid with −5.9 kcal/mol docking energy. Interactions were visualized using maestro and discovery studio programs.

Table 1.

Molecular docking results of carvacrol, oleanolic acid and ursolic acid.

		Chemical Structure
Compound Name	PubChem ID	2D	3D	Docking Energy (AutoDock)
Carvacrol	CID_10364			−4.0 kcal/mol
Oleanolic acid	CID_10494			−6.0 kcal/mol
Ursolic acid	CID_64945			−5.9 kcal/mol

In order to evaluate the affinity of the ligand molecules with the Mpro protein, we also studied the residues of protease involved in forming hydrogen bonds with the phytochemicals and the strength of these bonds. Amino acids residues participated in hydrophobic interactions between protease and ligand molecules were also investigated. Molecular interactions (hydrogen bonds and hydrophobic interactions) play a key role in giving shape and stabilizing the docking complexes (Wade & Goodford, 1989).

With AutoDock binding energy −4.0 kcal/mol, Gly143 amino acid was found to be involved in forming hydrogen bond with carvacrol (Figure 2). In addition to this, six other residues (Leu27, His41, Met49, Asn142, Cys145, Met165) are involved in forming hydrophobic interactions. Oleanolic acid compound make complex using AutoDock with binding energy −6.0 kcal/mol. The Gln189 amino acid formed the hydrogen bond, and two amino acids (Cys145, His163) were involved in hydrophobic interactions. Autodock binding energy of ursolic acid compound was −5.9 kcal/mol. The Ser46 amino acid at the active site of protease formed hydrogen bond. Nine amino acids (Thr24, Thr25, Thr26, Cys44, Thr45, Asn142, Gly143, Cys145, Glu166) participated in the formation of van der Waals interactions (Figure 2). Carvacrol having less binding energy as compared to oleanolic acid and ursolic acid molecules, the binding mode of interaction was found to be reasonably good. Previous studies reported a similar trend of the presence of binding pockets in main Mpro of SARS-CoV, which confirms our study (Aanouz et al., 2020; Khan et al., 2020; Li et al., 2020; Muralidharan et al., 2020). Aanouz et al. (2020) reported the β–Eudesmol molecule as a potential inhibitor of Mpro of COVID-19 based on its significant antibacterial and antiviral power well documented in literatures. However, this compound showed the low affinity interaction with Mpro and formed the hydrogen only with Thr111.

In our earlier study we have reported the ursolic acid carvacrol molecules as potential inhibitors of Small heat shock protein16.3 (sHSP16.3) of Mycobacterium tuberculosis (MTB) (Jee et al., 2018). Docking study of ursolic acid with sHSP16.3 formed two strong hydrogen bonds with Glu56 and Pro58, while carvacrol was found to be bound with only one hydrogen bond Glu92. Binding affinity of both compounds with sHSP16.3 confirms the present study results. From a biological point of view, reported phytochemicals, carvacrol, oleanolic acid and ursolic acid which are proposed as potential inhibitors of the Mpro having a significant antiviral activity with evident to bibliographical research and performed in-vitro studies (Gilling et al., 2014; Jesus et al., 2015; Kamalabadi et al., 2018; Khwaza et al., 2018; Marinelli et al., 2019; Tohmé et al., 2019).

Evaluation of drug likeness

Carvacrol, oleanolic acid and ursolic acid have respectively the following molecular weights (150.22, 456.70, 456.70) g/mol, all three molecules have molecular weight ≤ 500 g,/mol which follows the Lipinski’s rule. The carvacrol molecule has been found to have topological polar surface (TPSA) as 20.23 Å2, while oleanolic acid and ursolic acid molecules were found to have 57.53 Å2. According to the rule lowest TPSA values always produce good results; therefore, we noted that selected molecules are better behaved than the co-crystallized ligand (Daina et al., 2017; Lipinski et al., 2001). The (LogP) values of carvacrol, oleanolic acid and ursolic acid were found with the range of 2.82, 6.06, 5.94, respectively. Predicted LogP values depict that these molecules can be absorbed in the body. All three molecules present number of hydrogen bond donors: ≤5, a number of hydrogen bond acceptors ≤10 and also molar refractivity values between 48.01, 136.65, 136.91, respectively. The hydrogen bonding and molar refractivity calculations showed that these molecules validate the five Lipinski’s rule. These molecules also follow the Veber’s rule which denotes the oral bioavailability of drug molecules. Since in some CAD molecules cannot be synthesized, synthetic accessibility (SA) is a key aspect of drug designing. The carvacrol, oleanolic acid and ursolic acid molecules were found to have SA score 1, 6.08, 6.21, respectively. According to SA method molecules having score 1 is easy to synthesize it, on the other hand, score 10 represent difficulties to synthesize the molecules. All three molecules having less than 10 score, so they can be easily synthesized. The results of the ADME calculations are listed in the Table 2.

Table 2.

Results of the phytochemical molecules druglikeness properties.

Drug Likeness Properties	Carvacrol	Oleanolic acid	Ursolic acid
Molecular weight g/mol	150.22	456.70	456.70
Concensus Log Po/w	2.82	6.06	5.94
Num. H-bond acceptors	1	3	3
Num. H-bond donors	1	2	2
Molar Refractivity	48.01	136.65	136.91
Lipinski	Yes	Yes	Yes
Veber	Yes	Yes	Yes
Bioavailability score	0.55	0.56	0.56
Synthetic accessibility (SA)	1.00	6.08	6.21
TPSA (Ų)	20.23	57.53	57.53
No of rotatable bonds	1	1	1
Solubility	1.46e-01	3.45e-04	9.72e-04

MD Simulations

MD simulations are one of proven in-silico methods for obtaining dynamic data at atomic spatial resolution and picoseconds or finer temporal resolution (Benson & Daggett, 2012; Gajula et al., 2016). The main protease docked complexes with phytochemical compounds; carvacrol, ursolic acid and oleanolic acid simulation study done for 50 ns simulation period to analyze the stability of these docked phytochemical compounds in binding region of Mpro.

The Mpro with compound carvacrol (black) showed stable and constant range of RMSD between 0.23 nm to 0.30 nm and similar results showed by ursolic acid (blue), only slight changes showed in starting of the simulation. The oleanolic acid (red) showed stable RMSD between 4 ns to 20 ns at RMSD range between 0.21 nm to 0.38 nm, after 20 ns RMSD increases and constant at 0.41 nm (Figure 3). The differences of backbone RMSD with carvacrol, ursolic acid and oleanolic acid complex in protein backbone RMSD suggest that Mpro have conformational changes in oleanolic acid while no significant change shown in carvacrol and ursolic acid during 50 ns simulation.

Figure 3

RMSD analysis of protein backbone during MD simulation.

The compound carvacrol (red) has shown fluctuations at two different intervals on 50 ns time scale as it may has changed a conformation in binding region of Mpro. The first stable conformation between 6 ns to 29 ns and second stable conformation is between 34 ns to 50 ns. The RMSD remain constant at 0.3 nm and large fluctuation observed between 28 ns to 33 ns RMSD > 0.4 nm. As large fluctuation in compound, no effect on protein structure was found. Similar changes also seen in ursolic acid (magenta) RMSD fluctuations, it shown stable RMSD after 18 ns between 0.58 nm to 0.66 nm while in starting period between 4 ns to 17 ns, RMSD consistently fluctuated (Figure 4). It also not affects the fluctuation of protein backbone which depicts that binding region of protein has some fluctuation that causes the compound fluctuations during simulation. It may cause by large binding region and presence of loop at binding region.

Figure 4.

Analysis of RMSD of ligands and Mpro complexes; Mpro-carvacrol complex (Black), Mpro-oleanolic acid complex (Red), Mpro-ursolic acid complex (Blue).

The oleanolic acid (blue) showed higher but stable RMSD between 0.44 nm to 0.61 nm throughout simulation. Oleanolic acid binding affects the RMSD of protein backbone but it has stable till end of the simulation. Two similar conformations of ursolic and oleanolic acid showed different behaviour with protein binding region during simulation analysis. This observation we also confirm by their local changes in residues level by RMSF plot. Protein backbone with carvacrol (black) and ursolic acid (blue ) not affect the behaviour of protein or structure of protein during simulation while oleanolic acid (red) affects the structure of protein during simulation. RMSF confirms the changes in structure with binding of these phytochemical compounds (Figure 5).

Figure 5.

Analysis of RMSF of Cα during MD simulation; Mpro-carvacrol complex (Black), Mpro-oleanolic acid complex (Red), Mpro-ursolic acid complex (Blue).

Hydrogen bond analysis

Hydrogen bonding plays a crucial role in determining the binding strength in between ligands and protein. The carvacrol (red) and ursolic acid (black) have constant range of hydrogen bond between 2 to 4 in whole simulation while oleanolic acid showed changes in bonding. More hydrogen (>5) bonds between 0 to 12 ns, after 12 ns the hydrogen bond are <4 and last 10 ns, the hydrogen bond between 2 to 3. This may suggests that there is a conformational change around oleanolic acid in the binding site during simulation (Figure 6). Over all observation suggested that all three protein-ligand complexes are stable during simulation.

Figure 6.

Intermolecular hydrogen bonds between the ligands and Mpro protein; Mpro-carvacrol complex (Black), Mpro-oleanolic acid complex (Red), Mpro-ursolic acid complex (Blue).

Free binding energy analysis/Poisson − Boltzmann surface area (MM-PBSA)

The pre docking suggest the binding energy of complex while binding free energy (ΔG bind) analysis after simulation suggest the consistency of non-bonding interactions energy of binding region with compound also called post docking. In previous reports, it has been shown that 100–200 snapshots are enough to calculate the binding free energies (Khan et al., 2020; Sarma et al., 2020). The van der waal energy component compared with each complex, the carvacrol (−93.913 +/−9.776) has less effect on binding affinity while ursolic (−189.889 +/−12.027) and oleanolic acid (−197.509 +/−11.086) has strong binding affinity. Similar results are shown in binding energy, ursolic (−168.918 +/−13.703) and oleanolic acid (−158.999 +/−15.306) showed similar effect of binding energy while moderate in carvacrol (−82.781 +/−9.401). The electrostatic energy doesn’t show significant in ursolic and oleanolic acid while moderate in carvacrol. The polar solvation and SASA energy are shown moderate effects on binding energy component in each compound (Table 3).

Table 3.

Free binding energy calculations of each Mpro and ligand complexes.

Complexes	Binding energy (kJ/mol)	van der Waal energy (ΔEvdW)(kJ/mol)	Electrostattic energy (ΔEelec), (kJ/mol)	Polar solvation energy (ΔG polar) (kJ/mol)	SASA energy (kJ/mol)
Carvacrol	−82.781 +/− 9.401	−93.913 +/− 9.776	−37.778 +/− 6.738	58.079 +/− 6.712	−9.169 +/− 0.708
Oleanolic acid	−158.999 +/− 15.306	−197.509 +/− 11.086	−5.129 +/− 5.617	59.634 +/− 8.719	−15.995 +/− 0.944
Ursolic acid	−168.918 +/− 13.703	−189.889 +/− 12.027	−1.985 +/− 2.485	37.781 +/− 9.198	−14.825 +/− 1.721

Radius of gyration analysis

In order to determine the compactness of the system with the time, Rg was calculated, in which higher Rg values depict less compactness (more unfolded) with conformational entropy, while low Rg values explain high compactness with more stability in the structure (more folded). As evident from Figure 7(a), the simulation Rg values of all three cases reported as 2.05–2.15 nm. The fewer changes in the Rg value exhibited the stability of the protein in the complex, which showed less difference. Rg results revealed that the binding of these three molecules does not induce structural changes. The Rg values of all three protein-ligand complexes (Figure 7(a)) support their condensed architecture as well as size. In a recent study, Khan et al. (2020) calculated the Rg values for Remdesivir Saquinavir and Darunavir drugs (reported as inhibiting drugs against chymotrypsin-like protease of SARS-CoV-2). The average gyration score of Remdesivir was found to be 22.25 ± 0.1 Å, while Saquinavir and Darunavir durgs prompted to have average Rg scores of 22.32 ± 0.4 Å and 22.29 ± 0.6 Å, respectively.

Figure 7.

(a) Radius of gyration (Rg); Mpro-carvacrol complex (Black), Mpro-oleanolic acid complex (Blue), Mpro-ursolic acid complex (Red); (b) Solvent surface accessible area (SASA); Mpro-carvacrol complex (Black), Mpro-oleanolic acid complex (Red), Mpro-ursolic acid complex (Blue).

SASA analysis

Additionally, we also conducted the solvent accessible surface area (SASA) analysis of all of the three protein-ligand complexes. SASA is important to measure of the receptor exposed to the solvents during the simulation. The exposures of the hydrophobic residues by the binding of the ligand molecules contribute to the values of the SASA. The SASA value exhibited between 155–165 nm2 (Figure 7(b)) showed that the binding of ligands does not affect the folding of the protein. Khan et al. (2020) have performed the SASA analysis for 3CLpro and its protein-drug complexes. They reported that SASA values for the 3CLpro, 3CLpro-Paritaprevir complex and 3CLpro-Raltegravir complex were 133.89 nm2, 131.79 nm2 and 131.26 nm2, respectively. The results reported in the present study suggest that all three of the complexes were impressively stable after the binding of ligands to active sites of Mpro protein.

Conclusion

The Mpro protein has shown to be crucial and highly potent target for the inhibition of novel COVID-19. This study concludes three natural compounds (ursolic acid, carvacrol and oleanolic acid) as potential inhibitors against Mpro. Molecular docking analysis revealed that Carvacrol molecule having less binding energy as compared to other two molecules, oleanolic acid and ursolic acid. The binding mode of interaction was found to be reasonably good. MD simulations revealed that all three docking complexes showed stability at 50 ns. These inhibitors also fulfil the ADME parameters as well as Lipinski’s rule of five. All these reported compounds are natural and also commercially available for further in vivo/in vitro validations. The information generated from this present study may be utilized in future for the development of more phytochemical based therapeutics against COVID-19.