In silico investigation of phytoconstituents from Indian medicinal herb 'Tinospora cordifolia (giloy)' against SARS-CoV-2 (COVID-19) by molecular dynamics approach.

The recent appearance of COVID-19 virus has created a global crisis due to unavailability of any vaccine or drug that can effectively and deterministically work against it. Naturally, different possibilities (including herbal medicines having known therapeutic significance) have been explored by the scientists. The systematic scientific study (beginning with in silico study) of herbal medicines in particular and any drug in general is now possible as the structural components (proteins) of COVID-19 are already characterized. The main protease of COVID-19 virus is Mpro or 3CLpro which is a key CoV enzyme and an attractive drug target as it plays a pivotal role in mediating viral replication and transcription. In the present study, 3CLpro is used to study drug:3CLpro interactions and thus to investigate whether all or any of the main chemical constituents of Tinospora cordifolia (e.g. berberine (C20H18NO4), β-sitosterol (C29H50O), coline (C5H14NO), tetrahydropalmatine (C21H25NO4) and octacosanol (C28H58O)) can be used as an anti-viral drug against SARS-CoV-2. The in silico study performed using tools of network pharmacology, molecular docking including molecular dynamics have revealed that among all considered phytochemicals in Tinospora cordifolia, berberine can regulate 3CLpro protein's function due to its easy inhibition and thus can control viral replication. The selection of Tinospora cordifolia was motivated by the fact that the main constituents of it are known to be responsible for various antiviral activities and the treatment of jaundice, rheumatism, diabetes, etc. Communicated by Ramaswamy H. Sarma.

Introduction

In 2020, with the outbreak of a newly detected corona virus, the whole world is witnessiene">ng a paene">ndemic situatioene">n which has put n class="Species">human life in crisis as more than 15000000 numbers of active patients with 619150 people have already died due to this virus in last couple of months (https://covid19.who.int/). The viruses, which cause illness in animals and human, coronavirus is the member of that family of viruses. In human, common cold, Middle East Respiratory Syndrome (MERS) and Severe Acute Respiratory Syndrome (SARS) are the different types of respiratory infections which these viruses can cause (Beadling & Slifka, 2004; de Wit et al., 2016). COVID-19 is the most recently declared disease by WHO caused by a new member “novel coronavirus” of this virus family (https://www.who.int/, Guo et al., 2020). Before its outbreak appeared in Wuhan, China, in late December 2019 (Guan et al., 2020; Wang, Hu, et al., 2020), this new virus and disease were unknown to the human civilization. On 11 February 2020, the international committee on taxonomy of viruses declared the new novel coronavirus as “Severe Acute Respiratory Syndrome Coronavirus 2” (SARS-CoV-2) (Gorbalenya et al., 2020) (originally named 2019-nCoV). This virus is part of the Coronaviridae family of virus SARS-CoV having a genome sequences of 79.5% sequence matching (Graham et al., 2013; van der Hoek et al., 2004). The surface glycoproteins of this virus creates a crown-like appearance which can be visible under electron microscope. Crown like spike proteins are the important part of SARS-CoV-2 (Hendaus & Jomha, 2020). It recognizes some specific human proteins. The human proteins which coats inside of the nose and the cells of lungs are mainly can interact with the spike protein of the SARS-CoV-2 virus. After interaction the two proteins bind together. In the product form, the spike protein of CoV-2 changes its shape and causes the human receptor cell to engulf the virus. After entering human body, it then replicates itself and can infect neighbouring cells and tissues by damaging mainly lung, heart, brain cells and many other organs. After few weeks of SARS-CoV-2’s appearance in human society, its structural components were characterized. Results are also available on how the virus targets, interacts with human proteins and how we can fight against it (Jin et al., 2020), but till date we are searching for the answer to fight against this virus. Till now structural biologists have identified more than 160 structures of nine various kinds of SERS CoV-2 proteins by X-ray crystallography. In human, the principal target of vaccine or drug to act on is the spike protein that protrudes from the lipid shell of the virus. The infectivity of the virus would be much reduced and maybe even eliminated if and only if the human immune system are primed to recognize and counteract the spike protein of SERS CoV 2 (Khan, Zia, et al., 2020). COVID-19 has now become a pandemic worldwide affecting mostly all countries globally. It is an infectious disease with a high potential for transmission for human to human close contacts through the respiratory droplets (such as coughing) and by fomites that can propagate through air at a minimum distance of 1 m (Rothe et al., 2020). Research works performed till now suggests that maintaining a distance of more than 1 m between two individuals which is termed as ‘social distancing’ or ‘physical distancing’ along with proper hand-hygiene reduces the chance of being infected by COVID-19. Though there are many predictions about the airborne transmission of this disease, but until now no scientifically valid evidence is available (Doremalen et al., 2020). Fever, dry cough, tiredness, headache, diarrhoea, skin infection and less working of sensing organs like test less feeling of food, less sensing of smell etc are the most common symptoms of COVID-19. These symptoms are appeared usually in mild form and begin gradually afterwards. Few infected people even become asymptomatic or only have very mild symptoms. The main point of concern is that till now we don’t have any specific therapies for COVID-19. Research results regarding the treatment of COVID-19 are not very productive also (Weiss & Murdoch, 2020). As per current data of July 24th, 2020, over 15012731 cases of COVID-19 active patients have been confirmed worldwide and over 619150 death cases have been reported so far. Due to nonavailability of any vaccine or medicine, the supportive care and non-specific treatment to ameliorate the symptoms of the patient are the only options currently to treat COVID-19 (Salata et al., 2019). However, extensive research works leading to clinical trials of both western and traditional medicines are going on all over the world to develop vaccines and medicines to prevent and cure CoV-2 infections (Bhardwaj et al., 2020; Enmozhi et al., 2020; Jin et al., 2020; Rameez et al., 2020). To combat COVID-19, some preliminary research observations have indicated that the combinations of some clinically applied anti-malarial drugs (e.g. Chloroquine, Hydroxychloroquine) (Chang et al., 2016; Khan, Zia, et al., 2020) and anti-HIV vaccines can be used to treat COVID-19. Also some conventional FDA approved drugs (Remdesivir, Nelfinavir, Paritaprevir, Raltegravir, Praziquantel, etc.) are being tried as the potential drug against this CoV-2 and found with certain curative effect in vitro (Hendaus, 2020; Khan, Jha, et al., 2020; Wang, Cao, et al., 2020). But till no none of the clinically applied drug/vaccine response is found to be very encouraging. Also toxicity of the currently applied drugs remain an inevitable issue causing serious adverse effects to the patients (Grein et al., 2020; Kaisari & Borruat, 2020; Wang, Guo, et al., 2020). In short, no specific therapy, medicine or vacine for the effective treatment of COVID-19 has yet been reported. This unavailability of cure to this disease motivated us to investigate the possibility of inhibition of CoV-2 by some phytoconstituents available from some Indian medicinal plants. It is a well-known fact that in the treatment of pandemic and endemic diseases, Medicinal herb extracts have accumulated thousands of years’ experience (Jee et al., 2018; Kumar et al., 2018). To get a valid answer to act against COVID-19, various complementary and alternative vaccine, medicine and treatment methods are also needed for the managements of patients with SARS-CoV-2 infection. So research works with different phytoconstituents from various traditional medicinal herbs is certainly worth of investigation (Khaerunnisa et al., 2020) in this present emergency situation. To deal with this emergency for the present work, we have chosen an Indian medicinal plant: Guduchi Pippali (Giloy) or Tinospora cordifolia (Figure 1). Tinospora cordifolia is the member of the family of Menispermaceae. It is usually found in mostly Asian counties like India, Myanmar, Sri Lanka, and China. The plant is commonly used as a main component in many traditional Ayurvedic medicines. It is used as a parent medicine for the therapies against several common diseases like jaundice, rheumatism, urinary disorder, skin diseases, diabetes, anemia, inflammation, allergic condition, etc. (Sonkamble & Kamble, 2015). For the treatment of helminthiasis, heart diseases, leprosy, rheumatoid arthritis, etc the stem of the plant is very useful. It also supports the immune system by increasing the body's resistance to various infections, supports standard white blood cell structure, function, and levels (Sharma et al., 2019). All of the above mentioned pharmacological actions of Tinospora cordifolia originates from its chemical constituents of different parts like leaves, stem, root, flower, seed etc. The various categories of chemical constituents of different classes such as alkaloids, glycosides, steroids, phenolics, aliphatic compounds, polysaccharides essential oils, a mixture of fatty acids, and polysaccharides are present in different part of the plant body, including root and stem. Phytochemistry of all these are well documented in the literature (Sharma et al., 2019).

Figure 1.

Morphology of Tinospora cordifolia.

Estimation of pharmacological properties of small molecules when it interacts with target protein is as a crucial step for the drug discovery. We have selected a few small sized phytocoene">nstituents extracted compouene">nds from n class="Species">Tinospora cordifolia (Figure 2) (https://www.drugbank.ca/drugs). These phytoconstituents were selected after applying some proper virtual screening by evaluating their drug-likeness, pharmacokinetics and lipophilicity properties which are set of guidelines for identification of potential drug compounds. After virtual screening, we have evaluated their potential inhibition properties against CoV-2 main protease by the mechanism of molecular docking. After studying proper binding affinity of proposed drugs against the target receptor by the Molecular docking, finally we have verified the drug inhibition with target receptor by different interaction energies like hydrogen bonding, hydrophobic, coulombic etc by MD simulation approach to validate the applicability of proposed drug molecule as suggested by Anuj et al. (2020). In vivo estimation of drug molecules is time consuming and expensive so to deal with the present COVID crisis, in silico methods have become inevitable as we need an urgent solution of SARS-CoV-2 infection.

Figure 2.

Different phytoconstituents extracted compounds from Guduchi Pippali (Giloy) or Tinospora cordifolia as ligand molecules. [# from Gausswiew, rest from Drugbank].

Different phytoconstituents extracted compounds from Guduchi Pippali (Giloy) or Tinospora cordifolia as ligaene">nd molecules. [# from Gausswiew, rest from Drugbaene">nk].

Important pharmacological information of various chemical constituents have indicated the research community to explore computational molecular docking and MD methods for investigating the drug-protein interactions. Variety of simulation methods and databases have yet been used for the in-silico prediction of target drugs (Cheng et al., 2013; Boopathi et al., 2020). Many groups are presently working on different quantum mechanical simulation techniques to search varieties of interacted complex configuratioene">ns of differeene">nt orgaene">nic aene">nd iene">norgaene">nic systems aene">nd their iene">nteractioene">ns with differeene">nt eene">nviroene">nmeene">nts (Raene">na & Chowdhury, 2017). We have earlier studied the eene">nviroene">nmeene">ntal effects oene">n differeene">nt orgaene">nic probe molecules iene">n view of their applicatioene">ns iene">n various biomedical field (Siene">ngla & Chowdhury, 2014). In the preseene">nted study, we have reported some eene">ncouragiene">ng respoene">nses we obtaiene">ned iene">n terms of iene">nhibitioene">n poteene">ntials from some of the tested phytocoene">nstitueene">nts of curreene">nt probe system: n class="Species">Tinospora cordifolia. We strongly consider that the outputs of the present work will lead to some important insights into the development of alternative drugs for COVID-19.

Materials and methods

Potential target protein structures for SARS-CoV-2

Coronaviruses have positive-sense siene">ngle-straene">nded RNA. n class="Species">Coronavirus SARS-CoV-2, responsible for COVID −19 is member of β-coronavirus genus (Fehr & Perlman, 2015). CoV-2 contains genome size of ∼30 kilobases. It encodes for various non-structural and structural proteins. The structural protein components of SARS-CoV-2 comprise: envelope (E), membrane (M), spike (S) and nucleocapsid (N) (Gorden et al., 2020; Woo et al., 2005). Structure of SARS-CoV-2 has been identified quickly after it first appearance and its genomic sequence is already available (Wu et al., 2020) for us. As the virus has been discovered very recently so very few available immunological information (immunogenic epitopes eliciting antibody or T cell responses) about the virus are available until now. It is a fact that until now, no available therapeutics or approved and tested vaccines/medicine are there to act against any human-infecting coronaviruses. The newly invented SARS-CoV-2 is very closely related to SARS-CoV structures. So utilization of this known protein structure can be applicable for quick discovery of drugs against this newly appeared virus (Jiang et al., 2020). It is reported very recently that CoV-2 enters into host cells by the spike (S) glycoprotein and forms homo trimers protruding from the viral surface (Jin et al., 2020). This spike protein interacts strongly with the human ACE2 (angiotensin-converting enzyme 2) receptor (Velavan & Meyer, 2020). After entering the host cell, SARS-CoV-2 replicates itself through some cyclic processes. First it translates its genomic RNA (gRNA). Proteolysis of the translated polyprotein takes place with viral 3C-like proteinase. After that replication of gRNA takes place. The viral replication complex appears which consists of RNA dependent RNA polymerase (RdRp), helicase, 30-to-50 exonuclease, endoRNAse, and 20-O-ribose methyltransferase. Lastly the assembly of viral components takes place. The proteins (S) which are associated with replication, are the primary targets of post-entry treatment drugs or vaccine design which can suppress viral replication/infection or neutralizing antibodies (Abs) upon infection (Gao et al., 2020; Gupta et al., 2020; Khan, Zia, et al., 2020; Sarma et al., 2020). SARS-CoV-2 3C like proteinase is already predicted to bind with different FDA approved antiviral commercially available drugs like atazanavir, lopinavir, ritonavir, remdesivir, efavirenz and some other antiviral drugs which have a predicted affinity with clear efficacy value of K > 100 nM potency (Hendaus, 2020; Khan, Jha, et al., 2020). Prediction suggests that viral proteinase targeting drugs should act very efficiently on the viral replication process. In support of the prediction, it was reported by the docking study of some HIV proteinase inhibitors of the CoV proteinase that lopinavir, atazanavir and ritonavir may inhibit the CoV proteinase (Dayer et al., 2017). Like lopinavir, atazanavir and ritonavir there are some other case studies also with other inhibitor drugs like hydroxyxhloroquine, remidisivir etc., but until now there is no real evidence about whether these drugs can act efficiently as predicted against COVID-19. Traditional drug discovery usually takes years of research and trial. Now the whole world is dealing with an emergency situation with an urgent need of required drug against this virus to save the infected lives. An alternative approach to complement the existing procedure and to combat this problem can be to use the computer aided drug design (CADD) or computer-assisted structure-based drug design (SBDD). For the present work, some particular type of SBDD approaches have been implemented. Specifically, we tried here the in silico docking model and MD simulation approaches to search for medication of COVID-19 by using two most variable viral proteins (I, II) SARS-CoV-2 protease enzyme Mpro or 3CLpro as the receptor (Figure 3(a, b)).

Figure 3.

(a). Target variable viral proteins (I, II) SARS-CoV-2 protease enzyme Mpro or 3CLpro as the receptor. (b) Ramachandran plot for both the receptor protein.

(a). Target variable viral proteins (I, II) SARS-CoV-2 protease enzyme n class="Gene">Mpro or 3CLpro as the receptor. (b) Ramachandran plot for both the receptor protein.

3CLpro is a key to n class="Chemical">SERS-CoV enzyme. CoV 3CLpro is responsible for the maturation of itself and the subsequent maturation of the replicate polyproteins. This protein has a conserved catalytic domain. The protein can control functions of the virus like: replication and transcription processes. Such important functionalities of 3CLpro makes it an ideal target for drug development (Anand et al., 2003). Also due to its viral life cycle and the absence of closely related homologues in human, the 3CLpro protease has become most attractive receptor protein for antiviral drug design. Recently 3D structure of one of the 3CLpro like protease protein is reported by X-Ray crystallographic data (PDBID: 6LU7) (Jin et al., 2020). Protease proteins are essential to the transmission and virulence of the virus. For an effective therapy, by inhibiting these protease proteins with proper drug inhibitor, the severity of the infection can definitely be reduced. We have tried to check the inhibiting and binding possibilities of two of these natural 3CLpro like protease protein substrates in order to find the possibilities of effective drug design for COVID-19 in the present work.

Protein receptors preparation

In this work we have used two CoV-2 n class="Chemical">3CLpro main proteases as target of the potential drug molecules. The proteases used are: protease in complex with an inhibitor N3 (PDBID: 6LU7) (Jin et al., 2020; Zhang et al., 2020) referred to as I and the other one is protease bound to potent broad-spectrum non-covalent inhibitor X77 (PDBID: 6W63) referred to as II (Mesecar, 2020). The 3D structures of both were downloaded from the Protein Data Bank (https://www.rcsb.org/) and are shown in (Figure 3(a, b)) without added water and hydrogen. All the available required properties of these proteins are described in Table SD1. At the beginning of the present investigation, all protein structures have been cleaned by removing the existing lead components. Also the existing water molecules and ions have been removed from protein structures. After cleaning, Gasteiger charges of protein structures have been calculated and procedure of inclusion of polar hydrogens have been performed. Then subsequently the merging of non-polar and rotatable bonds were defined. All of the above mentioned steps have been performed by Auto Dock and MG Tools of Auto DockVina software (Trott & Olson, 2010). Lastly the inbuilt ligands were removed from the protein molecules using Discovery studio 2020 (Dassault Systèmes BIOVIA, 2017) and the cleaned and final proteins structures were saved in PDB format. Active amino acid residues which will be involved in mediating enzyme activity for protein-ligand interaction were predicted for each target protein using Ramachandran plot (Figure 3(a, b)) by visualizing the dihedral angles ψ against φ of amino acid residues. This two-dimensional plot shows the allowed and disfavoured values of ψ and φ. In the present work for our target 3CLpro proteins, each plot specifies localization on the A chain residues. This localization determines the quality of the protein structure. Good quality reflects the efficient and accurate docking possibilities. In the present work, the Ramachandran plot for both the target proteins show good localization of residues which are the indication of good probable docking results.

Ligand drug molecules preparations

For drug discovery, estimatiene">ng required iene">n vivo aene">nd iene">n silico pharmacological properties of proposed small sized molecules are reflected as a very essential step. In vivo estimatioene">n is expensive aene">nd time coene">nsumiene">ng aene">nd so iene">n present scenario iene">n silico methods have become certaiene">n. For iene">n silico estimatioene">n the process is knowene">n as virtual screeniene">ng. In virtual screeniene">ng different n class="Disease">Drug-likeness rules (a set of guidelines for the determination of structural properties of proposed drug compounds) are used for preliminary screening of required drug-like properties of a molecule. Some mostly industry applicable drug-likeness rules are: Lipinski's rule, MDDR-like rule, Veber’s rule, Ghose filter, Egan rule, Muegge rule etc. Out of these different methods (rules), one of the dominant method is Lipinski's rule of five (Ro5) which is known as “a rule of thumb”. It is used to evaluate if a proposed chemical compound possesses certain required pharmacological properties which would make the compound a potential candidate for orally active drug in humans (Lipinski, 2004). Candidate drugs that imitate Ro5, usually have an increased chance of reaching the market since they tend to have lower attrition rates during clinical trials. Ro5, with its filters (Molecular weight ≤500, H-bond donors ≤ 5, H-bond acceptor ≤ 10, MLOGP ≤ 4.15, molar refractivity should be between 30 and 140) helps in distinguishing between non-drug like and drug like molecules which helps to circumvent preclinical and clinical failures by early preclinical development of drugs. We have done the virtual screening of our proposed drug molecules by applying most of the above mentioned drug-likeness rules (Veber et al., 2002) including some other screening properties like polar surface area (Daina et al., 2017), physicochemical properties, pharmacokinetics, lipophilicity, synthetic accessibility etc (Table 1). ADME (Adsorption, Distribution, Metabolism and Excreation) is important technique to analyze the pharmacodynamics of the proposed drug molecule. SWISS-ADME software (https://www.swissadme.ch) helps us to identify proposed drug molecule by applying different virtual screening methods. Different components of lipophilicity (iLOGP, WLOGP, XLOGP3, MLOGP, Log Po/w), water solubility (Log S (SILICOS-IT)), pharmacokinetics (GI absorption, BBB permeant, P-gp substrate, Log (Kp)) also help us to preliminary test the suitable drug molecule. For choosing the proposed drug molecules through virtual screening, sometimes we found some of them have violated some selection rules like Ro5, Veber etc. In drug industry today there are several important drugs available in market who violate some likeness rules. Some very popular drugs like dabigitran etexilate, bromocriptine mesylate, olmesartan medoxomil, fosinapril, and reserpine etc revealed two Ro5 rule violations (Giacomini et al., 2010). Some very popular tyrosine kinase inhibitors (e.g. lapatinib and nilotinib) and HIV protease inhibitors (e.g. lopinavir, atazanavir, nelfinavir) do not follow the Ro5 criteria. Atazanavir exhibits three Ro5 violations. Keeping all these points in mind we have considered some natural compounds as our potential inhibitors for 3CLpro protease although some of them violate one or two criteria of Ro5 but validate other criteria. After thoroughly testing the proposed phytoconstituents by all mentioned virtual screening of drug-likeness properties we have tested the drug binding possibilities with target 3CLpro main proteases by molecular docking and further MD simulation mechanisms. Specifically, based on the positive responses from virtual screening, we have selected five phytoconstituents from Tinospora cordifolia as potential ligand drugs: berberine (C20H18NO4), choline (C5H14NO) and tetrahydropalmatine (C21H25NO4) from the group of alkaloids, β-sitosterol (C29H50O) from steroids, octacosanol (C28H58O) from aliphatic group. Details of the structure of these molecules were downloaded from Drug Bank in pdb format and are described in Figure 2 and Table 1(a–e) along with their various chemical, physical, pharmakokinetics and drug likeness properties obtained from SWISS ADME. All 2D and 3D view of the proposed ligands were obtained from Drugbank website. For octacosanol we have used Gaussian supported Gaussview to draw 3D structure (Figure 2). All of our proposed drug molecules have molecular weight less than 500 g/mol. All of them have TPSA values less than 50 Å2. All drug molecules have H-bond donors ≤ 5, H-bond acceptor ≤ 10. All of the proposed drug molecules have synthetic accessibility count between less than 10 so that they can be synthesized easily. Though β-sitosterol and octacosanol have violated some of the drug likeness properties, still we have consider these phytoconstituents as potential inhibitors due to their availability in drug industry. For example, octacosanol is used as a probable therapeutic agent for treating Parkinson's disease. It also possess antiaggregatory properties and cholesterol-lowering effects. It is used as a component of Ginsamin Power as a tablet from Biogrand Co., Ltd (https://www.drugbank.ca/). It has passed the clinical trial also for the treatment of high Cholesterol condition (Taylor et al., 2003). For Molecular docking study, the ligand file is required in pbdqt format. These ligand drug molecules have been saved in pdbqt format by Auto Dock Tools 1.5.6 (Michel, 1999).

Table 1.

Molecular configuration and drug likeliness properties of proposed ligand drug molecules for COVID-19 by SWISS ADME data.

Drug likeness properties	berberine (a)		choline (b)	β-sitosterol (c)	tetrahydropalmatine (d)	octacosanol (e)
	Physicochemical Properties
Molecular weight (gm/mol)	336.361		104.1708	414.718	355.434	410.7595
Num.H-bond acceptors	4		1	1	5	1
Num. H-bond donors	0		1	1	0	1
No of rotatable bonds	2		2	6	4	26
Molar Refractivity	94.87		30	133.23	103.99	137.87
Topological polar surface area TPSA (Å2)	40.80 Ų		2 0.23 Ų	20.23 Ų	40.16 Ų	20.23 Ų
	Lipophilicity
Log Po/w (iLOGP)	−0.00		−2.41	5.05	3.69	7.20
Log Po/w (XLOGP3)	3.62		−0.40	9.34	3.24	13.61
Log Po/w (WLOGP)	2.19		−0.32	8.02	2.52	10.14
Log Po/w (MLOGP)	2.53		−3.46	6.73	2.20	7.07
Log Po/w (SILICOS-IT)	7.24		−0.57	7.04	3.75	10.96
Concensus Log Po/w	9.80		−1.38	7.24	3.08	9.80
	Druglikeness
Lipinski	yes		yes	yes	yes	yes(1 violation)
Veber	yes		yes	yes	yes	yes(1 violation)
Ghose	yes		Partly yes (2 violation)	No	yes	No
Egan	yes		yes	yes (1 violation)	yes	yes(1 violation)
Muegge	yes		yes (1 violation)	Partly yes (2 violation)	yes	No
Bioavailability score	0.55		0.55	0.55	0.55	0.55
Synthetic accessibility (SA)	3.14		1	6.30	3.59	3.72
	Pharmacokinetics
GI absorption	High		Yes	Low	High	Low
BBB permeant	Yes		No	No	Yes	No
P-gp substrate	Yes		No	No	Yes	Yes
Log Kp (skin permeation)	−5.78 cm/s		−7.22 cm/s	−2.20 cm/s	−6.17 cm/s	0.86 cm/s
	Water Solubility
Log S (SILICOS-IT)	−5.92		−1.26	−6.19	−5.87	−10.53
Solubility (mg/mL)	4.00e − 04		5.74e + 00	2.69e − 04	4.81e − 04	1.20e − 08
Toxicity	Rat LD50: 2.7834 mol/kg		Oral rat LD50: 3400 mg/kg	Not Available	Not Available	Not Available

Molecular configuration aene">nd drug likeliene">ness properties of proposed ligaene">nd drug molecules for n class="Disease">COVID-19 by SWISS ADME data.

Molecular docking and visualization

Molecular docking predicts the potential drug-target interactions by energy minimization and binding energy calculations. Interaction between small molecules (ligand) and protein receptor (may be an enzyme) creates the possibility of inhibition of ligand towards receptor enzyme. So output of molecular docking can demonstrate the feasibility of aene">ny biochemical reactioene">n. It caene">n also predict aene">n optimized orientatioene">n of ligaene">nd towards the target receptor. By different biene">ndiene">ng modes of ligaene">nd with target receptor, dockiene">ng result may provide a raw material for the n class="Species">rational drug designing. In combination with scoring function, Molecular docking method is used for finding more potent, selective and efficient drug candidates. The docking algorithm performs through some cyclical processes. The cyclic process by specific scoring functions helps to identify the suitable ligand conformation. After reaching to a minimum energy, the process converges to a solution (Yuriev et al., 2011). Docking mechanism uses inhibitor drug’s binding properties with nucleic acid of target protein. The binding property establishes a correlation between drug’s molecular structure and cytotoxicity. Molecular docking can also predict whether the ligand/drug is docked with the receptor protein/DNA or not. Docking algorithm gives output in the form of quantitative predictions of binding energetics by providing the rankings of docked compounds. The ranking is based on the binding affinity of ligand: receptor complexes. The best confirmation of the ligand: receptor complex is that or those which has/have lower binding energy. After getting positive prediction, further experimental procedures are applied to the development of new drug molecule. Docking-based studies on the suggested inhibitors onto the protease of CoV-2 through Receptor: ligand docking analysis was performed using AutoDock Vina (Trott & Olson, 2010). The macromolecule (protein) file is saved in pdbqt format by Auto Dock Tools (Trott & Olson, 2010) and ready to be used for docking. Similarly the ligand molecule has been saved in pdbqt format by Auto Dock Tools. For AutoDock Vina algorithm, for setting up the configuration file following parameters were used: (i) count of binding modes- 9; (ii) exhaustiveness − 8 and (iii) applied maximum energy difference- 3 kcal/mol. Ligand centered maps were generated and Gridbox center was set to coordinate x y, and z of residue position of the target protein respectively.

Out of all the possible poses (optimized ligand: protein complex structure) suggested by simulation according to the binding modes, the pose showing maximum hydrogen boene">nds aene">nd miene">nimum biene">ndiene">ng affiene">nity(kcal/mol), were chosen as the best ligaene">nd: proteiene">n complex structure formed by ligaene">nd-receptor iene">nteractioene">n. For choosiene">ng the best possible ligaene">nd: proteiene">n complex structure we have taken the help of root meaene">ns square deviatioene">n (RMSD) method. For RMSD evaluatioene">n oene">nly movable heavy atoms (ligaene">nd) were coene">nsidered to identify the best mode. Two variaene">nts of RMSD: RMSD/lb (lower bound) aene">nd RMSD/ub (upper bound) were used here for the whole simulatioene">n. The two variaene">nts differ iene">n how the atoms are matched iene">n the distaene">nce calculatioene">ns. For the distaene">nce for upper bound is coene">nsidered when the ligaene">nd has no symmetry as RMSD/ub matches each atom iene">n a coene">nformatioene">n with itself iene">n the other coene">nformatioene">n igene">noriene">ng aene">ny symmetry. We have used RMSD/ub as the ligaene">nd molecules have no specific symmetry. Aene">nother criteria of choosiene">ng best ligaene">nd: proteiene">n pose is identifyiene">ng the types aene">nd number of boene">ndiene">ng between them. The metabolite which makes maximum number of H-boene">nds, hydrophobic boene">nds with the receptor proteiene">n mostly show better capability to form ligaene">nd: proteiene">n complex formatioene">n by iene">ntermolecular iene">nteractioene">n between iene">nhibitor aene">nd receptor. Different output poses were aene">nalysed iene">n Disn class="Species">covery Studio visualizer 2020 version 20.1.0.19295 (Dassault Systèmes BIOVIA, 2017) for the formation of non-bonded hydrogen bonds. The best pose structure was analysed also by their binding affinity, inhibition constants and other supporting interactions.

MD Simulations

At atomistic level, MD simulation results helped us to investigate the structural dynamics of receptor CoV-2 protease (n class="Chemical">3CLpro) upon binding with small ligand (proposed drug) molecule. LINUX based platform ‘‘GROMACS 5.1 Package’’ (Berendsen et al., 1995) was used for determination of thermodynamics stability of the proposed ligand: protein complex. MD simulations were performed using latest CHARMM36 all atom force field (Soteras et al., 2016). Before simulation, separate topologies were prepared for receptor and ligand by different external tools (Gajula et al., 2016). For our ligand, before creating the topology first we have optimized the ligand structure by Gaussian 9.0 by Density Functional Theory (DFT) with the basis set 6.31G (d,p) (Becke, 1997; Frisch, 2004). With the optimized structure we have generated the ligand topology with the CHARMM General Force Field (CGenFF) program (Vanommeslaeghe & MacKerell, 2012). In CHARMM all-atom force field, all H atoms of ligand are explicitly represented. For topology we have used the cgenff_charmm2gmx.py (http://mackerell.umaryland.edu/charmm_ff.shtml#gromacs). To perform simulations in aqueous solution, we have used a well known water model: TIP3P. During solvation process the bare protein and protein: ligand complex were solvated in the cubic box having specific boundary conditions (with 10 Å buffer distance) with volume as 893,000 Å3. As per procedure 4Na+ ions were added for electrically neutralizing the probe system. Before simulation we have done the energy minimization on the system to sort out any bad starting structures and also to minimize solute structure in vacuum before introducing solvent molecules. The Protein ligand complex system was equilibrated under suitable simulation parameters consistent for our selected CHARMM General Force Field for building the energy minimized solvated system. The steepest descent algorithm has been used for energy minimization of the system with varying time (ps) for 500,000 iteration steps. For reducing the steric clashes, the applied algorithm has a cut-off up to 1000 kJmol−1 (Anuj et al., 2020). The MD simulation has started with the minimization of the system. The minimization has been achieved at two phases each having 500,000 steps. In first part, equilibration was obtained having each step of 2 fs with a boundary condition of constant number of particles (N), volume (V), and temperature (T). We have used 10000-ps NVT equilibration. In second phase, the equilibration was achieved under the pressure of 1 atmosphere 298 K. Here the boundary condition maintained as constant NPT (particle numbers, pressure, temperature). Such condition is known as isothermal-isobaric ensemble. Energy minimization data and output of MD simulation data helped us to visualize several thermodynamic parameters like potential energy, kinetic energy, total energy, temperature, density progression, radius of gyration, RMSD, RMSF, SASA, etc. of the bare host protein and ligand: receptor complex. For equilibration step computation, for covalent bond constraints LINCS algorithm was applied. To quantify strength of interaction between ligand and protein it is useful to compute the nonbonded interaction energy between two species rather than calculating free energy of the system. For the nonbonded energy calculation we have used short range Lennard-Jones and Coulomb interaction energies. For these calculation we have maintained a 1.4 nm radius cut-off. For long range electrostatics calculation Particle Mesh Ewald (PME) method was used having 1.6 Å of Fourier grid spacing. Berendsen temperature coupling method (V-rescale) has been used for maintaining inside box temperature. NPT equilibration has been achieved by Parrinello-Rahman pressure coupling method. The time parameter has been fixed to 10000 ps with each step of 0.002 fs for the final step of simulation. After final step of each simulation, trajectories were obtained. Obtained trajectories and results were analysed using the graphical tool Origin pro. Least-square fitting method was used to evaluate RMSD for protein backbone. Similarly, RMSF was obtained for protein Cα backbone. The compactness factor of protein was calculated by radius of gyration (R). According to rule if the protein or complex structure is stable then the radius of gyration should maintain a stable value. To compute total solvent accessible surface area (SASA) the tool sasa was used. SASA measures the area of receptor exposure to the solvents during the simulation process. Number of hydrogen bonds and distribution of intermolecular hydrogen bond lengths were calculated with maintaining 3.5 Å distance cut-off condition throughout the simulation.

Computational details

Dell Gen9 server with 8 Core I7 processors and 16 GB of RAM with GPU NVIDIA MX130 was used for the MD simulations and corresponding energy calculations.

Results and discussion

Virtual screening and molecular docking analysis

Each expected ligand drug molecule was docked to CoV- 2 3CL maiene">n protease (n class="Chemical">3CLPRO). We have to keep in mind that docking algorithm is nondeterministic in nature. The minimum of scoring function cannot be identified by docking algorithm for a fixed receptor: ligand complex conformation. Table SD1 shows the structure of possible ligands found in the active site pockets of 3CLPRO proteins I and II. Molecular docking is used to find out interaction of tested inhibitors: berberine, β-sitosterol, choline, tetrahydropalmatine and octacosanol with both 3CLPRO proteins I and II. The organic compound berberine (C20H18NO4) can be extracted from Hydrastis canadensis L., Berberidaceae. It is also known as protoberberine alkaloids. Berberine is usually found in the roots, rhizomes, stems, and bark. One of the main source of berberine is Tinospora cordifolia. It can be found in other plants also. Berberine is used as a natural dye with a color index of 75160. Due to its strong yellow fluorescence, it is useful in histology for staining heparin in mast cells. Berberine is relatively toxic parenterally, but still it is used as oral drug for the treatment of various parasitic and fungal infections. It is also used a drug for antidiarrheal diseases (Dewick, 2009). β-sitosterol (C29H50O) has a structural similarity with cholesterol. It is used to reduce cholesterol levels in the body. It is used as a constituents in many drugs for reduction of swelling. Being a steroid, β-sitosterol is a precursor of anabolic steroid boldenone (Wang et al., 2008). Octacosanol (C28H58O) is basically an alcohol “straight-chain aliphatic 28-carbon primary fatty alcohol” which is used health industry as a nutritional supplement. It has been studied as a potential therapeutic agent for the treatment of Parkinson's disease (Rudkowska et al., 2008). Octacosanol is reported to possess cholesterol-lowering effects, antiaggregatory properties and ergogenic properties. Tetrahydropalmatine (C21H25NO4) has a role as an adrenergic agent, a non-narcotic analgesic and a dopaminergic antagonist. Tetrahydropalmatine is under investigation in clinical trial for the Treatment of Schizophrenia (Taylor et al., 2003). Choline (C5H14NO) is a nutrient that supports various bodily functions, including cellular growth and metabolism. It is also used for DNA synthesis and nervous system maintenance (Corbin & Zeisel, 2012).

Candidate drugs that confirm to the Ro5 and other drug likeness rules etc. used to have lower attrition rates duriene">ng cliene">nical trials aene">nd hence have showene">n their stroene">ng caene">ndidature as potential drugs reachiene">ng the market (Table 1a–e). n class="Chemical">Berberine, choline and tetrahydropalmatine were satisfying the most prescribed virtual screening properties. Auto Dock Vina, the docking software uses free-energy scoring function based on a linear regression analysis to identify the best binding mode prediction for protein: ligand complexation. It identifies the best possible pose structure from a set of diverse protein-ligand complexes with a highest negative value of binding energy and lowest value of inhibition constant. Docking simulation results have helped us to find out the best possible interacted protein: ligand poses of different proposed inhibitor drugs. For the first tested inhibitor berberine with 3CLpro protein I (6LU7), berberine: 6LU7 complex revealed 9 different poses based on their binding affinity, drieding energy, dipole moment, inhibition constant etc. According to molecular docking simulation protocol the pose with highest negative values of binding energy, lowest value of complex driending energy usually indicates maximum binding affinity for the ligand: receptor complex formation (Table 2).

Table 2.

Binding mode of each ligand:protein (Berberine: 6LU7) complex using molecular docking for Protein (I).

Pose	Binding affinity (G) (kcal/mol)	Hydrogen bonded interaction number	Drieding energy (ligand)	Drieding energy (protein)	Dipolemoment of ligand(debye)
1	−7.3	0	414.419	6,044.79	1.676
2	−7.3	2	415.725	6,045.91	1.677
3	−7.1	0	412.076	6,042.69	1.675
4	−7.0	0	414.361	6,044.95	1.68
5	−7.0	0	415.549	6,045.57	1.672
6	−6.9	0	414.623	6,044.94	1.676
7	−6.9	0	414.409	6,044.77	1.68
8	−6.9	0	415.07	6,046.38	1.68
9	−6.8	0	415.635	6,046.05	1.675

Binding mode of each ligand:protein (Berberine: n class="Chemical">6LU7) complex using molecular docking for Protein (I).

For pose 2 we obtained the better interacted position for ligand: protein complex with the binding affinity of −7.3 kcal/mol which is lower than highly tested COVID 19 drug n class="Chemical">chloroquine (−6.29 kcal/mol). To verify that this the best docked site we have also computed the Dreiding energy of different poses which consider the Dreiding force field to calculate the energy of a specific complex structure by summing energy components like bond lengths, bond angles, dihedral angles (Mayo et al., 1990). The Dreiding energy values for individual protein and ligand are different but when they form complex the Dreiding energy becomes minimum for most favourable structure. We obtained the lowest value of energy (415.725) for the pose 2 of the complex structure (Table 3). To validate the better interaction between berberine and 6LU7, we have computed the inhibition constant (ki) which is an indication of how potent berberine as inhibitor towards 6LU7. It is the concentration to produce half maximum inhibition. To evaluate the inhibition constant we have used relation as where G is the binding affinity, R the universal constant and T (298 K) the room temperature.

Table 3.

Interaction detail for different ligands: berberine, β-sitosterol, choline, tetrahydropalmatine and octacosanol with receptor protein I (6LU7).

Ligand	Best Binding affinity (kcal/mole)	Hydrogen bonded interaction (protein donor: ligand acceptor, distance in Å)	Hydrophobic interaction (protein donor: ligand acceptor, distance in Å)	Dipole moment of ligand (debye)	Drieding energy between protein and ligand	Inhibition constant (M)
berberine	−7.3	(A:THR25:HG1 - :UNK0:O, 2.11198)(A:SER46:HG - :UNK0:O, 2.88545)	(:UNK0 - A:HIS163, 3.61032)	1.67	415.725	4.4 × 10−6
β-sitosterol	−7.1	(:UNK0:O - A:ARG188:O, 3.12606)	(:MET165 -:UNK0, 4.54304)	1.761	1280.56	6.16 × 10−6
choline	−3.4	(:UNK0:O - A:TYR54:OH, 2.8027)(:UNK0:C - A:MET49:O, 3.71379)(:UNK0:C- A:GLN189:OE1, 3.61393)	(:UNK0:C - A:HIS41, 3.8048)	5.374	104.198	3.2 × 10−3
tetrahydropalmatine	−6.4		(:UNK0:C - A:MET165, 4.73002)(:UNK0:C - A:LEU167, 4.89871)(:UNK0:C - A:PRO168, 5.01258)(:UNK0 - A:CYS145, 5.27402)(:UNK0 - A:MET165, 4.58548)	2.688	373.677	2.01 × 10−5
octacosanol	−6.6	(:UNK0:O - A:LEU141:O, 3.37559)(:UNK0:C - A:HIS163:NE2, 3.62548)	(:UNK0 - A:ALA191, 4.8791)	1.257	414.094	1.43 × 10−5

Interaction detail for different ligands: berberine, β-sitosterol, n class="Chemical">choline, tetrahydropalmatine and octacosanol with receptor protein I (6LU7).

For the best pose (2) structure of berberine: n class="Chemical">6LU7 complex, the ki value obtained at room temperature (298 K) as 4.4 x10−6 M which proves the higher affinity of berberine towards receptor 6LU7 (Table 3). The computed ki value is far lower than the toxicity dose range of berberine as proposed drug (Table 1) which validate the strong candidature of berberine as target drug to attach 6LU7. The berberine: 6LU7 interaction energy further verified by the dipole moment values of ligand and target (Table 3). The strong interaction for pose 2 was further verified by the number of weak nonbonded hydrogen bonded interactions and hydrophobic interactions present between protein and ligand in the optimized protein: ligand complex structure. Hydrogen bonded interaction and hydrophobic interaction due to their weak nature have significant roles in stabilizing the energetically-favored ligands in a suitable pocket of the environment of protein structures. It is a proven fact that weak interactions “hydrogen bonding and hydrophobic interactions” always stabilize the ligands at the target protein site by altering the binding affinity (Patil et al., 2010). We observed the presence of intermolecular hydrogen bonds and hydrophobic interaction between protein (Residues: THR25, SER46, HIS163) and berberine in its pose 2 which means better interaction between donor and acceptor moiety. We have repeated the molecular docking simulation of berberine with other 3Lpro protein II also and obtained the best ligand: protein complex structure with binding affinity of −7.7 kcal/mol and inhibition constant of 2.23 × 10−6M (Table SD3). For best poses of berberine: protein complexs, the donor–acceptor surface and different possible interactions are shown in Figure 4(a, b) in 3D and 2D view.

Figure 4.

Donor: acceptor surface and possible types of interactions in best pose structure obtained from molecular docking for berberine: 6LU7 and berberine: 6W63.

Donor: acceptor surface and possible types of interactions in best pose structure obtained from molecular docking for berberine: n class="Chemical">6LU7 and berberine: 6W63.

We have repeated same molecular docking approach for other ligand structures: β-sitosterol, octacosanol, n class="Chemical">Tetrahydropalmatine, Choline with protein I, protein II and have identified their best possible ligand: protein interaction pose position in terms of their best binding affinity value, Dreiding energy, dipole moment, inhibition constants, number of hydrogen bond, hydrophobic bond etc. and mentioned them all in Tables 3, SD3 and Figure 5(a, b) and Figure SD2 in 3D and 2D view. From the molecular docking results for every suggested drug molecule, with their best possible ligand: target complex structures i.e. for best pose, the donor–acceptor surface with their possible hydrogen bonding and hydrophobic interactions are shown in Figure (SD2a, b) in 3D and 2D view. Our screening identified that out of five possible bioactive ligand structures, berberine shows the best potentiality to inhibit with the CoV-2 enzymes (3Lpro): I and II by its best docking affinity compared to the other ligands against 3Lpro. Berberine: 3Lpro complex has shown less binding energy, minimum inhibition constant value as compared to other proposed drug molecules with its good binding mode of interactions. After berberine the binding affinity followed by β-sitosterol, octacosanol, Tetrahydropalmatine, Choline. Simulation results also revealed that all five proposed drug molecules showed good stability as a complex with the targeted 3Lpro protease. These inhibitors also satisfy the required drug likeness properties according to Ro5, Veber etc rules including their molar refractivities, pharmakokinetics, polar surface areas and logP values. All of the reported potential natural phytoconstituent drug compounds are commercially available in drug industry and so can be easily available for further in vivo/in vitro validations. We are hopeful that the information generated from this present study will definitely be utilized for the development of phytochemical based therapeutics against COVID-19. Since berberine shows the better possibility of inhibition towards 3CLpro protease we have further studied the applicability of berberine as potential drug molecule by applying rigorous MD simulation approach.

Donor: acceptor surface and possible types of interactions in best pose structures obtained from molecular docking for different ligands with COVID-19 protease eene">nzyme n class="Chemical">3CLpro (6LU7).

MD Simulation analysis

For obtaining dynamic data at atomic spatial resolution, MD simulation is a verified in silico method which can simulate in picoseconds/nanoseconds or further finer temporal steadfastness (Benson & Daggett, 2012; Gajula et al., 2016). The complex (berberine: n class="Chemical">6LU7) formed by main protein protease (6LU7) docked with phytochemical compound (berberine) has been studied under this simulation for 10000 ps to 100 ps period to analyze the stability of the studied structure. For MD simulation to work every structure (receptor and complex) have to be optimized by energy minimization process. For energetically minimized structure the potential energy should be minimum and negative with a maximum force value. We have obtained energetically minimized structures for both probe (protein and complex) systems (Figure 6). For both the cases we have obtained steady convergence of potential energy. The analyzation and comparison of the potential energy (Epot) of the stabilized structure of 6LU7 protein in its bare state and in berberine: 6LU7 complex form have been done rigorously. 6LU7 has an average energy of −1.27 × 106 ±56.7 (kJ mol−1) for Epot, while the complex has an averaged Epot of −0.2 × 106 ± 38.77 (kJ mol−1). With their lowest Epot values both systems at an energy minimum level and ready real Molecular Dynamics. To check further, optimized structures of 3CLpro protease and complex structures have been equilibrated by NVT and NPT ensembles to stabilized temperature, pressure, density, volume etc within a time scale trajectory of between 100 ps and 10000 ps. Analysing the temperature progression data (Figure SD4) with an running average value of 10 ps, it was observed that after starting the simulation process, the temperature of the system quickly reached the stable room temperature value (298 K) with maintaining the temperature stability throughout the whole equilibration process. This is the reason for bare protein system, a shorter equilibration period (on the order of 50 ps/100 ps) have been applied. The applicability of shorter equilibration period have been also verified by pressure and density data (Figure SD5) with same running average value of 10 ps. Also we have observed a very stable values of density and pressure over the period of time trajectory. From the stability of temperature, pressure and density values it was justified to conclude that the system is well-equilibrated with respect to temperature, pressure and density and are ready for MD simulation. We have studied various thermodynamic parameters of host protein and complex system by MD simulation data for various time trajectory. The important thermodynamic parameters we have studied are the RMSD, RMSF, inter-molecular H- bonds, R, SASA, potential energy, kinetic energy, total energy and various nonbonded interaction energies to understand the stability and possible conformational changes of the bare state of CLpro receptor protein and in the berberine: 6LU7 complex at each frame of the time resolved simulation trajectory. All MD simulated output of bare protein and protein: ligand complex for the time resolved trajectories have been presented in Table 4.

Figure 6.

Optimized geometry of receptor 6LU7 and berberine: 6LU7 complex.

Table 4.

MD simulation output of time resolved trajectory of 6LU7 in its bare state without any ligand and in the complex state with ligand berberine.

Serial no	Parameter	Bare 3CLpro protease (6LU7)		6LU7 + berberine complex
		Mean	Range	Mean	Range
1.	SR Columbic Interaction energy (kJ mol−1)	NA	NA	−51.7714 ± 3.8	20–−80
2.	SR LJ Interaction energy (kJ mol−1)	NA	NA	−130.261 ± 1.1	−80–−160
3.	Average Interaction Energy energy (kJ mol−1)	NA	NA	−186.95 ± 5.1
4.	RMSD (nm)	0.12	0.08–0.16	0.16	0.08–0.18
5.	Inter H-Bonds	NA	NA	0.026	0–2
	Intra H-Bonds	630	614–653	213	196–220
6.	Radius of gyration	2.25 ± 0.01	2.25–2.26	2.25 ± 0.01	2.25–2.26
7.	SASA (nm2)	22	19–26	148	146–156
9.	Potential Energy (kJ mol−1)	−1.27 × 106 ± 56.7	−7.3 × 105–−1.3 × 106	−0.2 × 106 ± 38.77	−2.4 × 104–−0.2 × 106
10.	Total Energy (kJ mol−1)	−1.17 × 106 ± 0.01 × 106	−1.3 × 106–−1.0 × 106	−1.8 × 10 5 ± 0.01 × 105	−1.8 × 105–−1.7 × 105

Optimized geometry of receptor 6LU7 aene">nd n class="Chemical">berberine: 6LU7 complex.

MD simulation output of time resolved trajectory of 6LU7 iene">n its bare state without aene">ny ligaene">nd aene">nd iene">n the complex state with ligaene">nd n class="Chemical">berberine.

Throughout the whole time resolved trajectory the total energy of both bare protein (6LU7) aene">nd n class="Chemical">berberine: 6LU7 complex remain stable having average value of −1.8 × 10 5 kJ mol−1 with a minor fluctuations between −1.8 × 10 5 kJ mol−1 and −1.7 × 105 kJ mol−1 which validate the good stability of both bare protein and its complex structures for the whole time trajectory (Figure SD 6,7). The variation of total energy of the complex due to its kinetic energy and potential energy components is shown in Figure 7.

Figure 7.

Total energy of berberine: 6LU7 Complex structure in time scale 1000 ps.

Total energy of berberine: n class="Chemical">6LU7 Complex structure in time scale 1000 ps.

Protein backbone RMSD (Figure 8) of bare protein was observed to be fluctuating between 0.08 and 0.16 nm with an average value of 0.12 nm. In presence of compound berberine after complex formatioene">n n class="Chemical">berberine: 6LU7, RMSD showed a stable and constant range between 0.08 nm to 0.18 nm with an average value of 0.16 nm. The obtained average RMSD for complex is comparatively slightly higher than the RMSD of bare protein backbone. The complex RMSD data and bare background protein have been compared with the crystal structures for both the configuration. With full time scale range the obtained RMSD variation levels up to ∼0.1 nm for both bare protein and protein: ligand complex indicating the stability of both the structures. Subtle differences between the plots of gaseous structure and crystal structure indicated minor structural difference at t = 0 ps. This minor difference in structures is very much expected since for energy-minimized structures position restraints can never be 100% perfect. Comparing the differences of protein backbone RMSDs for bare CLpro protein structure with berberine: protein complex structure suggest that CLpro has shown no significant change in presence of berberine during time resolved simulation. Whereas the RMSD of complex (berberine: 6LU7) has indicated several fluctuations at different time intervals on 10000 ps time scale. These fluctuations may arise due to the changed conformation in binding region of CLpro (Figure SD8). In the full time scale intervals, the first stable structure was observed between 1645 ps to 2860 ps, second one is observed between 3574 ps to 6372 ps and the third one is between 6523 ps to 8753 ps. Throughout the whole time trajectory, in the stable region the RMSD remain constant at 0.5 nm whereas for large fluctuation regions the observed RMSD was more than 0.6 nm. The fluctuations were observed between 0 ps to 1644 ps, 2861 ps to 3573 ps and 6373 ps to 1000 ps.

Figure 8.

Root mean square Deviation (RMSD) backbone graphs of 6LU7 in its bare state and in complex with berberine. Complex structure is compared with its crystal form also.

Root mean square Deviation (RMSD) backbone graphs of 6LU7 iene">n its bare state aene">nd iene">n complex with n class="Chemical">berberine. Complex structure is compared with its crystal form also.

The above observation also confirmed by the RMSF plots of Cα for both the bare protein and its complex where local changes were observed with variation of residues numbers (Figure 9). For bare protein 6LU7, RMSF value for each residue was simulated aene">nd averaged for the entire time trajectory aene">nd identificatioene">n has been doene">ne for more flexible regioene">ns respoene">nsible for the coene">nformatioene">nal chaene">nges. In simulated RMSF graph of bare n class="Chemical">6LU7, major fluctuations in peaks have been observed with the residues positions between 50 and 80 and the tail region of the protein with the residue number around 280. The fluctuations in RMSF values were observed between 0.13 nm and 0.23 nm. For berberine: 6LU7 complex, the major fluctuations in peaks were observed initially around residue number at 50 and in between140–200 values, with fluctuations around 0.23 nm. The RMSF data suggest that all active and flexible residues of 6LU7 get stabilized due to binding of berberine for the complex formation, which hints the possibility of inhibitory activity of the ligand towards receptor protein without affecting the protein structure. From the detailed RMSD and RMSF analysis we can conclude that protein backbone of 6LU7 in presence of berberine was not affected much and that why the behaviour of protein or structure of protein during simulation remain mostly unaltered (Figure 9).

Figure 9.

Graphs of Root mean square fluctuations (RMSF) of 6LU7 in its bare state and in berberine: 6LU7 complex.

Graphs of Root mean square fluctuations (RMSF) of 6LU7 iene">n its bare state aene">nd iene">n n class="Chemical">berberine: 6LU7 complex.

To get an idea about the binding strength in between ligand and protein for protein: ligand complex, hydrogen boene">ndiene">ng plays a crucial role. Ligaene">nd n class="Chemical">berberine has constant range of intermolecular hydrogen bonding with receptor protein between 0 to 2 in throughout the whole simulation process with an average value ±0.26. This result validate the nonexistence of any type of conformational change around berberine during the complexation in the binding site throughout simulation process (Figure 10). The intermolecular hydrogen bond number computed through MD simulation over the full time scale trajectory perfectly matches with the docking results also. Over all observation suggested that the berberine: 6LU7 complex is stable during simulation.

Figure 10.

Intermolecular hydrogen bond numbers between berberine and 6LU7 in their complex (berberine: 6LU7) form for the full scale time trajectory.

Intermolecular hydrogen boene">nd numbers between n class="Chemical">berberine and 6LU7 in their complex (berberine: 6LU7) form for the full scale time trajectory.

The radius of gyratioene">n (Rg) measurement of a proteiene">n ligaene">nd complex or basic backboene">ne proteiene">n are used to determiene">ne the compactness of the system with the time. For stably folded proteiene">n structures, R value always like to maiene">ntaiene">n a relatively steady value for full time scale. For unfolded proteiene">ns, R always chaene">nges over the time. Agaiene">n for less compactness iene">n the structures with coene">nformatioene">nal entropy exhibits higher R values, while high compactness with more stability iene">n the structure exhibits a low R value. For the present case, R values of both complex aene">nd base proteiene">n structures are reported iene">n between 2.25 aene">nd 2.26 nm with aene">n average value of 2.225 nm (Figure 11) which exhibited the stability of both bare proteiene">n aene">nd proteiene">n: ligaene">nd complex form. Small variatioene">n iene">n R values also revealed that the biene">ndiene">ng of n class="Chemical">berberine towards binding site of 6LU7 does not induce any change in the parent CLpro protein structure and establish the condensed stable architecture, size of the berberine: 6LU7 complex form (Figure 10). For the current paper we have shown the values over the time scale of 100 ps. The three major anti COVID-19 drugs Remdesivir Saquinavir and Darunavir also reported as showing the average R score values as 2.2 ± 0.1 nm (Khan, Zia, et al., 2020) which are exactly matching with our computed data for our proposed drug berberine.

Figure 11.

Total radius of gyration of bare 6LU7 and berberine:6LU7 for the time trajectory 1000 ps.

Total radius of gyration of bare n class="Chemical">6LU7 and berberine:6LU7 for the time trajectory 1000 ps.

Further to validate the applicability of berberine as proposed drug for n class="Disease">COVID-19 we have computed the SASA of our proposed berberine: 6LU7 complex structure. SASA measures the area of receptor exposure to the solvents during the simulation process. We obtained the SASA value of receptor protein between 19 − 26 nm2. The exposures of the hydrophobic part of receptor residues due to the binding of the ligand molecule towards receptor always contribute to SASA value. In the present work for both bare 6LU7 and berberine: 6LU7 complex, the computed values of SASA values were observed between 146 and156 nm2 which justified that ligand binding does not affect the folding of the receptor protein very much. So we may conclude that berberine: 6LU7 complex is very stable structure due to the binding of berberine with 6LU7 CLpro protein.

To quantify the strength of the interaction between berberine aene">nd n class="Chemical">CLpro protease, it is always useful to compute the short-range nonbonded interaction energy between these two species which is far better than free energy or a binding energy. For protein ligand system, for free or binding energy simulation techniques, the mostly used force fields are not parametrized and so the obtained energies are in most cases become actually physically not meaningful. Whereas for nonbonded interaction energy simulation, the used force fields are parametrized to specifically target quantum mechanical interaction energies with water, so it is intrinsically balanced against meaningful quantities, and so the nonbonded interaction energy shows meaningful interpretation. For interaction energy computation we are interested in Coulombic short range protein ligand interaction energy terms and Lennard Jones short range protein ligand interaction energy terms. The total interaction energy is obtained by using these two energy values.

The obtained total interaction energy for the berberine: n class="Chemical">6LU7 complex is shown by the 3D Figure 12. We have obtained total interaction energy for the berberine: 6LU7 conformation in every 1000 ps time gap. The average short-range Coulombic interaction energy component compared with complex 51.7714 ± 3.8 kJ mol−1 has less effect on binding affinity, whereas the average short-range Lennard-Jones energy is −130.261 ± 1.1 kJ mol−1 has shown strong effect on binding affinity. The total interaction energy after propagating the error according to the standard formula for addition of two quantities appeared as −186.95 ± 5.1 kJ mol−1. The low region of the 3D view of the time trajectory of interaction energy perfectly matches with the RMSD stable regions of the complex form and high region matches with the fluctuations of the RMSD time scale trajectory regions. All the above mentioned MD simulations data validated the good inhibition possibility of berberine towards the target CLpro protease 6LU7.

Figure 12.

Total nonbonded interaction energy between 6LU7 and berberine in their complex form for the time trajectory 1000 ps. Total energy is the combination of Coulombic and Lennard Jones interaction energies.

Total nonbonded interaction energy between 6LU7 aene">nd n class="Chemical">berberine in their complex form for the time trajectory 1000 ps. Total energy is the combination of Coulombic and Lennard Jones interaction energies.

Conclusion

Traditional medicinal plants are considered as a great source of input for the treatment of various diseases. For current study, we have studied the applicability of antiviral potential of some phytoconstituents extracted from a specific medicinal plant: Tinospora cordifolia agaiene">nst n class="Disease">SARS-CoV-2 infections. Tinospora cordifolia is the source of various type of bioactive compounds including alkaloids, steroids, glycosides, aliphatics. Spotlight of the current study is to find an essential drug for the COVID19 disease using the antiviral activity of these compounds. Summarizing the docking and MD simulation results and interpreting the various thermodynamics parameters like temperature, pressure, density, root mean square deviation, root mean square fluctuations, inter-molecular hydrogen bonds, solvent accessible surface area, radius of gyration, potential energy, kinetic energy, total energy and interaction energies for non bonded interactions like hydrogen bonding, hydrophobic interactions, short range interactions with their analysis plots helps us to identify the best inhibitor potential drug for SERS CoV-2. The best binding docking poses of virtually screened phytochemicals berberine, β-sitosterol, octacosanol, tetrahydropalmatine, choline from the parent herb Tinospora cordifolia with that of the selected 3CLpro targets I, II of main protease enzymes suggest the promising potential of these molecules to be used as raw drug material. Molecular docking and MD analysis have revealed that berberine having less binding energy and higher nonbonded interaction capability as compared to other molecules with its good binding mode of interactions. Simulation results have also revealed that berberine: 3CLpro docked complex established better stability and can act as better inhibitor towards the CoV-2 protein compared to other inhibitors. All the reported phytochemical compounds in the present work are natural and commercially available as drugs. The present work suggests that among all the proposed phytochemicals extracted from Tinospora cordifolia, berberine has established its strong candidature to serve as potential inhibitors in regulating the 3CLpro protein's function and further controlling against viral replication. Our results are expected to grab the attention of the researchers in the field of new drug discovery against SERS-CoV-2 for which till date no specific drugs or vaccines are available. The approach adopted here is general in nature and similar approach may be used to investigate the potential applications of other medicinal herbs and available drugs against COVID-19. In addition, the study may be extend to more precise investigation of protein-drug interaction using different quantum mechanical simulation methods like TDDFT, HF with modern tools like ONIOM which will help us to identify optimized and energetically favored drug:target-protein complex structures. However, such a task is computationally demanding and requires more computational resources and time. In a future communication, we would like to extend our work in that direction, so that if a clinical trial of the drug molecules studied here is aimed the same can be started with more confidence. Finally, we conclude this paper with a word caution that before using any outcome of an in silico study, proper in-vivo and in-vitro rigorous research works are to be performed.

Click here for additional data file.

Click here for additional data file.