key: cord-1008253-j0qbz74g
authors: Garg, Saksham; Anand, Ashutosh; Lamba, Yash; Roy, Arpita
title: Molecular docking analysis of selected phytochemicals against SARS-CoV-2 M(pro) receptor
date: 2020-10-16
journal: Vegetos
DOI: 10.1007/s42535-020-00162-1
sha: a9b9a850d02985800db04174d2a027c63212d72b
doc_id: 1008253
cord_uid: j0qbz74g

Presently world is on a war with the novel coronavirus and with no immediate treatments available the scourge caused by the SARS-CoV-2 is increasing day by day. A lot of researches are going on for the potential drug candidate that could help the healthcare system in this fight. Plants are a natural data bank of bioactive compounds. Many phytochemicals are being studied for various ailments including cancer, bacterial and viral infections, etc. The present study aims to screen 38 bioactive compounds from 5 selected plants viz., Azadirachta indica, Curcuma longa, Zingiber officinale, Ocimum basilicum and Panax ginseng against SARS-CoV-2. Lipinski’s rule was taken as the foundation for initial screening. Shortlisted compounds were subjected to molecular docking study with M(pro) receptor present in SARS-CoV-2. The study identified that gedunin, epoxyazadiradione, nimbin and ginsenosides have potential to inhibit M(pro) activity and their binding energies are − 9.51 kcal/mol, − 8.47 kcal/mol, − 8.66 kcal/mol and − 9.63 kcal/mol respectively. Based on bioavailability radar studies gedunin and epoxyazadiradione are the two most potent compounds which are used for molecular dynamics simulation studies. Molecular dynamics studies showed that gedunin is more potent than epoxyazadiradione. To find the effectiveness and to propose the exact mechanism, in-vitro studies can be further performed on gedunin.

As per World Health Organization (WHO), on August 8, 2020, there were over 18.9 million confirmed cases of COVID-19 worldwide with the total death toll of over 709 thousand. Following the USA and Brazil, India has the thirdhighest confirmed COVID-19 cases with toll going over 1.7 million cases (https ://www.who.int/emerg encie s/disea ses/ novel -coron aviru s-2019/situa tion-repor ts/). In December 2019, an outbreak of pneumonia cases in Wuhan, Hubei Province, PRC (People's republic of China) was recorded. The unknown pathogen entity was soon identified as a novel coronavirus. Therefore, pneumonia caused by the virus was termed as Novel Coronavirus Infected Pneumonia (NCIP). Initially, the virus was called 2019-nCoV, which was later changed to SARS-CoV-2 by the International Committee on the Taxonomy of Viruses. The taxonomy was soon followed by the declaration of a pandemic on March 11, 2020 by the WHO .

Belonging to the family of Coronaviridae, coronaviruses are enveloped viruses having non-segmented RNA as their genetic material. The 2019-nCoV was isolated and identified from bronchoalveolar lavage fluid samples from a patient in Wuhan district on January 3, 2020. Full genome sequencing and phylogenetic analysis of the virus reported that 2019-nCoV is a distinct clade of beta-coronaviruses and related to SARS-CoV and MERS-CoV . It was also reported that 2019-nCoV shares 89.1% similarity with SARS-CoV . Hence, it was renamed to SARS-CoV-2. The structure of M pro from SARS-CoV-2 was given by Jin et al. (2020) . M pro (6LU7) is an inclusive and important part of the viral replication machinery and hence, making it a potential target receptor to combat against COVID-19 pandemic ).

The most common symptoms experienced by the patient on the onset of disease are cough, fever, and myalgia or fatigue while less common symptoms are headache, sputum production, haemoptysis, and diarrhoea . Organ dysfunction such as acute respiratory distress syndrome (ARDS), shock, acute cardiac injury, or even death can occur in severe cases . The SARS-CoV-2 virus has also been shown to disrupt the normal immune responses of the body which further prompts an uncontrolled inflammatory response in severe cases. Hence, the patients with severe cases manifest lymphopenia, lymphocyte activation and dysfunction, granulocyte and monocyte abnormalities, and an increase in immunoglobin G (IgG) .

The statistics shows that despite of the measures taken by different countries the daily number of cases are still rising. To control the spread there is an urgent need for a novel or repurposed drug to combat the pandemic. Bioinformatics significantly reduces the time and monetary costs that go into designing experiments, their execution, and laboratory trials. Molecular docking applies an innovative method by combining physiochemical principles with the complex scientific calculation algorithms. Docking assists in the drug development based on the characteristic interaction between the receptor and ligand/drug molecule. It helps industry to focus on the compounds with high potential against the particular target .

Polysaccharides derived from leaves of Azadirachta indica show anti-viral properties against poliovirus, antibovine herpes virus type 1, and duck plague virus (Kumar and Navaratnam 2013) . Rhizome extracted from Curcuma longa has inhibitory activity against influenza A neuraminidases (Dao et al. 2012 ) and H5N1 avian influenza virus (Sornpet et al. 2017) . Fresh Zingiber officinale was effective against the human respiratory syncytial virus (Chang et al. 2013) and chikungunya virus (Kaushik et al. 2020) . Ethanolic extracts and purified compounds from Ocimum basilicum exhibit a wide variety of anti-viral activities against hepatitis B antigen, HSV-1, and ADV-8 to name a few. Panax ginseng has been proven to inhibit the replication of human gammaherpesviruses (Kang et al. 2017) . Therefore, in this study, we have selected these five plants i.e. Azadirachta india (Indian lilac/Neem), Curcuma longa (turmeric), Zingiber officinale (Ginger), Ocimum basilicum (Basil) and Panax ginseng (Asian ginseng) to find a potential therapeutic natural compound against M pro receptor of SARS-CoV-2. The results obtained would help in developing an effective treatment for COVID-19 from plant source.

Main protease of SARS-CoV-2 was used in this study. The 3-dimensional structure was extracted in PDB format from the RCSB PDB data repository. PDB id given to the structure is 6LU7 and primarily structure is a homodimer having two A chains composed of 306 amino acids and N3 molecule acting as its inhibitor (Fig. 1) .

A total of 38 bioactive compounds from five different plants including Azadirachta indica, Curcuma longa, Zingiber officinale, Ocimum basilicum and Panax ginseng were selected as ligands and structures were obtained from PubChem databank in.sdf format (Table 1) . For the docking purpose, all the ligands were converted into.pdb file format using Biovia Discovery Studio Visualizer.

For the initial screening purposes, a web-based tool named SwissADME (https ://www.swiss adme.ch/) was used to eliminate a few compounds according to Lipinski's rule of five parameters. For a compound to qualify as ligand it should have < 500 Da molecular weight, a high lipophilicity i.e. value of Log P being less than 5, hydrogen bond acceptors being less than 10 and H-bond donors less than 5. Any compound with 2 or more violations was ruled out for further study (Lipinski 2004) . 

To obtain protein-ligand docked complex Autodock 4.2 was utilized. The downloaded structure of 6LU7 and each ligand was optimized prior to docking. From the protein 3D structure, water molecules and the inhibitor N3 molecule were removed. Addition of polar hydrogen bonds, Kollman charges and Gasteiger charges summed up the protein and ligand optimization. A grid box of 60 × 60 × 60 was prepared around the binding site of the protein with 0.375 Å spacing. Genetic algorithm was set as the search parameter and output was handled in Lamarckian GA run and docking log file (DLG) were obtained for further analysis of binding energy. The analysis of DLG file revealed a total of 10 conformations for each ligand. The conformation with highest negative binding energy was selected and docked complex was converted to a 2D structure to examine the interactions formed at binding site of 6LU7 with ligand.

Drug-likeliness of drug candidates with binding energy less than the control one was analyzed in a comprehensive way taking 6 physiochemical properties into consideration and forming a bioavailability radar using the SwissADME tool (https ://www.swiss adme.ch/). Six parameters considered: solubility, size, polarity, lipophilicity, flexibility and saturation. The pink shaded region defines the optimal values of the 6 parameters and deviation from these parameters on a large scale is suggestive of the ligand not being orally bioavailable (Diana et al. 2017). 

The receptor-ligand complex obtained after docking were subjected to molecular dynamic simulation using the Desmond-Maestro module 2020. The software provides high performance algorithms as its default settings and all the default setting were used to obtain high speed and precise results. The docked complex was first submerged in a TIP3P water model in an orthorhombic shape. The whole system was neutralized by adding 3 sodium ions at 0.15 M concentration. All the atoms were aligned using optimized potentials for liquid simulations-AA (OPLS-AA) 2005 force field. NVT was used as the ensemble class with SHAKE/RAT-TLE algorithm to limit the moment of the covelntly bonded atoms. The simulation was started at 300 K and 1 bar pressure for 50 ns. RESPA integrator was utilized to combine all the parameters of the dynamic simulation. A trajectory for 50 ns was set to show in 1000 frames to analyze the dynamic nature of the interaction and component stability.

Lipinski's rule of five was applied to estimate the drug likeliness of the all selected 38 candidates. This comparative method helps us to rule out few compounds according to their physiochemical properties. Compounds violating two or more parameters were out listed and rest of the compounds were considered to be ligands for the docking study. Out of 38 phytochemicals, 7 compounds: azadirachtin, oleoresin, epigallocatechin gallate, theaflavin digallate, hesperidin, neohesperidin and diosmin violated more than 1 parameter henceforth, remaining 32 compounds were subjected to docking studies (Table 1) .

All the filtered ligands from the ADME analysis were subjected to molecular docking analysis. Molecular docking is an essential computational tool in the drug discovery domain. It is done to further select the potential compounds and study the bond formation in the protein-ligand complex at the binding site. Figure 1a represents all 10 residues namely: THR24, THR26, ASN142, CYS145, PHE140, HIS163, HIS164, GLY143, GLU166, HIS172 which are present in at the active site of the M pro protein. N3 (native inhibitor) was taken as a control and comparative study of the docking results of all 31 ligands (Table 2 ) with the control revealed that four compounds having better binding energy as compared with the binding energy of N3 (-8.15 kcal/mol) (Fig. 6 ).

Among the 10 conformation of Nimbin, − 8.66 kcal/ mol was the least binding energy obtained. Five different types of interaction were observed including van der waals, H-bond, alkyl, pi-alkyl and carbon hydrogen bond (Fig. 2) . His163 and Cys145 forming the conventional H-bond while Pro168 and Met165 were engaged with a pi-alkyl and alkyl bond respectively. Glu166 and Leu167 were interacting with the ligand using carbon hydrogen bond and remaining residues weakly interact with the ligand via van der waals bond formation.

Gedunin-M pro complex (Fig. 3) had − 9.51 kcal/mol as the minimum binding energy. A sum total of 7 types of bond formation was observed. Glu166 formed convention H-bond, pi-anion bond and carbon hydrogen bond with gedunin. His41 forms a pi-sigma as well as pi-alkyl bonds while Cys145, Met49, and Leu27 formed alkyl bonds. His172 and Phe140 both formed carbon hydrogen bond. His143 contributed to stabilization by forming an additional pi-alkyl bond remaining all the residues were attracted to the ligand by van der waals bond.

Epoxyazadiradione was the other compound with binding energy greater than of control. It showed a binding energy of − 8.47 kcal/mol. Four types of interaction can be observed. Along with numerous residues involved in weak van der waals interaction, Cys145, Gly143 and Thr26 forms H-bond. Cys145 and Met165 formed alkyl bonds and Thr25 forms a carbon hydrogen bond (Fig. 4) .

Ginsenosides-M pro complex showed the minimum binding energy of -9.63 kcal/mol among all the conformations and ligands. A total of five different type of stabilizing interactions were observed. A conventional H-bond formation was done by Thr26 and Glu166 residues. His41 formed pi-sigma and pi-alkyl with 3 atoms of the ligand. Leu27, Cys145, Met165 and Met49 helped stabilizing complex via alkyl bond formation. Nine more residues can be observed around the ginsenosides interacting via van der waals forces (Fig. 5) .

M pro belongs to a protease class of proteins and this particular protease is seen to play essential role in replication process of many viruses (Chen 2020) . The main protease of SARS-CoV-2 was recently positioned in the pdb database as 6LU7 and provided the first potent target for the drug development process. The structural comparison revealed a high percentage (96%) of similarity between the SARS-CoV-2 and SARS virus main protease. The catalytic site is conversed in SARS-CoV-2 main protease and it is formed of two amino acids: His41 and Cys145 forming a catalytic dyad (Mirza and Froeyen 2020) . Catalytic dyad serves its purpose in maturation of the virus thus, making it an active target site in the protein (Chang 2010) . It can be observed that except nimbin all the other three ligands interact with the catalytic dyad residues in different fashion. Among the three, gedunin and ginsenosides had direct bond formation with both residues while epoxyazadiradione forming H-bond with only Cys145 and having weak Van der wall interaction with His41. Thus, all the four compounds can hold potential and might help in the process of finding the right cure to the ailment (Fig. 6 ).

Further analysis included four selected ligands viz., nimbin, gedunin, epoxyazadiradione and ginsenosides. A more illustrated and comprehensive study was done using bioavailability radar. Bioavailability radar is a descriptive tool to look into the drug-likeliness of the ligands based on six physiochemical properties. Study found gedunin and epoxyazadiradione to be orally bioavailable as the ligand radar entirely fit the pink shaded area. Nimbin was not an orally bioavailable ligand as it disobeyed the size parameter. Along with nimbin, ginsenosides having poor solubility in water and high value of lipophilicity was also predicted as not orally bioavailable candidate (Fig. 7) .

Two ligands (which completely satisfy bioavailability radar) i.e., gedunin and epoxyazadiradione were subjected to molecular dynamic simulation for 50 ns. − 117,519 kcal/ mol and − 117,476 kcal/mol which obtained potential energies for gedunin-6LU7 and epoxyazadiradione-6LU7 respectively. Total energy was another parameter 

The mean Root-mean-square deviation (RMSD) of gedunin-M pro complex and epoxtazadiradione-M pro complex was 2.017Å and 2.023Å respectively. Figure 8a indicates that RMSD for the gedunin-M pro complex was increasing for a brief period of time, however, the complex can be seen stabilizing toward the end of the simulation. The RMSD for epoxyazadiradione-M pro complex can be seen increasing throughout the simulation period. Small crust and troughs in RMSD can be seen scattered throughout the trajectory but they are due to the internal vibrations of the system. Root mean square fluctuation (RMSF) plot obtained assisted in understanding the stability of residues in the protein. Residues packed to form secondary structures like betasheet and alpha-helix confer less fluctuation thus have less RMSF while the loops and terminal residues have higher RMSF. Figure 8b reveals that both complex's residues follow the above stated phenomenon. However, most of the residues binding with epoxyazadiradione were observed to have higher RMSF value suggesting instability at the binding void.

The average radius of gyration (R g ) for gedunin-M pro complex and epoxyazadiradione-M pro complex were 4.099Å and 4.0813Å respectively. From Fig. 8c it can be seen that R g value for gedunin-M pro complex underwent some minor fluctuations, up to 0.1Å, which could be possible due to the packaging of the complex while epoxyazadiradione-M pro complex underwent some fluctuations, up to 0.15Å. Towards the end of the simulation, gedunin-M pro complex attained R g value near to its average value. Epoxyazadiradione-M pro complex also obtained a value near its average R g but it experiences more frequent fluctuations throughout the simulation.

The mean solvent-accessible surface area (SASA) for gedunin-M pro complex and epoxyazadiradione-M pro complex came out to be 171.412 ± 30.815 Å 2 and 269.354 ± 49.952 Å 2 respectively. Evident from Fig. 8d , the gedunin-M pro complex had attained a value about the mean value of SASA in the middle of the simulation. Epoxyazadiradione-M pro complex showed a very high peak in SASA towards the end of the simulation signifying interaction of hydrophobic and internal residues with the solvent which destabilizes M pro and is suggestive of greater chances of unfolding of the protein structure.

By analyzing the interactions of gedunin-M pro complex, all residues observed in docking can be seen interacting with the ligand. Few more additional residues were found interacting with gedunin line THR26, THR24, ASN142, etc. THR26 and GLU166 were interacting with the ligand 100% of the simulation time with two types of bonds namely, H-bond and water bridges thereby, holding the ligand in the binding pocket. Other residues like GLN189, SER144, HIS172 formed H-bond but for a very brief period of time.

Other types of interactions were also visible accounting for additional stability of the ligand in the binding void (Fig. 9a) . Study of interaction of epoxyazadiradione-M pro revealed that all the residues observed in docking results showed interaction with the ligand. THR26 and HIS41 were forming H-bond for interaction. THR24 formed around 40% of the time waster bridges. Other residues like MET49, GLN189, GLN306, GLU166 and many more added to the stability of the ligand in the void. Comparing the total number of contacts, the protein molecule was able to form number of bonds with gedunin as compared to epoxyazadiradione. During the 50 ns simulation time with gedunin only one instance was observed devoid of any contacts while with epoxyazadiradione lot of instances can be observed having zero contacts (Fig. 9b) . Computing and analyzing secondary structural analysis can be used to understand the protein packing and characteristics of folding of the protein with different ligands. Table 3 shows the percentage number of secondary structures of 6LU7 protein with gedunin epoxyazadiradione ligands. With gedunin 40.23% of total residues formed secondary structure while binding of epoxyazadiradione with the protein led to dropping of number to 40.09%.

In the present study 38 compounds were selected from five plants. These compounds were screened using Lipinski's rule of five and determined drug-likelihood of the compound. 31 compounds were drug-likeable which were subjected to molecular docking. Docking results were compared with the native inhibitor of M pro i.e., N3 and four ligands were found more potent in terms of binding energy. All four ligands were studies using bioavailability radar. Our results proposed gedunin and epoxyazadiradione to be potent inhibitors of the 6LU7 along with nimbin but nimbin was not found orally bioavailable. Similarly, ginsenosides showed the best docking result but not the bioavailability radar result. Further molecular dynamics showed that gedunin is more potent than epoxyazadiradione. To find the effectiveness and to propose the exact mechanism in-vitro studies can be encouraged on gedunin, epoxyazadiradione and nimbin to find a potent cure for the COVID-19.

Quaternary structure of the SARS coronavirus main protease

Fresh ginger (Zingiber officinale) has anti-viral activity against human respiratory syncytial virus in human respiratory tract cell lines

Pathogenicity and transmissibility of 2019-nCoV-A quick overview and comparison with other emerging viruses

SwissADME: a free web tool to evaluate pharmacokinetics, drug-likeness and medicinal chemistry friendliness of small molecules

Curcuminoids from Curcuma longa and their inhibitory activities on influenza A neuraminidases

In silico screening of potential Chinese herbal medicine against COVID-19 by targeting SARS-CoV-2 3CLpro and angiotensin converting enzyme II using molecular docking

Structure of Mpro from SARS-CoV-2 and discovery of its inhibitors

Antiviral activity of 20(R)-ginsenoside Rh2 against murine gammaherpesvirus

Anti-viral activity of Zingiber officinale (Ginger) ingredients against the Chikungunya virus

Neem (Azadirachta indica): Prehistory to contemporary medicinal uses to humankind

Lead-and drug-like compounds: the rule-of-five revolution

Structural elucidation of SARS-CoV-2 vital proteins: Computational methods reveal potential drug candidates against main protease, Nsp12 polymerase and Nsp13 helicase

Antiviral activity of five Asian medicinal pant crude extracts against highly pathogenic H5N1 avian influenza virus

Clinical characteristics of 138 hospitalized patients with 2019 novel Coronavirus-infected pneumonia in Wuhan

A new coronavirus associated with human respiratory disease in China

COVID-19: immunopathogenesis and Immunotherapeutics