Curcumin to inhibit binding of spike glycoprotein to ACE2 receptors: computational modelling, simulations, and ADMET studies to explore curcuminoids against novel SARS-CoV-2 targets

The recent emergence of the novel coronavirus (SARS-CoV-2) has raised global concern as it is declared a pandemic by the WHO. However, to date, there is no current regimen to mitigate the molecular pathogenesis of SARS-CoV-2 virus. Curcuminoids, bioactive ingredients present in Curcuma longa (turmeric), are known to exhibit diverse pharmacological properties. To the best of our understanding to date, SARS-CoV-2 uses angiotensin-converting enzyme 2 (ACE2) for the host cellular entry. This is mediated via proteins of SARS-CoV-2, especially the spike glycoprotein receptor binding domain. Accordingly, our primary objective is to thwart virus replication and binding to the host system, leading us to probe curcuminoids efficiency towards key surface drug target proteins using the computational biology paradigm approach. Specifically, fourteen natural curcuminoids were studied for their possibility of inhibiting SARS-CoV-2. We studied their in silico properties towards SARS-CoV-2 target proteins by homology modelling, ADME, drug-likeness, toxicity predictions, docking molecular dynamics simulations and MM-PBSA free energy estimation. Among the curcuminoids docked to the receptor binding domain of SARS-CoV-2 spike glycoprotein, the keto and enol forms of curcumin form strong hydrogen bond interactions with ACE2 binding residues Q493, T501, Y505, Y489 and Q498. Molecular dynamics simulations, free energy binding and interaction energy validated the interaction and stability of the docked keto and enol forms of curcumin.


Introduction
The world has recently witnessed an outbreak of a potentially lethal coronavirus, nCovid-19, an urgent public health issue with an increase in morbidity and mortality. Coronaviruses belonging to the family Coronaviridae signicantly threaten human health, and other species. Recent ndings claim that the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is evolutionarily closely related to bat coronavirus and SARS-CoV with a homology of $95% and $70%, respectively. 1-3 SARS-CoV-2 is a rapidly transmissible virus that tends to change its genetic material to survive in various environmental conditions. Structurally, this virus contains a nucleocapsid (N), where its genome is packed inside a helical capsid. A membrane protein (M) and small envelope protein play a major role in the virus assembly, and the spike glycoprotein or S protein plays a key role in host entry by harboring certain crucial amino acid residues of human angiotensin-converting enzyme 2 (hACE2). The spike glycoprotein forms a large protrusion on the surface of the virus, exhibiting a crown-like appearance, and hence the name coronaviruses, 4-8 as shown in Fig. 1.
Similar to SARS-CoV, SARS-CoV-2 also facilitates its entry into human angiotensin-converting enzyme 2 (ACE-2) through the receptor-binding domain (RBD) of the spike glycoprotein. Angiotensin-converting enzyme 2 (ACE-2), is an enzyme located SARS-CoV-2 spike glycoprotein RBD is considered a new therapeutic intervention. 11 However, in the current pandemic situation, traditional drug discovery or vaccine development is a daunting and time-consuming task, which can be offset by computational-aided drug design.
Moreover, instead of designing new lead compounds, drug repurposing or the use of natural products can be an alternative for the development of new antiviral compounds, which will enhance the speed of research in this area. In the current study, we focused on the natural plant Curcuma longa, commonly known as turmeric, a perennial herbaceous rhizomatous plant belonging to the ginger family Zingiberaceae, which is widely used in India. 12 Curcuma longa chemical constitutes are widely used for treating various ailments and possess a wide variety of therapeutic properties including antiviral, 13 analgesic, 14 antimicrobial, 15 antiproliferative, 16 and anti-inammatory 17 activity (Fig. 2).
Among the various chemical compounds, curcumin has gained importance among researchers because its compounds have been exhibited activity against viruses such as the human immunodeciency virus (HIV), dengue virus, herpes simplex virus (HSV), hepatitis virus, inuenza A virus (IAV), and Ebola virus. [18][19][20][21] Recently, reports on computational drug design demonstrated the therapeutic potential of curcumin as a dual inhibitory agent acting on S-protein and ACE-2. 22 Thus, the main rationale of this study was to search for potential curcumin derivatives against these two drug targets, i.e. the spike glycoprotein RDB region and envelope protein. Besides, the pharmacokinetics, pharmacodynamics and drug-likeness of all the curcumin derivatives were studied. Finally, a molecular dynamics simulation study was performed to understand the stability of the protein-ligand complex with time, and a CHARMm interaction analysis and free energy binding using MM-PBSA demonstrated the stability of curcumin complexed with the spike glycoprotein.

Homology modelling
Homology models are the backbone of structural biology, where a structure without NMR or X-ray crystallography data can be modelled using various machine learning algorithms. Moreover, it is a reliable and widely used technique for predicting unknown protein structures. The experimental protein structures of the novel coronavirus SARS-CoV-2 belonging to the genus Betacoronavirus are limited in structural databases. Therefore, the template modelling concept was implemented tsuences were scanned in the TMpred tool to identify the TMalpha helices among the three query sequences, and only the envelope protein showed transmembrane helices (Fig. 3A). This journal is © The Royal Society of Chemistry 2020 RSC Adv., 2020, 10, 31385-31399 | 31387 We searched for templates for each query sequence using BLAST search against PDB databases. According to the BLAST search of the spike glycoprotein_RBD (333-526) region, it shows 100% template coverage with the RBD region of SARS-CoV-2 (PDB ID 6M0J_E) with a resolution of 2.45 A. Also, the homology of bat coronavirus RaTG13 (PDB ID 6ZGF_A), SARS-CoV protein (PDB ID 3SCI_E), SARS-CoV BJ01 (PDB ID 5X58_A) and SARS-CoV-2 is about 90.13%, 73.54%, and 73.09% identity with query coverage of 100%, respectively. In contrast, the homology of the MERS S-protein (PDB ID 4L72_B) and SARS-CoV-2 S-protein RBD is about 24% identity and the query coverage rate is 44%. The BLAST search results of the query sequence to the PDB structure alignment is depicted in Fig. 4.
Subsequently, the BLAST search of the envelope protein (QHS34548) of SARS-CoV-2 showed the homology of SARS-CoV PDB ID's 2MM4_A and 5X29_A. SARS-CoV envelope protein 5X29_A shows maximum query coverage of 82% with 88.71% identity and 2MM4_A of SARS-CoV sequence alignment with a query of only about 77% with 91.38% identity. However, the aligning query sequence with individual templates of SARS-CoV using the MODELLER algorithm showed 83% sequence similarity and 81.5% sequence identity (Fig. 3B). Hence, template 5X29_A was used for homology modelling of the SARS-CoV-2 envelope protein. However, no template match was found with MERS-CoV. On the other hand, the membrane protein sequence for template structure identication showed no signicant similarity results with the default scoring matrix of BLASTp. Consequently, other scoring matrices were introduced to check the structural similarity of the sequence, excluding an iterative process that showed no signicant similarity to the query. Hence, the modelled envelope protein (Fig. 3C) and available X-ray crystallography structure 6M0J_E spike glyco-protein_RBD of SARS-CoV-2 were used for further studies.

Structure validation process
The modelled envelope protein structure from the SARS-CoV structural parameters was studied for validation purposes. The best-modeled structure using MODELLER was initially chosen based on the lower discrete optimized protein energy (DOPE) 23 of À4812.83 kcal mol À1 , and the higher total PDF (probability density function) of 165.148 shows 92.98 ERRAT quality factor. The Procheck tool evaluates the residue-by-residue stereochemical quality, structure geometry and distribution of Phi and Psi angle amino acids in favoured, allowed and generously allowed region of a modelled protein. 24 The Ramachandran plot (RP) showed 91.8% amino acids in the core region or most favoured regions and 8.2% residues were in the additional allowed regions with no outliers (Fig. 5A). Most of the amino acids in the favoured region signify that the modelled structure of the envelope protein is reliable, and can be equally compared with NMR structure quality. QMEANBrane is a unique structure quality assessing SWISS-MODEL tool explicitly used for modelled transmembrane protein investigation. 25 The envelope protein probed in SWISS- MODEL was within the transmembrane insertion energy, and also satised statistical potentials in the naturally occurring oligomeric state (Fig. 5B). The template (5X29) and modelled structure superimposition are shown in Fig. 5C, where the overall structure root mean square deviation is <2.5 A. The quality of the structure was further conrmed using the RAMPAGE tool. In  This journal is © The Royal Society of Chemistry 2020 RSC Adv., 2020, 10, 31385-31399 | 31389 addition, ProSA structure analysis 26 was performed for the modelled protein and template structure, where the Z-scores of 0.54 and 0.59, respectively, located in the NMR structure region further conrm the modelled protein structure quality (Fig. 6).

Computational pharmacoinformatic study on curcumin derivatives
ADMET, TOPKAT and drug-likeness assessment for lead molecules can provide insight into their quantitative structureproperty relationship. Thus, the models used in the in silico pharmacokinetic study were obtained through a quantitative structure-activity relationship (QSAR). Furthermore, they reduce the research time, biological waste and cost. ADMET was initially used to a screen a large library of compounds. Descriptors such as the blood-brain barrier (BBB) penetration, hepatotoxicity and CYP2D6 enzyme are high priority models followed by Ames mutagenicity and carcinogenicity. All the curcuminoid compounds were tested using the ADMET model. The pharmacokinetic analysis results are shown in (Fig. 7). According to the ADMET plot Alog p_98 vs. PSA, it can be  observed that compound 71315012 (curcumin glucuronide) is outside of the ellipses due to its low human intestinal absorption (HIA) and undened BBB penetration, whereas some molecules such as 5318039 (hexahydrocurcumin), 11068834 (octahydrocurcumin) and hexahydrocurcumin are near to 99% condence ellipses of the BBB, and the remaining compounds are inside 95% and 99% condence ellipses of HIA. Also, the curcuminoids compounds exhibit good membrane permeability, except curcumin glucuronide, which has a polar surface area 27 (PSA) of <140 A 2 . Overall, the curcuminoids compounds are non-toxic to hepatic cells, non-inhibitors of a metabolic enzyme (CPYD26), and exhibit very high to medium penetration across the BBB except a few molecules (Table 1). Similarly, the in silico pharmacodynamic models of the curcumin derivatives are free from carcinogens and mutagens. Quantitative-structure toxicity relationship (QSTR)-based toxicity parameters such as rat (TD50), rat_oral (LD50), rat_inhalation, fathead minnow (LC50) and Daphnia (EC50) are summarized in Table 2. With only one RO5 violation allowed, Alog P of #5, molecular weight of #500, hydrogen bond acceptors (HBA) of #10 and hydrogen bond donor (HBD) of #5 are also known as Lipinski's rule of 5 or Pzer's rule of ve, which is used to evaluate the drug-likeness property of lead molecules. 28 Compounds that conform to these R05 rules may be pharmacologically and biologically active, in addition to Veber's rule of rotatable bonds (RB) of <10 and PSA of <140 A 2 (ref. 29) for orally bioactive compounds. Accordingly, all the compounds except curcumin glucuronide follow both Lipinski's and Veber's rule. Overall, 13 compounds possess oral bioavailability by obeying the drug-likeness rule tabulated in (Table 3).

Biological signicance of receptor-ligand docking
Molecular docking is a key tool to understand the mode of binding of a compound to the active site of the target proteins.  Therefore, the present study aimed to nd the binding interaction of curcuminoid derivatives to key target residues of SARS-CoV-2 through docking studies using the CDOCKER algorithm, a grid-based CHARMm simulation tool. Among the  14 compounds docked to the site of RBM, only two forms of curcumin, namely, its keto and enol forms, interact to form strong hydrogen bond interaction with the key residues Q493, N501, Y505, Y489, and Q498 with the CDOCKER docking score of À20.753 kcal mol À1 and À16.8067 kcal mol À1 (Fig. 8A and B), respectively. Thus, the receptor-ligand interaction study indicated that both forms of the curcumin pharmacophore have the ability to interact with the spike glycoprotein, anchoring the residues of ACE2. Also, blocking these residues possibly does not facilitate its interaction with the human ACE2 receptor, and thereby viral infection can be controlled.
A small envelope protein (E-protein) is another drug target protein that plays three functional roles in SARS-CoV-2, including viral assembly, 34 pathogenesis of the virus 35 and release of virions. 36 Targeting the E-protein functions aborts the viral assembly process, and consequently the formation of immature virions. 8 On contrary, the exact function of the envelope protein is still enigmatic because various studies have shown that even without E-protein, the virus utilizes accessory proteins to form its core structure, but the virus efficiency is reduced by a hundred-to a thousand-fold during morphogenesis. 37 During docking, it was observed that curcumin keto (À19.174 kcal mol À1 ) and enol (À14.115 kcal mol À1 ) interact with the TM alpha helix of the Eprotein to form hydrogen and hydrophobic interactions ( Fig. 9A and B), respectively. Thus, curcumin may be a candidate compound for treating SARS-CoV-2.

Receptor-ligand interaction pharmacophore generation for SAR analysis of curcumin derivatives
A pharmacophore is a key tool to derive a structure-activity relationship (SAR) to distinguish the active functional groups of chemical compounds that are responsible for their biological activity. Curcumin is already known for its pleiotropism of pharmacological activity. 38 Previously, the remarkable effect of dietary phenolic curcumin (1,3-dicarbonyl group) was reported for HIV drug targets, including integrase, Tat-mediated transactivation of the HIV-LTR and protease. 39 In the current study, SAR was probed for two forms of curcumin, i.e. its keto form (1,3-dicarbonyl group) and enol form (1,3-keto-enol). To identify the probable functional moieties of curcuminoids, the receptor-ligand interaction complex was used to enumerate the pharmacophore. The pharmacophore interaction features generated for keto and enol curcumin are depicted in Fig. 10. It can clearly be observed that the keto form of curcumin shares  Uniquely, the keto form of curcumin containing the ]O chemical group showed a tendency to bind with the active site amino acid in both drug targets. In contrast, in the enol form,  the ]O group is replaced with OH, which showed no binding interaction with the receptor sites. This reveals that the keto groups attached to the curcumin structure can be considered the important and crucial functional pharmacophore. Similarly, the enol form shares common pharmacophore hydrophobe aromatics and hydrophobes with the drug targets. Furthermore, to conrm this, the signicance of substructure elimination of the -OCH 3 moiety was studied. Accordingly, the removal of the OCH 3 moiety in the curcumins showed no docking or interaction pharmacophore generation, which implies that together with the main core of curcumin, the OCH 3 functional group is also required for binding interaction with the drug targets of SARS-CoV-2.

Binding energy and stability of the complex
Understanding the affinity of a compound to its target protein/ receptors is the main objective of the structure-based drug design and drug discovery process. Specically, the binding of a lead molecule to a receptor or signalling protein may alter its biological activity, and thus molecular structure recognition is considered a fundamental component in the virtual screening of drugs.
This can easily be achieved by means of calculating the binding energy of a docked receptor-ligand complex. The binding energy of the S-glycoprotein and envelope protein with the keto and enol forms of curcumin was computed using the equation (energybinding ¼ energycomplex À (energyligand À energyreceptor)), 40 where the negative energy of the binding complex shows the strength of the protein-ligand interaction ( Table 4). A molecular dynamics study was carried out to study the stability of the receptor-ligand interaction. At the end of the study, binding site conformational changes were observed for the SARS-CoV-2 drug target proteins. The S-glycoprotein and envelope protein bound with the keto form of curcumin showed less secondary structural conformational changes compared to that with the enol form ( Fig. 11 and 12), respectively. Further, signicant exibility and departure from the initial structure from molecular dynamics were estimated using the root mean square deviation (RMSD).
According to Fig. 13A, it can be observed that the RMSD deviation of the S-glycoprotein_keto deviations is within 2.1 A, also the stability of the complex was maintained throughout the molecular dynamics study. In contrast, the S-glycoprotein_enol complex gradually changed to maintain its structural stability, This journal is © The Royal Society of Chemistry 2020 RSC Adv., 2020, 10, 31385-31399 | 31395 but the deviation did not exceed 2.5 A. Subsequently, the conformational uctuation of the envelope protein with the keto and enol curcumin complexes showed a sharp peak at the beginning of the molecular dynamics. Conversely, aer 0.2 ns, the conformational changes of both curcumin complexes were less up to the end of the molecular dynamics study. Another time-dependent analysis is the radius of gyration (R g ), which is used to measure the compactness of proteins. Fig. 13B illustrates that the S-glycoprotein complexes with the keto and enol form showed negligible changes in protein structural folding, with R g values in the range of 17.6 to 18. In contrast, the R g value of the envelope protein with the enol form of curcumin changed over time, but at the end of the study, the structural compactness of the protein was retained, although the envelope protein complex with the keto form of curcumin showed less uctuation and was stable throughout the molecular dynamics study. Thus, overall, the complex stability is relatively higher for the keto-bound complex than the enol complex.

Free and interaction energy parameter analysis for receptor-ligand complex
The combination approach of molecular mechanics energies with Poisson-Boltzmann surface area (MM/PBSA) is used to estimate the free energy of the binding of small ligands to biological macromolecules. The binding free energies for the two keto-curcumin-S-protein RBD and enol-curcumin-S-protein RBD complexes were estimated using the MM-PBSA method. The calculated binding free energy was À4.067 kcal mol À1 for enol-curcumin. Dongling et al. 41 described that the favorable interaction of a complex can be measured in terms of electrostatic interaction. The analysis of the free energy components showed the that the binding free energy of À7.006 kcal mol À1 and more favorable electrostatic interaction (Fig. 14) for the keto-curcumin-S-protein RBD complex are the main reason for the higher affinity of keto-curcumin compared to enolcurcumin. On the other hand, the van der Waals interaction is more favorable for both the keto-curcumin and enolcurcumin conformations by À21.671 kcal mol À1 and À21.168 kcal mol À1 , respectively. Finally, the interaction energy of the complex was calculated using the CHARMm energy, which exhibited a 0.98-fold change difference for keto to enol curcumin. Comparatively, all the energy values are more favorable for both forms of curcumin with slight fold changes in energy values.

Conclusion
SARS-CoV-2 has resulted in a devastating pandemic with global concern; however, present therapies in virology fail to prevent its effects. Currently, there is exigency in identifying novel leads with anti-viral properties to impede viral pathogenesis in the host system. Thus, two important curcuminoids of turmeric, i.e., its curcumin keto and enol forms, were demonstrated to be complementary to bind with the S-glycoprotein and envelope protein of SARS-CoV-2. However, the keto form of curcumin is more favourable for both these drug targets considering its docking score, binding energy and molecular dynamics simulation. Thus, this study indicates that surface proteins are key drug target proteins of SARS-CoV-2, and probably curcumin blocks essential biologically active drug target residues, thereby attenuating the viral infection. Hence, this computational biology approach identies curcumin as a drug candidate for further investigation in treating SARS-CoV-2. However, this was an initial study to identify the active pharmacophore of the compound and its binding site structural complementary for two drug targets. Thus, in the future, we aim to perform largescale molecular dynamics, and in vitro and in vivo experiments to conrm the efficacy of curcumin against SARS-CoV-2.

Homology modelling
The coronavirus protein sequences reported in the Indian State of Kerala such as the envelope protein (QHS34548), membrane protein (QIA98586) and spike glycoprotein_RBD region (QIA98583) structures not reported in the PDB databases were determined using the template-based modelling technique. All the sequences were downloaded from the National Genomics Data Center (https://bigd.big.ac.cn/ncov/). Two approaches of a basic local alignment search tool for protein (BLAST_p) were performed, where one was a direct hit on a server (https:// blast.ncbi.nlm.nih.gov/Blast.cgi) and the other technique was through the BLAST_NCBI search protocol in Biovia DS2019. The template structures were selected by considering maximum identity and query coverage with less positive and E-values, where template hits with SARS-CoV/MERS-CoV species were given preference. Query sequences were aligned with the template structure followed by build homology modelling using MODELLER 9.17v9. 42

Modelled structure validation
The homology model protein was subjected to quality assessment in various structural validation servers and tools. In Discovery Studio, the best model structure was selected based on three parameters, including superimposing the best model structures with the PDB template structure to calculate the root mean square deviation using the align and superimpose protein protocol. Statistically, a lower DOPE score and high PDF total energy represent the best quality of model structures with stable conformations. Besides, external web servers such as ERRAT (https://servicesn.mbi.ucla.edu/ERRAT/), Procheck (https:// servicesn.mbi.ucla.edu/PROCHECK/), and Rampage (http:// mordred.bioc.cam.ac.uk/$rapper/rampage.php) were used to assess the structure quality. Some structures that showed more outliers in the Ramachandran plot were further optimized by the energy minimization technique provided with an RMS gradient of 0.1 kcal mol À1 A À2 to prevent bad steric contact of atoms. 43 Besides, a deciency in side chain and loop amino acids was processed using side and loop renement by CHARMm simulation and force eld. 44 The best-validated structures were used for further docking analysis.

Chemical structure preparation
The natural form of 14 curcumin derivative compounds were retrieved from the PubChem compound and substance database (https://pubchem.ncbi.nlm.nih.gov/). All compounds were prepared using the ligand protocol to generate various 3D conformations, isomers, and tautomers of the compounds to remove duplicity and to x bad valances.

ADMET, drug-likeness and toxicity predictions
The curcumin derivative was submitted to the ADMET and TOPKAT tools of small molecule protocol for the in silico pharmacokinetics and pharmacodynamic studies. The pK a study included parameters such as human intestinal absorption (HIA), aqueous solubility, blood-brain-barrier penetration (BBB), cytochrome CYP2D6 inhibition, plasma protein binding (PPB) and hepatotoxicity. In the pharmacodynamic study, we included animal models such as NTP rodent carcinogenicity, rat_oral LD50, rat_TD50, fathead minnow LC50, Daphnia EC50, rat inhalational LC50 and Ames mutagenicity. All these models were developed and validated based on a quantitative-structure toxicity relationship (QSTR).

Molecular docking and binding energy
Molecular docking is a lock and key process to identify compounds with complementary structures to drug target proteins. The CDOCKER algorithm is a grid-based and molecular dynamic simulation-implemented docking protocol, which was utilized for docking compounds to the binding site of the protein. The spike glycoprotein and envelope proteins are two novel drug targets for the dreadful covid-19 virus, which are the binding sites identied using a receptor-based cavity tool in Discovery Studio 2019. However, the spike glycoprotein (Sglycoprotein) RBD region binding to ACE2 receptor interactions were analyzed using the analyze protein interface tool in macromolecules to spot hydrogen, hydrophobic and other interactions, which was considered as the binding site for compound docking (Fig. 15). Both proteins were placed in an equal grid spacing in right angles in the 3D direction of the input site sphere coordinates with the radius of À38.36X, 31.32Y, 3.33Z, 20.1 A and 37.37X, 33.46Y, 17.81Z, 26.873 A for the RBD region of the S-glycoprotein and envelope protein, respectively. The docked complexes were further analyzed and validated via the negative CDOCKER energy and receptorligand interactions. Further, to understand and quantify the strength of the interactions between a ligand and protein, the binding energies between all the docked poses of the ligand with receptor were studied by utilizing the binding energy protocol with CHARMm. This is a crucial process in drug discovery and lead optimization. The best pose with the lowest binding energy was considered for the molecular dynamics simulation.

Molecular dynamics and simulation
The top compounds with the best pose in receptor-ligand interaction for the two different drug target proteins aer docking validation were selected for molecular dynamics simulation. The whole system was subjected to CHARMm force eld to satisfy bonded and non-bonded interactions, followed by standard dynamics cascade in a ve-step simulation protocol. Initially, two 500 steps of minimization were performed using the steepest descent and conjugate gradient. 45 The minimized complexes were gradually driven from 50 K to a nal target temperature of 300 K, followed by equilibration simulations. Finally, production was performed for 1000 ps for four complexes. The leapfrog dynamics integrator and shake constraint were introduced throughout the molecular dynamics simulation to study the bonded and non-bonded interaction. The molecular dynamics trajectory was determined for complex stability and time-dependent analysis (RMSD and radius of gyration (R g )) using the Biovia Discovery Studio 2019 analyze trajectory protocol (Dassault Systèmes). 46

CHARMm interaction energy and molecular mechanics-Poisson-Boltzmann surface area (MM-PBSA) analysis
This study was focused on the impact of curcumoinds on the stability of the S-protein. Thus, the complex structure of ketocurcumin-S-protein RBD and enol-curcumin-S-protein RBD was used as a starting point for calculating the binding free energies. The 1000 ps molecular dynamics simulation was carried out using CHARMm. In CHARMm, molecular dynamics simulations are performed using a classical mechanics approach, in which Newton's equations of motion are integrated for all atoms in the system. 47 Thus, for each MDsimulated complex, we calculated the free binding energy (DG binding ) values for the 500 conformations of the MD trajectory. During the simulation, one conformation or snapshot was saved every 2 ps up to 1000 ps, and the nal DG binding is the average of 500 conformations of the receptor-ligand complex. In the MM-PBSA method, the free energy of the protein-ligand binding (DG binding ) is obtained from the difference between the free energies of the protein-ligand complex (G complex ) and the unbound receptor/protein (G protein ) and ligand (G ligand ) as follows: The output conformations of the molecular dynamics simulation were further sampled to study the interaction energy between sets of atoms across all conformations using CHARMm. The interface interacting residues of the receptor and ligand were selected as two atom sets together with a cut-off distance of 12-10 A to calculate the non-bonded interactions as follows: where, DE Int ¼interaction energy of ligand-protein, DE vdW ¼ van der Waals energy of ligand-protein, and DE elec ¼ electrostatic energy of ligand-protein.

Funding
This research project has not received any funding from any agencies.

Conflicts of interest
The authors conrm that this article content has no conicts of interest.