Molecular dynamics simulation perception study of the binding affinity performance for main protease of SARS-CoV-2.

Like common cold and flu, SARC-CoV-2 virus spreads by droplets of sneezes or coughs which virus affects people of various age groups. Today, this virus is almost distributed all over the world. Since binding process plays a crucial role between host and receptor, therefore, we studied the molecules intended toward inhibition process through molecular docking and molecular dynamics simulation process. From the molecular docking study, it is noteworthy that remdesivir shows better binding affinity toward the main protease of SARS-CoV2 compared to other studied drugs. Within studied phytochemicals, carnosic acid shows better binding poses toward main protease of SARS-CoV2 among studied phytochemicals. The amino acid residues GLN110 and PHE294 were almost found in all the studied interactions of drugs and phytochemicals with main protease of SARS-CoV-2. Furthermore, the results show a larger contribution of the Van der Waals energies as compared to others like electrostatic energies suggesting that ligands at the binding pocket are predominantly stabilized by hydrophobic interactions. The conformational change during ligand binding was predicted from Gibbs free energy landscape analysis through molecular dynamics simulation. We observed that, there were two main free energy basins for both docked carnosic acid complex and for docked remdesivir complex, only one main free energy basin was found in the global free energy minimum region.Communicated by Ramaswamy H. Sarma.

Introduction

Severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2) formerly referred to as novel coronavirus-19 (nCoV-19). World Health Organization (WHO) has announced that coronary virus disease 2019 (COVID-19) is a pandemic. SARS-CoV-2 causes medical problems such as pneumonia, sore throat, headache, nausea, runny nose, diarrhea, severe respiratory illness, multi-organ failure and occasionally death. Like common cold and flu this virus spreads by droplets when a person sneezes or coughs. This virus affects people of any age groups. Today, this virus is almost distributed all over the world. Fever, cough, headache and shortness of breath are the symptoms of this virus (Cheung et al., 2020). Severe symptoms such as bronchitis or pneumonia can be present in people with very weekly immunity (Zhao et al., 2020). Corona virus reaches the cell through the enzyme-converting angiotensin II (ACE2) spike receptor binding domain (Spike-RBD) of surface spike glycoprotein (S protein) (Mizuiri, 2015). ACE2 exists in epithelial cells of lungs, liver, blood vessels and heart. The ACE2 expresses itself deeply in cardiovascular and kidney systems (Yan et al., 2020). Subsequently, the level of ACE2 expression in patients with diabetic nephropathy or complex renal diseases is low. However, specific angiotensin-converting enzyme (ACE) inhibitors, thiazolidinediones and various steroids that simulate or boost the level of ACE2 expression are treated in patients with diabetic and renal disease (Mizuiri, 2015). Conversely, the ACE inhibitors and angiotensin receptor antagonists (ARB) are the most prescribed antihypertensive medications. In particular, ACE2 's increased level of expression will facilitate infection with a novel corona virus disease (COVID 19) (Guo et al., 2020). Therefore, most patients with diabetic nephropathy, renal failure and hypertension experience serious problems with COVID-19.

Using structure-based virtual screening, Elmezayen et al. (2020) reported drugs, such as Talampicillin, Lurasidone, Rubitecan and Loprazolam for COVID-19. Recently identified polyphenols from green tea by Ghosh et al. (2020) can be used as potential inhibitors against SARS CoV-2 Mpro and are promising drug candidates for COVID-19 therapy.

COVID-19 patients are treated with a number of antiviral drugs including ritonavir, umifenovir, lopinavir, ribavirin, darunavir, chloroquine and hydroxychloroquine (Dong et al., 2020). SARS-CoV-2 binds the spike receptor binding domain (Spike-RBD) to the target cells of the host via angiotensin converting enzyme 2 (ACE2) (Yan et al., 2020). The positive strand single-stranded RNA (ssRNA) joins the cell after binding to the host cell, which is a genetic material of SARS-CoV-2. This can function as both messenger RNA (mRNA) and genome. Host ribosome can directly transform the viral ssRNA inside the host cell into protein. The first protein which serves as functions for genome replication is expressed after infection. The positive strand of the viral genome to viral replication complexes formed in conjunction with viral protease. The complex includes proteins of both viral proteins and proteins of origin from host cells. The RNA-dependent RNA polymerase or RNA replicase, which is an enzyme, encodes all virus genome positive strand ssRNA, on the other hand. This will help catalyze RNA replication from the RNA template. RNAs and viral proteins are then combined into a new virus, and this process goes on until some agent destroys the virus cycle (Zhu et al., 2013).

In vitro studies suggested that chloroquine shows different influenza viruses (Yan et al., 2013), HIV (Boelaert et al., 2001; Romanelli et al., 2004), dengue virus (Farias et al., 2014), Zika virus (Delvecchio et al., 2016), and it is closely related to SARS-CoV-1 as a versatile drug molecule. The Endo-lysosomal pathway for entry of SARS-CoV-1 cells (Burkard et al., 2014). Previous work indicates chloroquine may be an effective therapeutic action against SARS-CoV-2 (Devaux et al., 2020). The antiviral treatment of COVID-19 such as lopinavir and ritonavir combination drug therapy, as these medications suppresses the HIV protease. Protease is a viral enzyme that allows, when forming a new virus chain, to cleave the long protein chain within the host. Previous research (Keshtkar-Jahromi & Bavari, 2020) has indicated that, in conjunction with hydroxychloroquine, lopinavir-ritonavir has an inhibitory effect against SARS-CoV-2. The chloroquine also supports the antiviral activity by interfering with the virus fusion cycle by increasing the pH and even altering the cellular receptor glycosylation mechanism in the case of corona viruses (Savarino et al., 2003). Chloroquine can interact with lopinavir-ritonavir, which may result in delayed ventricular re-polarization, which means QT interval prolongation (assesses the heart’s electronics properties). However, hydroxychloroquine is less toxic and more soluble than chloroquine (Plantone & Koudriavtseva, 2018), and it is recommended for the treatment of COVID-19 patients instead of chloroquine (Sahraei et al., 2020). Umifenovir (Arbidol) is a commonly used antiviral drug that inhibits virus fusion with host cell and also inhibits viral genome replication (Zhai et al., 2020). Interferon alpha is mainly used as an antiviral medication for the treatment of hepatitis and also prevents the SARS-CoV replication cycle (Dong et al., 2020; Stockman et al., 2006). In comparison, darunavir (Prozekal) is an antiviral drug that is used for treating AIDS. The darunavir primarily inhibits the viral protease in order to immature the virus or to stop new virus synthesis in the host cells. So this drug may be a possible COVID-19 inhibitor (Zhu et al., 2020). The risk of both acute respiratory distress syndrome (ARDS) and SARS is minimized by ribavirin (Chu et al., 2004; Yuan et al., 2020). A commonly used antiparasitic drug is from WHO's list of essential medicines niclosamide. It used to treat infections of humans with tapeworms (Andrews et al., 1982). It is anti-cestodal agent against tapeworm infection by inhibiting or inducing oxidative phosphorylation processes of adenosine triphosphate production in mitochondria and also controls multiple biological pathways (Alasadi et al., 2018). Niclosamide has an inhibitory activity on SARS-CoV-2 and viral replication of Middle East Respiratory Syndrome (MERS-COV) (Gassen et al., 2019; Wen et al., 2007).It can be used for potential therapeutic treatment against COVID-19 together with its antiviral properties (Xu et al., 2020). Caly et al. (2020) reported in a recent in vitro analysis that ivermectin is a SARS-CoV-2 inhibitor. Previously ivermectin was used as an antiparasitic drug. The ivermectin effects on ∼5000 fold reduction in SARS-CoV-2 viral RNA within 48 h were identified from the in vitro analysis, thus inhibiting viral replication. It's anti-SARS-CoV-2 antiviral properties may be used to treat COVID-19 (Caly et al., 2020). Patients with positive COVID-19 are treated with different combinations of antiviral drugs in the absence of any specific treatment regimen (Yousefifard et al., 2020).

Apigenin is a flavonoid that has antiviral activity against viral infection in the diet. The interaction of viral RNA with the host is disrupted and the host factors are targeted to stop viral replication (Zhang et al., 2014). Drymaritin is an active compound of the drymaria diandra plant, a type of alkaloid. It displays anti-HIV activity, and the activity was previously investigated in the cell line of H9 lymphocytes. The HIV virus has been infected with cells of the H9 lymphocytes and causes replication of HIV (Hsieh et al., 2004). A natural compound is carnosic acid, which demonstrates anti-viral activity against viral infection. This natural compound should be healthy for humans to eat. Carnosic acid significantly inhibits the human respiratory syncytial virus and is a possible therapeutic compound against it (Shin et al., 2013). It effectively defends the expression of viral replication and viral genes without affecting cell viability. It also shows activity against anti-cancer and antioxidants. Isoobtusitin is a natural compound which has shown the inhibitory effect on poliovirus and substituted the viral replication within the host. But isoobtusitin prevents just about 10 percent of viral replication in the case of HIV viruses (Curini et al., 2006). Ellagic acid is a flavonoid molecule occurring naturally. This demonstrates anti-hepatitis B virus activity. Ellagic acid inhibits hepatitis B viral protein (HBeAg) secretion in HepG 2.2.2.15 cells. Ellagic acid is a possible candidate for treatment against hepatitis B virus (Shin et al., 2005). Morin is a natural flavonoid group compound. This demonstrates an antiviral activity against infection with the herpes simplex virus. Within the host, Morin prevents the viral replication of herpes simplex (Bunyapraphatsara et al., 2000). Scutellarein is a natural flavonoid-group compound. It is potential SARS-CoV inhibitor. It interferes with SARS-CoV helicase protein's function as an ATPase. It has not been shown to exhibit any cytotoxicity against normal epithelial cells in in vitro examinations (Yu et al., 2012). In active phytochemicals like naringenin, morusin, oxymatrine, honokiol (Leung et al., 2020) were reported against COVID-19 and Mani et al. reported quercetin 7-rhamnoside, myricetin, tryptanthrin natural product derived phytochemicals as potential inhibitor against coronavirus (Mani et al., 2020).

Currently, COVID-19 is not available for immediate care or diagnosis. The purpose of this study is to determine the binding energies of different antiviral drugs and phytochemicals with main protease of SARS-CoV-2. Here we, taken some food and drug administration (FDA) approved drugs darunavir, hydroxychloroquine, interferon alpha, lopinavir, niclosamide, ribavirin, ritonavir, umifenovir, ivermectin and remdesivir. We've also chosen phytochemicals like apigenin, drymaritin, carnosic acid, isoobtusitin, ellagic acid, morin, scutellarein and triterpenoids. We performed molecular dynamics simulations (100 ns production run) for studying the dynamical behavior of better conformation.

Protein and ligand preparation protocol

The three-dimensional structure of main protease was retrieved from the protein data bank database with PDB ID 6Y2E (Kumar et al., 2020; Zhang et al., 2020). We used AutoDockTools (ADT) in our in silico analysis to prepare receptors and ligands. In the docking process (Bissoyi et al., 2017; Chand et al., 2018; Panda et al., 2019; Sahu et al., 2019; Sahu & Pattanayak, 2019), one molecule was taken as a receptor and another considered as a ligand, where we took a grid box for receptor. In this study the main protease protein's three-dimensional structure was taken as a receptor and various drugs (darunavir, hydroxychloroquine, interferon alpha, lopinavir, niclosamide, ribavirin, ritonavir, umifenovir, ivermectin and remdesivir) and phytochemicals (apigenin, drymaritin, carnosic acid, isoobtusitin, ellagic acid, morin, scutellarein and triterpenoids) are considered as ligand. During the docking we used gradient optimization algorithm to perform docking process (Meng et al., 2011). For the study of docking, the three-dimensional structure of main protease protein was considered. The non-essential water molecules including the hetero atoms were extracted before continuing with the docking test. Once all the poor touch was removed, we attached Kollman charges and hydrogen atoms to the receptor’s polar contacts. Since SARS-CoV-2 main protease(PDB ID: 6Y2E) proteins perform their functions through interactions with phytochemicals and drugs, the correct identification of the protein-ligand binding site plays an important role in the process of drug discovery. We studied SARS-CoV-2 main protease through CASTp 3.0 (Computed Atlas of Surface Topology of Proteins) online software for the prediction of binding pockets (Tian et al., 2018). It is noted that the alpha shape principle of computational geometry was applied by CASTp to detect and measure pockets. We found a total number of 36 binding pocket numbers in the main protease of SARS-CoV-2. We only selected only 10 best binding pockets based on the surface area and volume of the pockets. The grid box size 126 × 126 × 126 and the centers were −16.52, −26.112, 17.524 in the direction of x, y, z, respectively, was selected for 6Y2E in the docking process. The better binding energies obtained from docking process are found in the 2nd binding pocket (Pocket ID of CASTp result), that covering the residues GLN107, PRO108, GLY109, GLN110, PRO132, ILE200, THR201, VAL202, ASN203, GLU240, PRO241, ASP245, HIS246, ILE249, THR292 and PRO293. The above said all antiviral drug molecules and phytocompounds are downloaded from Pubchem database. The three-dimensional structures have been downloaded, and then converted to format .pdbqt using AutoDockTool. We also tested the drug likeness for the studied phytochemicals by computationally using the Swiss ADME server (Daina et al., 2017) and DataWarrior software package (Sander et al., 2015). We used Discovery Studio Visualizing tools for visualization and binding site analysis (Visualizer, D.S., 2013; Version 2.0. 1.7347, Accelrys Software Inc., San Diego).We performed molecular dynamics simulations (Al-Obaidi et al., 2020; Elmezayen et al., 2020) using the Gromacs program package (Berendsen et al., 1995; Van Der Spoel et al., 2005) to evaluate the dynamical properties and binding efficiency of the above listed complexes. We eliminate all conflicts from the studied protein and docked complexes to form the simulations method. All the simulations systems analyzed were used using all-atom force field CHARMM 36 (Best et al., 2012) for 100 ns trajectory. In the cubic shell, with simple point charge (SPC/E) water model, the three systems were solved (Berendsen et al., 1987). To maintain neutral ionization, we added four Na+ ions. To avoid the edge effect of periodic boundary conditions was employed of our systems. In order to minimize protein energy and docked complexes, the steepest descent algorithm was used and the long-range electrostatic force has been served using the Particle Mesh Ewald (PME) method (Ewald, 1921). After minimization, three additional steps were applied during simulation, such as heating, equilibration and production process. Constant number, volume and temperature (NVT) were introduced at 300 K and 1.0 atm followed by constant number of particles, pressure and temperature (NPT) for 100 ps. Implementation of the LINCS algorithm (Hess et al., 1997) has regulated all the covalent bonds. Pressure coupling of Parrinello-Rahman (Martonák et al., 2003) and temperature coupling (Berendsen et al., 1984) was maintained with constant pressure of 1 bar and 0.1ps, respectively. For each, 1.0 nm was chosen for the short-range cut-off of the columbic and Van der Waals interaction. During simulation, we were following a time step with 2 fs. For production simulations, all of the MD trajectories were saved per 10 ps with 100 ns.

Result and discussions

Molecular docking

AutoDock Vina, was used to assess the binding affinity of the phytochemicals and drugs to the main protease. On the basis of various interactions, binding energy was found. Here in this study, we carried out docking for darunavir, hydroxychloroquine, interferon alpha, lopinavir, niclosamide, ribavirin, ritonavir, umifenovir, ivermectin and remdesivir drugs and apigenin, drymaritin, carnosic acid, isoobtusitin, ellagic acid, morin, scutellarein and triterpenoids phytochemicals with main protease. Hydrogen bonding, hydrophobic and electrostatic interactions occurred in nearly all drugs and phytochemical cases. The protease protein's binding energy to the phytochemical, apigenin was −7.3 kcal/mol. THR111 and GLN110 residues include in the H-bonding with bond lengths 2.12 to 4.08 Å. The amino acid residues ILE249 and PRO293 interact through hydrophobic contact with bond lengths 4.77 and 5.29 Å, respectively. Protease protein binding energy to the phytocompound drymaritin is −6.6 kcal/mol. The residue GLN110 was involved in the H-bonding with a bond length of 3.85 Å. The main protease's amino acid residue PHE294 has established hydrophobic contacts with varying bond lengths between 3.85 and 5.07 Å. Diagrammatic representation of docking interaction of apigenin and drymaritin with main protease of SARS-CoV-2 shown in Figure 1.

Figure 1.

Diagrammatic representation of docking interaction of (a) apigenin (b) drymaritin with main protease of SARS-CoV-2. The ball and stick and stick models represented to ligands and interact amino acid residues, respectively. The pink dotted and green dotted lines represented the hydrophobic contacts and hydrogen bonds, respectively.

Binding energy of the main protease with the phytocompound carnosic acid was −7.9 Å. The residue, GLN110 has included in H-bonding with a bond length of 3.86 Å. Amino acid residue PHE294 contains the hydrophobic contacts with bond lengths ranging from 3.95 to 4.88 Å. The amino acid residue PRO293 was involved in hydrophobic contact with bond length of 4.28444 Å. Protease protein binding energy to phytochemical isoobtusitin was −5.8 kcal/mol. The residue ARG40 is linked in H-bonding with bond lengths from 2.89726 to 3.31188 Å, whereas ARG188 is involved in hydrophobic contacts with bond length of 4.10143 Å. Amino acid residue CYS85 of the main protease has participated in hydrophobic contacts with varying bond lengths of 4.48 to 5.37 Å. The amino acid residue TYR54 of the main protease has participated in hydrophobic contacts with varying bond lengths of 5.01 to 5.17 Å. Schematic representation molecular interactions of carnosic acid and isoobtusitin with main protease of SARS-CoV-2 shown in Figure 2.

Figure 2.

Illustration of docking interaction of (a) carnosic acid (b) isoobtusitin with main protease of SARS-CoV-2. The ball and stick and stick models are represented to ligands and interact amino acid residues, respectively. The pink dotted and green dotted lines are represented the hydrophobic contacts and hydrogen bonds, respectively.

The binding energy of protease protein with the phytochemical, ellagic acid was found −7.8 kcal/mol. The amino acid residue GLN189 of protease protein has participated through H-bonding with bond length of 2.66 Å. Amino acid residue LEU141 was accompanied by H-bonding with varying bond lengths of 2.29 to 2.51 Å, whereas SER144 was also accompanied with H-bonding with variable bond lengths in the rage of 2.18 to 3.07 Å. GLY143 was participated through by H-bonding with bond length of 3.24 Å. CYS145 was accompanied by H-bonding with varying bond lengths of 3.26 to 4.06 Å. The residue, GLU166 was participated in H-bonding with bond lengths of 3.18 Å. Further, the amino acid residue of the main protease, ASN142 was also bonded through H-bonding with bond length of 3.57 Å. The binding energy of protease protein with the phytochemical, morin was −6.7 kcal/mol. THR25 is participated in H-bonding interaction with a bond length of 2.14 Å. Amino acid residue CYS44 is participated in H-bonding interaction with bond length of 2.91 Å where as CYS145 residue is accompanied in hydrophobic contact with bond length of 5.13 Å. The amino acid residues, CYS44, GLU166 and SER144 have participated through H-bonding with varying bond lengths of 2.91, 2.78 and 3.08 Å. and MET49 has participated through hydrophobic contact with bond length of 5.36 Å. The schematic representation is shown in Figure 3 of the molecular interactions of ellagic acid and morin with the main SARS-CoV-2 protease.

Figure 3.

Representation of docking interaction of (a) ellagic acid (b) morin with main protease of SARS-CoV-2. The ball and stick and stick models are represented to ligands and interact amino acid residues, respectively. The pink dotted and green dotted lines are represented the hydrophobic contacts and hydrogen bonds, respectively.

The amino acid residue of protease protein THR111 was participated through H-bonding with a bond length of 2.76 Å during interaction of scutellarein with main protease. Similarly the amino acid residue, ASN151 was bonded through H-bonding with variable bond lengths of 2.03 to 2.56 Å. Furthermore, the residue, GLN110 was accompanied through H-bonding with variable bond lengths in the range of 3.78 to 4.10 Å. The residue PHE294 was interacted with hydrophobic contact with bond lengths of 5.05 to 5.67 Å. The residues ILE249 and PRO293 have interacted through hydrophobic contact with bond lengths of 4.73 and 5.27 Å, respectively. The schematic representation of binding conformation was shown in Figure 4(a). The binding energy of protease protein with the phytochemical, scutellarein was found −7.5 kcal/mol. The binding energy of protease protein with the phytochemical, triterpenoids was found −7.8 kcal/mol. The amino acid residues of protease protein such as ASP245 and GLN110 accompanied via H-bonding with bond lengths of 2.39 and 3.09 Å, respectively. The residue ILE249 was accompanied by hydrophobic contact with varying bond lengths in the range of 4.05 to 5.11 Å. The VAL202 residue was followed by hydrophobic interaction with a bond length of 4.39 Å where the PRO293 residue interacts with a bond length of 5.40 Å. The amino acid residue HIS246 was interacted through hydrophobic contact with bond lengths in the range of 4.83 to 5.15 Å with different atom of triterpenoids. The molecular interactions were shown in Figure 4(b).

Figure 4.

Interpretation of docking interaction of (a) scutellarein (b) triterpenoids with main protease of SARS-CoV-2. The ball and stick and stick models represented to ligands and interact amino acid residues, respectively. The pink dotted and green dotted lines represented the hydrophobic contacts and hydrogen bonds, respectively.

The amino acid residues of protease protein such as CYS145, HIS164, THR25 and GLU166 were accompanied through H-bonding with bond lengths of 2.56, 2.82, 3.18 and 3.05 Å, respectively. The amino acid residue MET49 was interacted through hydrophobic contact with bond length of 4.84 Å. The binding energy of protease protein with the drug, darunavir was found −7.0 kcal/mol. The binding energy of protease protein with the drug hydroxychloroquine was found −5.1 kcal/mol. The amino acid residues of protease protein such as GLU288, LYS5 and GLN127 were interacted through H-bonding with bond lengths of 2.21, 3.26 and 3.77 Å, respectively. The amino acid residue, PHE291 was accompanied with hydrophobic contacts with variable bond lengths of 5.16 to 5.19 Å where as PHE3 was interacted with a bond length of 4.59 Å. The amino acid residue, LEU282 was accompanied with hydrophobic contacts with bond length of 3.72 Å. Hydrophobic contacts with variable bond lengths of 4.15 to 5.18 Å were followed by the amino acid residue, LYS5. TRP207 is accompanied through hydrophobic contacts of 4.79 Å. Diagrammatic representation of docking interaction of darunavir and hydroxychloroquine with main protease of SARS-CoV-2 shown in Figure 5.

Figure 5.

Diagrammatic representation of docking interaction of (a) darunavir (b) hydroxychloroquine with main protease of SARS-CoV-2. The ball and stick and stick models are represented to ligands and interact amino acid residues, respectively. The pink dotted and green dotted lines are represented the hydrophobic contacts and hydrogen bonds, respectively.

The binding energy of protease protein with the drug, IFN-alpha was found −6.9 kcal/mol. The amino acid residues of protease protein such as SER158, ASP153 and GLN110 were interacted through H-bonding with bond lengths of 2.63, 3.43 and 3.71 Å, respectively. The amino acid residue, PHE294 was interacted through hydrophobic contact with bond length of 5.26 Å. Amino acid residue VAL104 was interacted through hydrophobic contact with bond length of 4.81 Å. The binding energy of protease protein with the drug, lopinavir was found −7.5 kcal/mol. The amino acid residue GLN110 interacted through H-bonding with bond lengths of 2.80 to 3.16 Å. The amino acid residues, ILE249, PHE294, HIS246, ILE249, ILE200, VAL202, VAL104, ILE106 and HIS246 are accompanied with hydrophobic contact with variable bond lengths of 3.55, 3.69, 5.30, 4.68, 5.18, 5.01, 4.94, 4.38, 5.39 and 5.30 Å, respectively. The docking interaction of ifn-alpha and lopinavir with the main SARS-CoV-2 protease was shown in Figure 6.

Figure 6.

Representation of docking interaction of (a) ifn-alpha (b) lopinavir with main protease of SARS-CoV-2. The ball and stick and stick models represented to ligands and interact amino acid residues, respectively. The pink dotted and green dotted lines represented the hydrophobic contacts and hydrogen bonds, respectively.

The binding energy of protease protein with the drug, niclosamide was −7.0 kcal/mol. The amino acid residue GLN110 interacted through H-bonding with bond length 3.07 Å. The amino acid residue PHE294 was associated with hydrophobic contact with a bond length of 4.47 Å. The amino acid residues ILE200 and ILE249 accompanied with hydrophobic contact with variable bond lengths of 5.46 and 4.98 Å, respectively. Amino acid residue, VAL202 was accompanied with hydrophobic contact with variable bond lengths in the range of 4.80 to 5.48 Å. Residues of amino acid PRO293 followed by hydrophobic contact with a bond length of 4.71 Å. The protease binding energy of the ribavirin drug was found to be −6.0 kcal/mol. The PHE140, GLU166, SER144, GLU166 and ASN142 amino acid residues were followed by H-bonding with bond lengths of 1.93, 3.02, 2.96, 3.36 and 3.61 Å, respectively. The residue, CYS145 was accompanied by hydrophobic contact with bond length of 4.78 Å. Docking interaction of niclosamide and ribavirin with main protease of SARS-CoV-2 shown in Figure 7.

Figure 7.

Illustration of docking interaction of (a) niclosamide (b) ribavirin with main protease of SARS-CoV-2. The ball and stick and stick models are represented to ligands and interact amino acid residues, respectively. The pink dotted and green dotted lines are represented the hydrophobic contacts and hydrogen bonds, respectively.

The binding energy of protease with the drug, ritonavir was found to be −7.3 kcal/mol. ASP197 was accompanied with electrostatic interaction with bond length of 5.17 Å .It was accompanied with H-bonding with a bond length of 3.67 Å. THR199 was accompanied with H-bonding interactions with bond length in the range of 2.84 to 3.13 Å. TYR239 was participated in H-bonding with a bond length of 2.90 Å. Hydrophobic interaction was followed by LEU286 and LEU287 with bond lengths of 4.72 and 3.51 Å, respectively. The residue, VAL171 was accompanied in hydrophobic contact in the rage of 3.68 to 4.72 Å. ALA194 and LYS137 were accompanied through hydrophobic contact with bond lengths of 4.86 and 4.29 Å, respectively. The binding energy of protease protein with the drug, umifenovir was found to be −6.3 kcal/mol. GLN110 was accompanied with H-bonding with bond length of 2.87 Å. VAL202, ILE249 and PRO293 were accompanied with hydrophobic contact with variable bond lengths of 5.47, 5.08 and 4.77 Å, respectively. The molecular interactions of ritonavir and umifenovir with main protease of SARS-CoV-2 was shown in Figure 8.

Figure 8.

Diagrammatic depiction of docking interaction of (a) ritonavir (b) umifenovir with main protease of SARS-CoV-2. The ball and stick and stick models represented to ligands and interact amino acid residues, respectively. The pink dotted and green dotted lines represented the hydrophobic contacts and hydrogen bonds, respectively.

The binding energy of protease protein with the drug, remdesivir was found to be −7.9 kcal/mol. The residues, SER158, ASP295, GLN110, THR292, PHE294, ARG298 and ASN151 were accompanied with H-bonding with bond lengths of 2.62, 2.73, 3.25, 2.93, 3.32, 3.25 and 4.04 Å, respectively. ILE200, VAL104, ILE249 and HIS246 were accompanied with hydrophobic contact with variable bond lengths of 4.62, 5.46, 4.75 and 4.77 Å, respectively. The binding energy of protease protein with the drug, ivermectin was found to be −7.2 kcal/mol. GLN110 and ARG105 were accompanied with H-bonding with bond length of 3.36 and 3.39 Å, respectively. PHE294 was associated with hydrophobic contact with variable bond length in the range of 3.52 to 5.16 Å. The amino acid residue, ILE249, was accompanied with hydrophobic contact with variable bond lengths of 5.10 to 5.11 Å. VAL104 participated through hydrophobic contact with bond length of 3.63 Å. The schematic representation of remdesivir and ivermectin with main protease of SARS-CoV-2 shown in Figure 9.

Figure 9.

Illustration of docking interaction of (a) remdesivir (b) ivermectin with main protease of SARS-CoV-2. The ball and stick and stick models represented to ligands and interact amino acid residues, respectively. The pink dotted and green dotted lines represented the hydrophobic contacts and hydrogen bonds, respectively.

From the molecular docking study, it is noteworthy that remdesivir shows better binding affinity toward the main protease of SARS-CoV2 compared to other studied drugs. Within studied phytochemicals, carnosic acid shows better binding poses toward main protease of SARS-CoV2 among studied phytochemicals. The amino acid residues GLN110 and PHE294 were almost found in all the studied interactions of drugs and phytochemicals with main protease of SARS-CoV-2. Diagrammatic representation of docking interaction of honokiol, morusin, myricetin and naringenin with main protease of SARS-CoV-2 is shown in Figures S1–S4, respectively. For oxymatrine, quercetin 7-rhamnoside and tryptanthrin was shown in Figures S5–S7, respectively. GLN 110 residue found in interaction between apigenin, drymaritin, carnosic acid, scutellarein, tryptanthrin, lopinavir, remdesivir and ivermectin with enzyme. ILE249 common residue found in interaction between apigenin, scutellarein, lopinavir, remdesivir, niclosamide and ivermectin with enzyme. PRO293 residue was found in apigenin, carnosic acid and scutellarein. Both PHE 294 and GLN 110 residues were found both complex of remdesivir and carnosic acid with enzyme. Binding energy and different residues involved in the molecular interactions of phytochemicals and drugs with main protease of SARS-CoV-2 are given in the Table 1.

Table 1.

Binding energy and different residues involved in the molecular interactions of phytochemicals and drugs with main protease of SARS-CoV-2.

Phytochemicals / Drugs	Chemical compounds	Binding energy (kcal/mol)	Binding residues
Phytochemicals	Apigenin	−7.3	GLN110, THR111, ILE249, PRO293
	Drymaritin	−6.6	GLN110, PHE294
	Carnosic acid	−7.9	GLN110, PRO293, PHE294
	Isoobtusitin	−5.8	ARG40, PRO52, TYR54 , CYS85, ARG188
	Ellagic_acid	−7.8	LEU141, ASN142, GLY143, SER144, CYS145, GLU166, GLN189
	Morin	−6.7	THR25, CYS44, MET49, SER144, CYS145, GLU166
	Scutellarein	−7.5	GLN110, THR111, ASN151, ILE249, PRO293, PHE294
	Triterpenoids	−7.8	GLN110, VAL202, ASP245, HIS246 , ILE249, PRO293
	Honokiol	−5.9	LEU282,SER284,LYS5, GLU288LYS5, PHE3, TRP207, PHE291
	Morusin	−7.5	LEU287, THR199, LEU286, LEU272,LYS137
	Myricetin	−7.2	THR25,CYS44, GLU166, HIS41, CYS145, MET49
	Naringenin	−7.0	THR25, GLY143, GLU166, CYS44HIS41, THR25, MET165
	Oxymatrine	−6.7	CYS145, MET49, MET165, HIS41
	Quercetin 7-rhamnoside	−7.5	THR25,GLU166, HIS164,MET165
	Tryptanthrin	−7.2	THR111,THR292 ,GLN110
Drugs	Darunavir	−7.0	THR25, MET49, CYS145, HIS164, GLU166
	Hydroxychloroquine	−5.1	PHE3, LYS5, GLN127, TRP207, LEU282GLU288, PHE291
	IFN_alpha	−6.9	VAL104, GLN110, ASP153, SER158, PHE294
	Lopinavir	−7.5	VAL104, ILE106, GLN110, ILE200, VAL202, HIS246, ILE249, PHE294
	Niclosamide	−7.0	GLN110, ILE200, VAL202, ILE249, PRO293, PHE294
	Ribavirin	−6.0	PHE140, ASN142, SER144, CYS145, GLU166
	Ritonavir	−7.3	LYS137, VAL171, ALA194, ASP197, THR199, ASN238, TYR239, LEU286, LEU287
	Umifenovir	−6.3	GLN110, VAL202, ILE249, PRO293, PHE294
	Remdesivir	−7.9	VAL104, GLN110, ASN151, SER158, ILE200, HIS246, ILE249, THR292, PHE294, ASP295, ARG298
	Ivermectin	−7.2	VAL104, ARG105, `GLN110, ILE249, PHE294

ADME (Absorption, Distribution, Metabolism and Excretion) drug likeness properties

The toxicological and physicochemical properties of above said phytocompounds were predicted by using OSIRIS DataWarrior V5.2.1, whereas drug likeness and pharmacokinetics properties were studied by Swiss ADME online server. The physiochemical and potential ADME properties commonly assessed with Lipinski’s rules. Accordingly we were studied the different properties of all the above said phytocompounds. The properties were analyzed for each of the phytocompounds in the terms of molecular weight, partition coefficient (Log P), water solubility coefficient (Log S), hydrogen bond donors, hydrogen bond acceptor, polar surface area (PSA) and number of rotatable bonds (RB). All of the tested molecules met the specification of drug likeness properties, which may be considered for candidates lead molecule to inhibit the main protease for SARS-CoV-2 shown in Table 2.

Table 2.

Predicted drug likeness properties of studied phytochemicals.

Name of phytochemicals	Pubchem Id	Molecular formula	Molecular weight (g/mol)	PSA (Ų)	No. of RB	Log P	Log S	HBA	HBD
Apigenin	5280443	C15H10O5	270.24	90.90	1	2.11	−3.94	5	3
Drymaritin	11687449	C15H10N2O2	250.25	43.60	1	2.30	−3.28	3	0
Carnosic acid	65126	C20H28O4	332.43	77.76	2	3.80	−5.03	4	3
Isoobtusitin	5482812	C15H16O5	276.28	68.90	4	2.70	−3.55	5	1
Ellagic acid	5281855	C14H6O8	302.19	141.34	0	1.00	−2.94	8	4
Morin	5281670	C15H10O7	302.24	131.36	1	1.20	−3.16	7	5
Scutellarein	5281697	C15H10O6	286.24	111.13	1	1.81	−3.79	6	4
Triterpenoids	71597391	C29H44O5	472.66	97.99	1	3.94	−5.19	5	4
Honokiol	72303	C18H18O2	266.33	40.46	5	4.19	−4.74	2	2
Morusin	5281671	C25H24O6	420.45	100.13	3	4.35	−6.11	6	3
Myricetin	5281672	C15H10O8	318.24	151.59	1	0.79	−3.01	8	6
Naringenin	932	C15H12O5	272.25	86.99	1	1.84	−3.49	5	3
Oxymatrine	114850	C15H24N2O2	264.36	49.74	0	0.90	−2.18	2	0
Quercetin7-rhamnoside	5748601	C21H20O11	448.38	190.28	3	0.62	−3.93	11	7
Tryptanthrin	73549	C15H8N2O2	248.24	51.96	0	2.16	−3.29	3	0

The molecular weight of all the studied phytocompounds are below the 500 Daltons. The Log P value (Vraka et al., 2017) is a parameter for absorption or express the hydrophobicity of a molecule, the higher hydrophobicity have more capability to penetrate the plasma membrane of cell. The value of water solubility coefficient or Log S is an important parameter to analyze pharmacokinetic behavior of lead molecule distribution and absorption. The Log S value was limits within −4.5 to −1 (Trapani et al., 2005). The PSA, which is related to absorption, is calculated to be less than 140 Å2. The number of RB is always less than 10, the minimum RB shows the best structure confirmation. All the above studied parameters for the above said phytocompounds were analyzed with Lipinski’s rules, which express the drug likeness of each of the phytocompounds. The limits of the Lipinski’s rules are considered as the amounts of hydrogen bond acceptors are should be smaller than 10, numbers of hydrogen bond donors are should be less than 5, molecular weight should be less than 500 Daltons and Log P value should be within limits 1 to 5. The quantity of drug to be absorbed by the gastrointestinal tract is determined by GI, here all the phytochemicals are highly absorbed by the gastrointestinal tract. Followed by we are studied the capability of phytocompound to cross the blood-brain-barrier (BBB). The BBB results are indicating that the drymaritin and isoobtusitin are capable to cross the blood-brain-barrier and these two phytocompounds are may also be taken for neurological disorders. But the other studied phytocompounds are not capable to cross the blood-brain-barrier, which are may be selected for the lead compounds. The capability of the phytocompounds to permit the skin is determined by Log Kp value. The most − ve value indicates the more skin permeability. Also we are studied the competence of the phytocompounds to cross through the biomolecules membranes by cytochrome P450 (CYP). The predicted results reveal that all the studied phytochemicals do not show any hazardous of tumorigenicity reproductive effects and irritant effect but apigenin and morin are showing mutagenicity in nature. From all the predicted pharmacokinetics and toxicological properties, overall data reveals that all the studied four phytocompounds are may be considered as lead compounds inhibit the tumorigenesis. But apigenin and morin were could not be considered as an inhibitor of SARS-COV-2 because these two compounds are showing high mutagenicity.

Molecular dynamic simulation

Molecular dynamic simulation was employed to this study to predict the conformational changes and flexibility of docked complexes and main protease of SARS-CoV-2. The dynamic behavior of best binding confirmations of drug (remdesivir) and phytocompound (carnosic acid) with main protease of SARS-CoV-2 were studied by using Gromacs suit. We studied properties such as root mean square deviation (RMSD), root mean square fluctuation (RMSF), gyration radius (Rg), solvent accessible surface area (SASA), number of hydrogen bond analyzes, principal component analysis (PCA), free energy landscape and Poisson-Boltzmann surface molecular mechanics (MM/PBSA).

The root mean square deviation (RMSD) was predicted by using ‘gmx rms’ utility .The RMSD was calculated for backbone atoms of SARS-CoV-2 main protease and above said docked complexes. The RMSD predictions were represents the stability of the compounds. It has been shown that the initial points of deviations of all the three systems are started from ∼ 0.09 nm. For main protease the maximum deviation was traced ∼ 0.26 nm at ∼ 43 ns and the average value is ∼ 0.19 nm. Similarly the maximum deviation for carnosic acid docked complex was noted ∼ 0.25 nm at ∼ 46 ns whereas the average value is ∼ 0.17 nm. Meanwhile the maximum deviation for remdesivir docked complex was ∼ 0.19 nm at ∼ 65 ns and the average value is ∼0.14 nm. The deviations for main protease, carnosic acid docked complex and remdesivir docked complex were ∼ 0.21, ∼ 0.2 and ∼ 0.15 nm at end of the trajectory. If we see the pattern of fluctuation, all the three systems gained stability after ∼ 60 ns (shown in Figure 10).

Figure 10.

RMSD plot of backbone atoms with respect to time in water of SARS-COV-2 main protease (red), docked complex of remdesivir (blue) and docked complex of carnosic acid (green).

Then, we extended our study to observe the fluctuation behavior of different residues by using gmx rmsf protocol of Gromacs module. The RMSF prediction represented the variation of protein stability. As shown in Figure 11, the mode of fluctuation is more, after residue number 225 in all the three studied simulations systems. Key residues in all three structures analyzed for detailed analysis from 1 to 200 residue numbers showed the same fluctuating pattern. Key residues are THR25, SER46, TYR154, GLY170 and ASP197, which were taken for RMSF prediction. The fluctuation pattern for key residues THR25, SER46, TYR154, GLY170 and ASP197 of main protease are at ∼ 0.12 nm, ∼ 0.16 nm, ∼ 0.12 nm, ∼ 0.16 nm and ∼ 0.14 nm, respectively. Similarly as per above said manner for carnosic acid docked complex the fluctuation picks are at ∼ 0.12 nm, ∼ 0.17 nm, ∼ 0.12 nm, ∼ 0.13 nm and ∼ 0.14 nm, whereas the fluctuation picks of remdesivir docked complex are at ∼ 0.13 nm, ∼0.16 nm, ∼ 0.14 nm, ∼ 0.17 nm and ∼0.16 nm. From the overall results we observed that, key residues THR25, SER46, TYR154, GLY170 and ASP197 were more prominent during interactions. RMSF trajectories provide enough knowledge about the complex's stability. SARS-CoV-2 main protease backbone atoms showing less fluctuation than predicted, docked remdesivir complex and docked carnosic acid complex.

Figure 11.

The RMSF plot of backbone of SARS-CoV-2 main protease (red), docked complex of remdesivir (blue) and docked complex of carnosic acid (green).

Radius of gyration properties assessed to predict protein compactness behavior. The initial point of Rg value for main protease, carnosic acid docked complex and remdesivir docked complex were ∼ 2.17 nm, ∼ 2.19 nm and ∼2.19 nm, respectively. For main protease the maximum Rg value was ∼ 2.23 nm at ∼ 35 ns and the average value was ∼ 2.17 nm. Similarly the maximum Rg value for carnosic acid docked complex was ∼ 2.4 nm at ∼ 47 ns and average value was ∼ 2.19 nm, whereas maximum Rg value for remdesivir docked complex was ∼ 2.25 nm at ∼ 78 ns with average value ∼ 2.22 nm. Both of the docked complexes are maintaining an equilibrium style of fluctuation throughout the simulation. From the overall Rg prediction (Figure 12), we observed that the remdesivir docked complex has more Rg value and it is more flexible as compared to carnosic acid docked complex and main protease.

Figure 12.

The time evolution plot of Rg of SARS-COV-2 main protease (red), docked complex of remdesivir (blue) and docked complex of carnosic acid (green).

In order to analyze more extensive study about compactness behavior of backbone atoms for all the three simulations systems were predicted by using solvent accessible surface area (SASA). The SASA predictions are employed by using ‘gmx sasa’ protocol of Gromacs module. The initial surface area occupied by main protease, carnosic acid docked complex and remdesivir docked complex are ∼ 166.73 nm2, ∼ 168.45 nm2 and ∼167.21 nm2, respectively. The average surface area occupied by main protease, carnosic acid docked complex and remdesivir docked complex are ∼ 168.45 nm2, ∼ 168.64 nm2 and ∼ 169.3 nm2, respectively. From the overall results we observed that main protease has occupied minimum surface area and remdesivir docked complex occupied maximum surface area. It is noted that lower value of solvent accessible surface area shows more compactness. In our present study the main protease is more compactness and remdesivir docked complex is more flexible as shown in Figure 13.

Figure 13.

SASA plot of backbone of SARS-COV-2 main protease (red), docked complex of remdesivir (blue) and docked complex of carnosic acid (green).

In order to predict the hydrogen bond networks, this affects the flexibility of the protein. The number of hydrogen bonds present in protein is predicted by using ‘gmx hbond’ utility of Gromacs module. The initial number of hydrogen predicted for main protease, carnosic acid docked complex and remdesivir docked complex are ∼ 202, ∼ 214, ∼204, respectively, at the initial point of the trajectory. The average numbers of hydrogen bonds are found for main protease, carnosic acid docked complex and remdesivir docked complex are ∼ 225, ∼ 222, ∼ 220, respectively. From our hydrogen bond analysis (Figure 14) we observed that, up to 40 ns docked complex of carnosic acid was showing higher number of hydrogen bonding. Within 60 to 100 ns main protease was showing highest number of hydrogen bonding as compare to carnosic acid docked complex and remdesivir docked complex.

Figure 14.

Number of hydrogen bonds participate in intra molecular interaction of SARS-COV-2 main protease (red), docked complex of remdesivir (blue) and docked complex of carnosic acid (green) with respect to time.

Binding free energy analysis was employed to determine the binding capacity between ligands and receptor. It is noted that, binding free energy is the summation of all non-bonded interactions. Herein, we estimated the binding free energy (Das et al., 2020) of the main protease of SARS-COV-2 and docked complex of remdesivir as well as the docked complex of carnosic acid using the gmx mmpbsa tool (references).Within 100 NS, MD simulation, our binding energy analysis result, we found that the binding energy of remdesivir and carnosic acid are −98.858 and −47.326 kJ/mol, respectively. The binding energies have confirmed that the ligand remdesivir have shown a better affinity toward inhibition site of main protease. Thus, overall all the studied complexes are stable, but among all the complexes remdesivir with main protease are energetically more favorable and representing the stable complex. In docked complex of remdesivir we are getting energies for Van der Waals, electrostatic, polar salvation and SASA are −194.571 ± 13.915, −65.306 ± 15.059, 183.720 ± 26.385 and −22.701 ± 1.821 kJ/mol, respectively, where as in case of docked complex of carnosic acid the energies are −113.435 ± 10.446, −30.830 ± 8.724, 110.289 ± 12.763 and −13.351 ± 1.242 kJ/mol, respectively. From the overall results we found that the Van der Waals energy was higher (in terms of negative value) as compare to other energies in both the cases. The polar salvation energy was showing positive value in both the cases. The higher Van der Waals energy indicates more hydrophobic contact between remdesivir and carnosic acid with main protease. Electrostatic energy, Van der Waals energy, SASA energy, binding energy and polar solvation energy, with a standard deviation for the docked complex of remdesivir and carnosic acid with the main protease of SARS-COV-2 shown in Figure 15.

Figure 15.

Representations of Van der Waals energy, electrostatic energy, polar solvation energy, SASA energy and binding energy with their standard deviation values for the docked complex of (a) remdesivir and (b) carnosic acid with main protease of SARS-COV-2.

Analysis of collective movements of protein was performed through the PCA by using the gmx covar utility (Das et al., 2020; Havranek & Islam, 2020; Mahato & Fischer, 2018; Rane et al., 2020). Overall flexibility of studied trajectories for main protease, docked complex of remdesivir and docked complex of carnosic acid by trace of the covariance matrix was found 8.054, 7.683 and 8.3254 nm2, respectively. The lesser value of covariance matrix was found in case of docked complex of remdesivir as compare to other two systems. The eigen values plotted against the eigenvector index obtained from the of backbone atoms of SARS-COV-2 main protease, docked complex of remdesivir and docked complex of carnosic acid was shown in Figure 16.

Figure 16.

The eigen values plotted against the eigenvector index obtained from the of backbone atoms in water of SARS-COV-2 main protease (red), docked complex of remdesivir (blue) and docked complex of carnosic acid (green).

The conformational change during ligand binding was predicted by comparison the Gibbs free energy landscape analysis for PC1 and PC2 (Das et al., 2020; Rane et al., 2020). The free energy landscape was monitored by using gmx_sham utility of Gromacs module. The Gibbs free energy landscape (Figure 17) shows that the energy value varying from 0 to 14, 0 to 12 and 0 to 14 kJ/mol for main protease, docked complex of remdesivir and docked complex of carnosic acid, respectively. The regions of deep blue, light blue, green color signifies a lower, highly stable energy while red color signifies the higher energy conformation. In the free energy land scape plot, we obtained that, as there were two main free energy basins (Manoharan & Ghoshal, 2018) for both docked carnosic acid complex and docked remdesivir complex, only one main free energy basin was in the global free energy minimum region.

Figure 17.

Two dimensional figures of free energy landscape calculated for principal components PC1 and PC2 of native model of (a) main protease, (b) docked complex of remdesivir and (c) docked complex of carnosic acid, whereas three-dimensional figures of free energy landscape of (d) native model of the main protease, (e) docked complex of remdesivir and (f) docked complex of carnosic acid. The regions of deep blue, light blue, green color signifies a lower energy (highly stable) while red color signifies the higher energy conformation. The unit of energy is in kJ/mol.

Conclusions

In the present study, we have to assess the binding energies of darunavir, hydroxychloroquine, interferon alpha, lopinavir, niclosamide, ribavirin, ritonavir, umifenovir, ivermectin and remdesivir drugs and apigenin, drymaritin, carnosic acid, isoobtusitin, ellagic acid, morin, scutellarein, triterpenoids, naringenin, morusin, oxymatrine, honokiol, quercetin7-rhamnoside, myricetin and tryptanthrin phytochemicals with the main protease. Among photochemical, carnosic acid showed higher binding affinity, whereas within the studied drugs, remdesivir showed higher binding energy. The amino acid residues GLN110 and PHE294 were almost found in all the studied interactions of drugs and phytochemicals with the main protease of SARS-CoV-2. We studied Poisson-Boltzmann surface molecular mechanics (MM/PBSA) through molecular dynamics simulation. In the case of remdesivir, we found Van der Waals, electrostatic, polar solvation and SASA energies are −194.571 ± 13.915, −65.306 ± 15.059, 183.720 ± 26.385 and −22.701 ± 1.821 kJ/mol, respectively. The higher Van der Waals energy indicates more hydrophobic contact between remdesivir and carnosic acid with the main protease. The binding energies have confirmed that the ligand remdesivir have shown a better affinity toward inhibition site of the main protease. Potential use of those inhibitors as suitable candidates for SARS-CoV-2 drugs will be interesting research in experimental studies in future.

Click here for additional data file.