key: cord-1005022-tr0q2k1x
authors: Muhseen, Ziyad Tariq; Hameed, Alaa R.; Al-Hasani, Halah M.H.; Qamar, Muhammad Tahirul; Li, Guanglin
title: Promising terpenes as SARS-CoV-2 spike receptor-binding domain (RBD) attachment inhibitors to the human ACE2 receptor: Integrated computational approach
date: 2020-10-07
journal: J Mol Liq
DOI: 10.1016/j.molliq.2020.114493
sha: c2bf474cb4a6624d2619343445bef30dcdb3bbe6
doc_id: 1005022
cord_uid: tr0q2k1x

The spike protein receptor binding domain (S-RBD) is a necessary corona-viral protein for binding and entry of coronaviruses (COVs) into the host cells. Hence, it has emerged as an attractive antiviral drug target. Therefore, present study was aimed to target severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) S-RBD with novel bioactive compounds to retrieve potential candidates that could serve as anti-coronavirus disease 2019 (COVID-19) drugs. In this paper, computational approaches were employed, especially the structure-based virtual screening followed by molecular dynamics (MD) simulation as well as binding energy analysis for the computational identification of specific terpenes from the medicinal plants, which can block SARS-CoV-2 S-RBD binding to Human angiotensin-converting enzyme 2 (H-ACE2) and can act as potent anti-COVID-19 drugs after further advancements. The screening of focused terpenes inhibitors database composed of ~1000 compounds with reported therapeutic potential resulted in the identification of three candidate compounds, NPACT01552, NPACT01557 and NPACT00631. These three compounds established conserved interactions, which were further explored through all-atom MD simulations, free energy calculations, and a residual energy contribution estimated by MM-PB(GB)SA method. All these compounds showed stable conformation and interacted well with the hot-spot residues of SARS-CoV-2 S-RBD. Conclusively, the reported SARS-CoV-2 S-RBD specific terpenes could serve as seeds for developing potent anti-COVID-19 drugs. Importantly, the experimentally tested glycyrrhizin (NPACT00631) against SARS-CoV could be used further in the fast-track drug development process to help curb COVID-19.

Human civilization has faced many pathogenic outbreaks, in the past few decades.

Among these, the occurrence and reoccurrence of CoVs related outbreaks have caused fetal respiratory disorders [1, 2] . The positive-sense signal strand RNA CoVs have envelop/capsid around and poses high frequency of genomic mutation and recombination, hence, categorized as extremely evolving viruses [3, 4] . Around six species of human coronavirus are known to exist (HCoV-OC34, HCoV-229E, MERS-CoV, HCoV-NL63, SARS-CoV, and HCoV-HKU1) that are responsible for many respiratory disorders [3, 5, 6] . At present, the world is fronting a severe respiratory coronavirus disease 2019 pandemic mainly due to a newly evolved infectious CoV strain: SARS-CoV-2 [4, [7] [8] [9] . Up till now (September 25, 2020), COVID-19 has affected ~32.1M people worldwide and caused ~980 K deaths. This havoc has not stopped yet and still thousands of people are being affected on daily basis. Presently, none of the vaccine/drug is available in the markets for human use against COVID-19 [10] [11] [12] , while many of the drugs are in the trial process and a few successful trials have been reported as well, but all are on the very initial stage [13] . Among these drugs, remdesivir has emerged as a promising drug [14] that is now being trailed/tested on patients suffering from severe disease [14, 15] .

The SARS-CoV-2 RNA enters into the human cell by a dense glycosylated spike (S) protein [16] that is basically a fusion protein belonging from the trimeric class I. The S protein normally occurs in a metastable prefusion structure that can be modified substantially for the fusion of viral membrane into the host cell [16, 17] . The infection of COVID-19 is mainly directed through the interaction of the host cell receptor, CoV J o u r n a l P r e -p r o o f Journal Pre-proof enveloped affixed S protein and angiotensin converting 2 (ACE2) enzyme [18, 19] . The binding S-protein has further two subunits, S1: [the receptor-binding domain (RBD)] and S2 (carryout the union of viral and host cell membrane) [18, 19] . The subunit S2 of SARS-CoV-2 is conserved and ~99% similar to the bat SARS-CoV and human SARS-CoV. In contrast, the S1 subunit is evolved differently (only ~70% similarity with human and bat SARS-CoVs) [19, 20] . The affinity among S-RBD of SARS-CoV-2 and H-ACE2 is ten times higher than SARS-CoV S-RBD (year 2003), which infers that H-ACE2 is an explicit receptor accountable for the attachment of a virus to the host cell [20, 21] . The inter-species transmission of CoVs can be attributed to the SARS-CoV S-RBD residues (Leu472, Tyr442, Asn479, Thr487, and Asp480) that undergo natural selection in SARS-CoV-2 [21, 22] . The residues such as, Lys353 and Lys31, present on H-ACE2 are said to be liable for the virus-binding hotspot residues that ultimately bind S-protein [22, 23] .

The Lys31 of H-ACE2 receptor is composed of a salt bridge between Glu35 and Lys31, while hotspot 353 contains an alternative salt bridge between Asp38 and Lys353, bounded by hydrophobic surroundings [20, 23] . H-ACE2 is recognized in SARS-CoV-2 by its residue (Leu455 and Gln493) that are liable for desirable molecular interfaces with hotspot 31, thus improving the binding of a virus to H-ACE2. Furthermore, other main residues of S-protein offer more backing for Lys31 (Leu455, Ser494, Phe486,). In serine based residue (494) of SARS-CoV-2, strengthens the structural solidity of Lys353 [20, 24] . Due to its crucial function the S protein, it is important for designing a drug. The Phytochemicals act as secondary metabolites and active constituents in medicinal plants e.g., flavonoids, alkaloids, and terpenes peptides [25] [26] [27] . These biologically active molecules are parts of defense system in plants hence, in turn, reduce the infections [25, 28, 29] . As a part of the remedial function, these compounds forage, scavenge and obstruct the viral entry and DNA\RNA replication [28, 29] . Among these, terpenes are most common and represent a large category of secondary metabolites that are composed of 5-C isoprene units that are interlinked in a complex manner [30] . Most of the terpenoids are biologically active and are widely used worldwide in order to treat many ailments such as, derivatives of Taxol used to treat cancer artemisinin has antimalarial effects [31] .Moreover, for the FDA-approved antiviral drugs, numerous in-silico researches have been conducted to differentiate new phytochemical compounds against SARS-CoV-2/COVID-19 [4, [32] [33] [34] [35] [36] [37] [38] .

There is a need to conduct a comprehensive analysis of the interactions of H-ACE2 receptor and S-RBD of SARS-CoV-2 in order to establish the remedies against COVID-19. Therefore, present study was conducted using high throughput screening (HTS) coupled with various computational techniques, including structure-based molecular docking approach, MD simulation approach and free binding energy MM/PBSA analysis, for the identification of novel medicinal plant based terpenes that can be helpful against SARS-CoV-2 infection. The comprehensive procedure and flowchart of the current study is given as Figure 1 . The outcomes of this work are significant in the field of anti-COVID-19 terpenes drug efficacy and give a push to its progress as a drug against SARS-CoV-2.

For the development of Natural-based novel drug, a dataset of 1000 plant bioactive terpenes compounds with reported therapeutic potentials were downloaded from NPACT database [39] and MPD3 database [26] . All the compounds were optimized prior to docking using Autodock tools [40] . Polar hydrogens were added and gassteiger charges were computed of each compound in the library. Later, the compounds were saved in pdbqt format for docking purpose.

Molecular docking of the terpenes dataset was conducted by Autodock Vina [41] . The of SARS-CoV-2 S-RBD crystal structure bound with the ACE2 receptor protein [42] was attained from the PDB database (PDB ID 6LZG) for docking purpose. All of the crystalline compounds were detached from the RBD structure prior to docking. The RBD structure was prepared using AutoDock tools [40, 43] . Gassteiger charges were optimized and Polar hydrogens were added. The structure was saved in pdbqt format for docking.

The Grid box was created with the size of 45 Å x 45 Å x 45 Å and centered on the important residues of RBD domain which binds to the ACE2 receptor. 20 poses for all compounds were created and results were analyzed on the basis of docking scores (kcal/mol). were chosen and visually inspected using PyMol [44] .

The available bioinformatics tool SwissADME (http://www.swissadme.ch/index.php) [45] was used for finding drug-likeness attributes. Lipinski's rule of five [46] was used to analyze the properties such as; hydrogen bond donor (HBD), hydrogen bond acceptor (HBA), molecular weight (MW), and lipophilicity (log P). PreADMET (https://preadmet.bmdrc.kr/) and pkCSM (http://biosig.unimelb.edu.au/pkcsm/) servers were employed to determine the ADMET (metabolism, distribution, excretion, absorption, and toxicity) attributes of target molecules.

Molecular dynamics simulations based on categorical solvent were conducted using GROMACS 5.1.4 (GROningen Machine for Chemical Simulations) to measure the binding and stability attributes of complex [47] . CHARMM36 force field was applied to achieve topology attributes of protein [48] . For force field ligand and topology parameters, CHARMM General Force Field (CGenFF) was employed [49] . Ligandprotein complexes resulted after docking were solvated in the cubical box through TIP3 water model. The neutralization of the system was performed by the addition of a suitable number of CL-ions. Energy was minimized in all systems through 10000 steps of steepest descent approach. In order to balance/maintain the temperature of the system ranged from 0 to 300 K slowly for 250 ps, NVT ensemble was used. Moreover, simulation of the J o u r n a l P r e -p r o o f Journal Pre-proof system was performed under NPT ensemble (at 1.0 bar and 300 K) for 250 ps. In order to quantify the long-range interfaces, Particle Mesh Ewald was used while for short-range Electrostatics and Van der Waals, a distance (10A) cut off was used [50] . For restraining all bonds, the use of LINCS (Linear Constraints Solver) algorithm was done [51] . Finally, 50 ns making run was performed and trajectories were saved after each 2 fs for all complex separately. Root Mean square fluctuations (RMSF), root mean square deviation (RMSD), SASA analysis and Radius of gyration (Rg) were done using GROMCAS modules.

PCA/essential dynamics statistical tool was used to identify similarities and dissimilarities in a particular dataset and to detect patterns in a particular dataset [52] .

The main idea lies in reducing the conformational space into essential subspace with varying degrees of freedom that comprises most of the positional atomic fluctuations.

Thus, atomic fluctuations are expressed in the form of a covariance matrix and the transformation of the conformational space is done by diagonalization of the matrix. Only Cα atoms mean interatomic distance variation was applied to generate covariance matrix.

Gmx covar module was used to construct gmx anaeig and covariance matrix module was employed to project trajectories on eigenvectors.

MM/PBSA was used to determine the binding free energies of protein-ligands [53] [54] [55] . G_mmpbsa module of GROMACS v.5.1.4 was employed to approximate the binding energies of simulated complex compounds [56] . The binding free energy 

In order to obtain understanding about the interaction between H-ACE2 and SARS-CoV-2 S-RBD, first we analyzed the interface of S-RBD and H-ACE2 complex (Figure 2A) .

We found as mentioned in literature, the overall interface is arbitrated essentially by polar interactions. A lengthy loop section of the RBD spans the arch-shaped α1 helix of the ACE2 similar to a bridge. The α2 helix and a ring that links the β3 and β4 anti-parallel J o u r n a l P r e -p r o o f strands, stated as loop 3-4, of the ACE2 also make partial assistances to the synchronization of the RBD. Overall, 14 potent residues of S-RBD (Phe486, Asn487, Ala475, Glu484, Tyr489, Leu455, Gln493, Tyr453, Ser494, Tyr449, Gln498, Gly496, Gly446 and Asn501) were identified that were making significant binding with ACE2 and contributing the stability of complex ( Figure 2B ).

In the present study, a library of terpenes was screened against the SARS-CoV-2 S-RBD domain. The approach is on the basis of computationally fitting molecules for a target protein using 3D structure of the active site, subsequently the ranking of the docked compounds and interaction analysis. The residues located at the interface section of SARS-CoV-2 RBD-ACE2 were targeted in order to screen the terpene compounds using AutoDock Vina. The two step docking experiments were performed. In the first step of docking, the RBD was retained rigid and 1000 terpene compounds were docked and filtered on the basis of the docking score. In the second step of docking, top 100 compounds were selected and docked flexibly to the targeted sites. Only important interface residues of the RBD were treated flexibly in the flexible docking experiment.

After flexible docking, 50 top hits were visually inspected on the basis of the interactions.

Glycyrrhizin (NPACT00631) hits having energies -11 kcal/mol, -10.3 kcal/mol and -9.5 J o u r n a l P r e -p r o o f kcal/mol, respectively, were selected for more analysis. The structures and characteristics of three finally screened compounds are displayed in Table 1 .

Discovery studio software was used to prepare two-dimensional plot of the molecular interaction network of ligands with S-RBD while the docked poses for each of these molecules are portrayed in Figure 3 . Out of 1000 terpenes, NPACT01552, NPACT01557

and NPACT00631 showed highest binding energies (-11 kcal/mol, -10.3 kcal/mol and -9.5 kcal/mol, respectively) and interact with RBD remains through five, five and six hydrogen bonds, respectively. As mentioned above, RBD domain binds with the ACE2 receptor protein through these important interface residues Gly502, Gly496, Thr500, Table 2 summarizes the energy scores and interacting residues which interact with compounds of selected terpenes.

For further validation of the drug likeliness capability of selected compounds, all three terpene molecules were subjected to the Lipinski's rule of five: absorption, distribution, metabolism, excretion, and toxicity testing. These parameters were evaluated using different thresholds. All three hit compounds passed the threshold of drug ability ( Table   1 & 3).

Molecular dynamics (MD) simulation is the most widely used technique to assess the micro interaction between ligand structure and protein. Simulation of all complexes was done by 50 ns MD in order to access further its stability and dynamics. After MD simulations, RMSD, RMSF, Rg, and SASA were assessed.

Root mean square deviations of backbone atoms of all the systems were analyzed to observe the stability of systems. The RMSD calculation of a structure along its trajectory relative to the reference structure tells the conformational changes in the structure. The RMSD curve indicates that the RBD-Apo, NPACT01552, and NPACT01557 complexes remained stable throughout the simulation with the average RMSD 0.12, 0.11 and 0.13 nm, respectively. There were no significant differences observed in NPACT01552, and NPACT01557 complexes. As a whole, the results of RMSD signify that the binding of J o u r n a l P r e -p r o o f these two compounds with RBD of S-protein is steady. However, compound NPACT00631 was stable in the start, but after 20 ns this compound fluctuated up to 0.24 nm RMSD until the end of simulation, which may indicate that this compound is more flexible as can be seen from Figure 4A Complex NPACT00631 has a larger conformational flexibility as compared to the RBD-Apo, complex NPACT01552, and NPACT01557 as evident from Figure 4B .

The radius of gyration (Rg) depicts the system density, and eventually impacts the folding rate as well as proteins stability. Rg was used to assess the complexes compactness and it was revealed that Rg of all systems was steady with the RMSD of J o u r n a l P r e -p r o o f system. This depicts that protein was stable and compressed throughout the 50 ns time Figure 4C .

Solvent accessible surface area (SASA) is alternative method in order to retain the protein folding and stability. The calculated SASA values for wild type and mutants are displayed in Figure 4D . The average SASA values for RBD-Apo, NPACT01552, NPACT01557 and NPACT00631 were found to be 109, 108, 108 and 109 nm 2 , respectively signifying that there were no notable differences in the available area of all the systems throughout the simulation process.

All the simulation systems covered the same conformational space except complex NPACT00631. The trace of the covariance matrix and eigenvalues for RBD-Apo, NPACT01552, NPACT01557 and NPACT00631 was 1.47, 1.33, 1.47 and 2.84 nm 2 respectively. The eigenvalues for RBD-Apo, NPACT01552 and NPACT01557 were nearly same, which indicate that these systems did not undergo random fluctuations, which can also be observed from the 2D projection of trajectories in Figure 5 . In the 

To clarify the energetics of the binding of compounds NPACT01552, NPACT01557 and NPACT00631 with the S-RBD quantitatively, we calculated the MM/PBSA calculations. 

The energy contributions of the active residues to the total binding energies were also measured for each terpene-SRBD complex to comprehend the residues that contributed significantly to the interactions of inhibitor-protein. The values for binding energy of the residues at the interface that significantly contributed to the total binding energy are displayed in Figure 6 . Residues Leu455, Phe456, Ala475, Phe486 and Tyr 489 contributed significantly in the total binding energy.

The comprehensive knowledge of the viral receptor identification mechanism is a key to understand the pathogenic, infectious, and the host specific nature of COVID-19 that is further important for developing remedial therapies, drugs and antiviral cures [20] . Till now, none of the drug or therapy can be specifically prescribed against SARS-CoV-2, while the application/development of new drugs is expensive and requires a lots of trial J o u r n a l P r e -p r o o f steps, hence, time consuming [4, 65, 66] . More worse is the prediction of WHO that COVID-19 may become endemic, which is an alarm in the scientific world that it has become more crucial to develop drugs against it immediately [67] . In this regard, the computational methods are scientifically sound to design extremely specific inhibitor for viral proteins, as an antiviral drug development/discovery. Recently, structural modeling has predicted a few drugs against SARS-CoV-2 proteins along with many other antiviral drugs [2, 65, 66, 68] .

The life cycle of the virus is consisted of four main phases, in the first phase, virus is attached to the host cell by identifying the specific receptors. In the second phase, it is fused with the host cell membrane and viral genetic material (RNA) is transferred into the host cell, (in next stage) later on, the viral genetic material gets a hold of host genome and start making structural and non-structural proteins and finally the viral parts are assembled and phages are ejaculated from host cell that are ready to attack other host cells [20, 69] . In order to cure the infections caused by RNA viruses such as dengue virus, chikungunya virus, SARS, Ebola virus, Sindbis virus MERS, etc., it is necessary to target these viruses at any of the stage of their life cycle [20, 70] . The drugs that block the entry of the virus into the host cell or inhibit the viral replication have already been introduced in the market for chikungunya, dengue, and other RNA viruses [20, 71] .

Inhibition of the capsid synthesis can also be a key for preventing viral infection at a budding stage [20, 72] . Similarly, such drugs can be developed against SARS-CoV-2 for targeting/inhibition the viral growth at any stage of its life cycle. The envelop of SARS-CoV-2 contains S-protein that is liable for the interaction with ACE2 receptor on the host cells through RBD unit [19] . This is the basic and first step for the interaction of the J o u r n a l P r e -p r o o f SARS-CoV-2 with host cell and creation of infection. The antiviral drug that can target S-RBD-ACE2 interaction that is able to render the entry of virus into the host cell and hence, offer a fast resolution to combat SARS-CoV-2 inflammation [20, 73] .

Till now, there is no precise antiviral specific strategy for the treatment of COVID-19 (vaccine/drug), supportive care such as, supplementation with the combination of broad-spectrum antivirals, antibiotics, corticosteroids as well as convalescent plasma is widely employed [74] . More than 3/4 th quarter of the world's population is using herbal medicines [26, 28] . Out of the total 2,50,000 species of flora present on the earth ~80,000 poses medicinal attributes [25, 28] . Many bioactive compounds are the part of plant biomass such as alkaloids, peptides, terpenes and flavonoids [4, 29] . In all of these compounds, terpenes are the most abundant and largest group of 5-C secondary compounds/metabolites that consisted of isoprene subunits that can be interlinked with each other in 1000 s of ways [75] . The in silico research has proposed the screening of many novel phytochemical that can be used as inhibiting agent for many of key protease of SARS-CoV-2 [4, [32] [33] [34] [35] [36] [37] [38] .

In accordance with the results of 1000 focused terpenes with reported therapeutic potentials library screening, molecular docking was executed to obtain binding mode insights and crucial molecular interfaces of the designated terpenes and S-RBD protein obtained from SARS-CoV-2. Likewise, compounds were evaluated through MD approach, MM-PBSA and per-residue energy disintegration revealed that the relative significance of amino acid responsible for the binding reinforced under the influence of various molecules having binding free energy. NPACT01552, NPACT01557 and NPACT00631 were projected to be possible inhibitors against S-RBD of SARS-CoV-2 J o u r n a l P r e -p r o o f as these compounds raise the binding energies and as a result can interact and block the residues of S-RBD that recognize the attachment sites on the host cell. Furthermore, these complexes also exhibit successful ADMET attributes and these compounds are also sideeffect free/have very low risk. The other drugs having poor ADMET can produce toxic metabolites [50] or may not be able to cross cell membranes barriers [51] , and the potential solubilities were also approximated within the suitable range [52] . Many of the compounds were proposed as non-inhibitors for the isozymes, such as, cytochrome P450, but these drugs play a very crucial role in metabolism and enzymes inhibition that may leads to toxicity and deprived defecation in liver [53] . (Glycyrrhizin) has the energy score of −9.5 kcal/mol and is at number three.

NPACT01557 is extracted from Aralia Dasyphylla and well documented to have anticancer effects as well [77] . NPACT00631 is also known as Glycyrrhizin that is present in Glycyrrhiza glabra and has been successfully used to treat the HIV-1 and chronic hepatitis C in past [78] [79] [80] [81] [82] . Before this, Cinatl et al, [79] has determined the J o u r n a l P r e -p r o o f Tables   Table 1 Details of the terpenes selected as top hits. J o u r n a l P r e -p r o o f 

Severe Acute Respiratory Syndrome Corona virus 2; COVID-19: Coronavirus disease 2019; S: Spike; RBS: Receptor Binding Domain

Funding This work was supported by the grants from the National Natural Science Foundation of

the Program for New Century Excellent Talents in University (NCET-12-0896), and the Fundamental Research Funds for the Central Universities (GK201403004)

Travel medicine and infectious disease

Data curation, Validation, Writing-Reviewing and Editing. Muhammad Tahir ul Qamar: Software, Validation, Writing-Reviewing and Editing. Guanglin Li: Conceptualization, Methodology, Supervision

Authors would like to acknowledge Shaanxi Normal University, Xi'an, China for providing facilities for this study. 

J o u r n a l P r e -p r o o f

☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:J o u r n a l P r e -p r o o f