Next Article in Journal
N-Acetylaspartyl-Glutamate Metabolism in the Cingulated Cortices as a Biomarker of the Etiology in ASD: A 1H-MRS Model
Next Article in Special Issue
Combining Different Docking Engines and Consensus Strategies to Design and Validate Optimized Virtual Screening Protocols for the SARS-CoV-2 3CL Protease
Previous Article in Journal
PROTACs and Building Blocks: The 2D Chemical Space in Very Early Drug Discovery
Previous Article in Special Issue
Structural Impacts of Drug-Resistance Mutations Appearing in HIV-2 Protease
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Computational Determination of Potential Multiprotein Targeting Natural Compounds for Rational Drug Design Against SARS-COV-2

1
Key Laboratory of Ministry of Education for Medicinal Plant Resource and Natural Pharmaceutical Chemistry, Shaanxi Normal University, Xi’an 710062, China
2
School of Life Sciences, Shaanxi Normal University, Xi’an 710062, China
3
Department of Medical Laboratory Techniques, School of Life Sciences, Dijlah University College, Baghdad 00964, Iraq
4
Department of Biotechnology, College of Science, University of Diyala, Baqubah 32001, Iraq
5
Foundation University Medical College, Foundation University Islamabad, Islamabad 44000, Pakistan
*
Author to whom correspondence should be addressed.
Molecules 2021, 26(3), 674; https://doi.org/10.3390/molecules26030674
Submission received: 6 December 2020 / Revised: 18 January 2021 / Accepted: 21 January 2021 / Published: 28 January 2021

Abstract

:
SARS-CoV-2 caused the current COVID-19 pandemic and there is an urgent need to explore effective therapeutics that can inhibit enzymes that are imperative in virus reproduction. To this end, we computationally investigated the MPD3 phytochemical database along with the pool of reported natural antiviral compounds with potential to be used as anti-SARS-CoV-2. The docking results demonstrated glycyrrhizin followed by azadirachtanin, mycophenolic acid, kushenol-w and 6-azauridine, as potential candidates. Glycyrrhizin depicted very stable binding mode to the active pocket of the Mpro (binding energy, −8.7 kcal/mol), PLpro (binding energy, −7.9 kcal/mol), and Nucleocapsid (binding energy, −7.9 kcal/mol) enzymes. This compound showed binding with several key residues that are critical to natural substrate binding and functionality to all the receptors. To test docking prediction, the compound with each receptor was subjected to molecular dynamics simulation to characterize the molecule stability and decipher its possible mechanism of binding. Each complex concludes that the receptor dynamics are stable (Mpro (mean RMSD, 0.93 Å), PLpro (mean RMSD, 0.96 Å), and Nucleocapsid (mean RMSD, 3.48 Å)). Moreover, binding free energy analyses such as MMGB/PBSA and WaterSwap were run over selected trajectory snapshots to affirm intermolecular affinity in the complexes. Glycyrrhizin was rescored to form strong affinity complexes with the virus enzymes: Mpro (MMGBSA, −24.42 kcal/mol and MMPBSA, −10.80 kcal/mol), PLpro (MMGBSA, −48.69 kcal/mol and MMPBSA, −38.17 kcal/mol) and Nucleocapsid (MMGBSA, −30.05 kcal/mol and MMPBSA, −25.95 kcal/mol), were dominated mainly by vigorous van der Waals energy. Further affirmation was achieved by WaterSwap absolute binding free energy that concluded all the complexes in good equilibrium and stability (Mpro (mean, −22.44 kcal/mol), PLpro (mean, −25.46 kcal/mol), and Nucleocapsid (mean, −23.30 kcal/mol)). These promising findings substantially advance our understanding of how natural compounds could be shaped to counter SARS-CoV-2 infection.

1. Introduction

Coronaviruses (CoVs) cause infection of the upper respiratory tract in higher mammals and humans [1], and several outbreaks have been associated in the recent past with CoVs reported first time in the year 2002 as SARS, in 2012 as MERS, and in late 2019 as COVID-19 [2,3,4,5]. The recent pandemic of COVID-19 is caused by a relatively new strain named SARS-CoV-2 [6,7,8]. The virus origin is thought to be zoonotic, with potential of transmissibility between person-to-person, resulting in an exponential rise in the number of confirmed cases worldwide [9,10]. Through December 2020, more than 220 countries reported the virus, with more than 64 million individuals infected, and thousands are still getting infected each day. Approximately, the virus has a mortality rate between 5% to 10% [11,12]. Additionally, due to mandatory lockdowns, isolation, and quarantines, millions of lives have been disturbed. The pandemic also badly affected global health, society, and the economy, and these sectors are facing significant challenges [13]. Three vaccines (by Pfizer, Moderna, and AstraZeneca) are authorized by WHO for emergency use and are available to very limited populations. No specific anti-SARS-CoV-2 drugs are currently recommended for SARS-CoV-2 treatment, making the situation difficult to handle. Supportive therapeutics and preventative measures are being taken and are productive in managing the virus [14,15]. Various efforts to target critical proteins of SARS-CoV-2 pathogenesis, including Spike receptor-binding domain (RBD) [16,17,18], main protease (Mpro) [19], Nucleocapsid N terminal domain (NTD) [20], RNA-dependent RNA polymerase (RdRp) [21], papainlike protease (PLpro) [22], 2′-O-RiboseMethyltransferase [23], viral ion channel (E protein) [24], and angiotensin-converting-enzyme 2 receptor (ACE2) [25], are on the way. Targeting multiple pathogenesis specific proteins within a close network of interaction or dependent functionality would effectively propose effective drugs against the SARS-CoV-2 [26].
SARS-COV-2 Spike protein is key to the host cell infection pathway as it mediates ACE2 recognition, attachment, and fusion to the host cell [16]. The RBD of S1 subunit of the Spike trimer binds explicitly to the ACE2 receptor [27]. This RBD region is an attractive target for therapeutics as it contains conserved residues that are essential in binding to ACE2 [27]. The Mpro of coronaviruses has been studied thoroughly for drug making purposes. These are papainlike proteases involved in processing replicase enzymes [28]. It has 11 cleavage sites in 790 kD-long replicase lab polypeptide, demonstrating its prominent role in proteolytic processing [19,29]. High structural similarity and sequence identity are seen in Mpro from SARS-CoV-2 to that of the SARS-CoV Mpro. It comprises two catalytic domains: chymotrypsin and picornavirus 3C protease like domain. Each contains β-barrel that are six in number and are antiparallelly containing active diad H41 and C145 [30]. These proteases have emerged as essential drug targets as they have a crucial role in replication. Furthermore, inhibitors of Mpro are found to be significantly less cytotoxic as the protein share less similarity with human proteases [31]. Preliminary studies have suggested that HIV protease inhibitors, lopinavir/ritonavir, could be potentially used against SARS-CoV-2 [32]. Additionally, HIV protease inhibitor, Darunavir, and HCV protease inhibitor, Danoprevir, are under clinical studies and in vivo trials for the treatment of SARS-CoV-2 infection [33]. The PLpro enzyme is vital in processing the polypeptide to produce a functional replicase complex and aids in viral spreading [22]. PLpro also plays a role in evading host antiviral immune responses by cleaving proteinaceous modification on the host protein after the post-translation phase [34]. Thus, targeting this enzyme is useful in highlighting therapeutic strategies that can suppress the virus infection and prompt antiviral immunity. The N protein is significant in viral RNA replication and its packing into new virions, making this protein a good candidate for newer drug identification that is specific and biological active [20].
In silico screening of drugs using different computer-aided drug designing applications greatly accelerate the rational drug design process. Ultimately, this saves time, and extra cost goes into the experimentation of leads that fail in the drug discovery process [35,36,37,38,39,40]. In this investigation, we performed a blind docking approach, followed by molecular dynamics (MD) simulation coupled with binding free energy techniques that dissect the structural dynamics and energy basis of molecular recognition [41,42]. The MPD3 phytochemical database [43] along with a pool of natural antiviral compounds were used against multiple SARS-CoV-2 protein targets to understand their binding mechanism and put forward a hypothesis on how to further optimize these structures to enhance selectivity and maximize anti-SARS-CoV-2 biological potency [44,45,46,47]. A schematic summary of the methodology used in this work is provided in Figure 1. The study results might have potential applications in designing new leads against SARS-CoV-2, which can target its multiple proteins as depicted in this study.

2. Results and Discussion

2.1. Molecular Docking

Molecular docking is a modeling approach investigating how the receptors and ligands fit together and how the enzymes interact with the ligands [48,49,50,51,52]. Docking calculations were performed in triplicate, and the compound conformations were ranked according to the binding energy in kcal/mol. We used remdesivir as control in docking. The compounds ranked consistently on top with the each receptor and showed a stronger binding score compared to remdesivir were selected for the downward analysis. A general overview of the binding energy of the compounds against the receptors used is presented in Figure 2. The top compound complex with each receptor was generated and subjected first to visual inspections to decipher atomic level interaction and determine the binding conformation. The docking analysis demonstrated glycyrrhizin followed by azadirachtanin, mycophenolic acid, kushenol-w, and 6-azauridine as the best binders among the ~5000 compounds used in this study. The 2D structures of these compounds are presented in Figure 3. Glycyrrhizin also showed stable interactions with the hotspot residues of SARS-CoV-2 spike protein receptor binding domain (RBD) in our previous study [53]. Glycyrrhizin-docked complex of each SARS-CoV-2 protein can be explained separately.

2.1.1. Mpro–Glycyrrhizin Complex

The Mpro of SARS-CoV-2 is a crucial enzyme and attractive drug target because of its central role in virus transcription and replication [54]. The docking study reported glycyrrhizin again as the best binder among the compounds used to the substrate-binding site of the Mpro (Figure 4). As seen in the binding with other receptors, the compound (2S,3S,4S,5R,6R)-6-(((2S,3R,4S,5S,6S)-6-carboxy-2,4,5-trihydroxytetrahydro-2H-pyran-3-yl)oxy)-3,4,5-trihydroxytetrahydro-2H-pyran-2-carboxylic acid was revealed to contribute in significant hydrogen bonding and other weak interaction at the active pocket of Mpro. At the binding cavity, the compound engages Asn238 through multiple hydrogen interactions, as well as Asp289. The rest of the compound structure makes a network of hydrophobic interactions mainly dominated by van der Waals contacts. To elucidate further the binding specificity and affinity of the glycyrrhizin for the active pocket residues of Mpro, the interaction profile was compared and contrasted with that for the reported cocrystallized N3 inhibitor [55]. Very low similarity in the binding interaction profile between the compounds was noticed; however, because of the difference in the compound structure, size, and preferred binding site, the pocket residues in contact with glycyrrhizin are close to the N3. This difference in the binding interaction points to the different glycyrrhizin-binding mechanism, where the active moiety favors binding with the P5 binding pocket that is absent in the case of the Mpro–N3 complex. The residues, particularly Asp197 and Thr198, flanked the active site, and any molecule involved in binding with these residues interfere with the natural substrate-binding, thus affecting the enzyme functionality [56]. Additionally, the bulk of the glycyrrhizin structure favors interactions with Domain II and Domain III of the Mpro, in addition to flanking residues of the substrate-binding pocket, thus possibly affecting the dimerization of Domain I and Domain II and rendering the enzyme noncatalytic [57]. Similarly, Zhang et al. reported Mpro complex with an α-ketoamide inhibitor. The cocrystalized lead identified binds to the same substrate binding site reported in this study [28]. Morever, calpain inhibitors and GC-376 analogs are also confirmed to accommodate in the same functional pocket [58]. Beside these, many in silico studies have demonstrated the binding affinity of drug molecules to this active side of Mpro [33,59,60,61].

2.1.2. PLpro–Glycyrrhizin Complex

The PLpro enzyme of SARS-CoV-2 is implicated in viral polyproteins processing that generate a replicase complex and assist in virus spreading. The enzyme also plays a fundamental role in cleaving post-translational proteinaceous modifications present on the host protein as a mechanism to avoid antiviral host immune responses [22]. The docked complex between PLpro and glycyrrhizin highlighted the compound binding at the central palm catalytic cavity (Figure 5). Good binding of the compound-rich electronegative oxygen in the (2S,3S,4S,5R,6R)-6-(((2S,3R,4S,5S,6S)-6-carboxy-2,4,5-trihydroxytetrahydro-2H-pyran-3-yl)oxy)-3,4,5-trihydroxytetrahydro-2H-pyran-2-carboxylic acid at the docked site is the output of several strong hydrogen bond interactions: Gln174, Asp179, and Asn128. Besides these residues, the compounds moiety also formed van der Waals interaction, critical from a stability perspective. The remainder of the compound structure produced van der Waals contacts at this central cavity. The preferred binding of glycyrrhizin is at the central palm, sandwiching the finger and thumb domains, adjacent to the active substrate-binding pocket, which makes a strong bond with many vital catalytic residues. In contrast to the cocrystallized peptide inhibitor VIR251, which has a different conformation and binds to a different substrate cavity site, the glycyrrhizin-binding site is close to the VIR251 site [62]. In terms of interacting binding residues, the glycyrrhizin correlates more with the GRL0617 inhibitor of SAR-CoV-2 PLpro [63]. Further, the effect of conformational change of the BL2 loop upon glycyrrhizin binding is important to evaluate in future studies to disclose the glycyrrhizin recognition mechanism.
In literature, many inhibitors of coronaviruses PLpro are documented that include zinc conjugate inhibitors, naphthalene, and thiopurine derivatives, and natural products [64]. These molecules are known to interact with the active site residues reported in this study. Tanshinones are reported to show inhibition of deubiquitinase and proteolytic activitiy of SARS-CoV PLpro [65]; 8-(Trifluoromethyl)-9H-purin-6-amine is a reversible noncovalent inhibitor, whereas N-Ethylmaleimide (NEM) modifies SARS-CoV PLpro Cys [63]. Moreover, 6-mercaptopurine (6MP) and 6-thioguanine (6TG) are slow and competitive inhibitors that form hydrogen bonds with catalytic residues of the SAR-CoV PLpro [66]. Several in silico studies also demonstrated a range of compounds that interfere with the functional site of SARS-CoV-2 PLpro [67,68,69,70].

2.1.3. Nucleocapsid–Glycyrrhizin Complex

The SARS-CoV-2 N protein is an RNA binding protein and offers several functions of viral transcription and replication [20]. It particularly plays a pivotal role in helical ribonucleoprotein packing during RNA genome packing, regulating RNA replication, and modulating infected cell metabolism. Blocking of this protein could lead to blocking viral replication, and thus an attractive target for drug development. The compound glycyrrhizin was found to prefer docking at the loop region 1 at the junction between the β-sheet core and β-hairpin (Figure 6). The molecule is aligned perfectly along the cavity volume where its (2S,3S,4S,5R,6R)-6-(((2S,3R,4S,5S,6S)-6-carboxy-2,4,5-trihydroxytetrahydro-2H-pyran-3-yl)oxy)-3,4,5-trihydroxytetrahydro-2H-pyran-2-carboxylic acid part is connected to the β3 and β4 sheets of the β-hairpin. Here, this chemical moiety is involved in hydrogen bonding with Thr92, Arg94, and Arg89, and van der Waals contact with Arg90 and Ala91. The (2S,4aS,6aS,6bR,8aS,12aS,12bR,14bR)-2,4a,6a,6b,9,9,12a-heptamethyl-13-oxo-1,2,3,4,4a,5,6,6a,6b,7,8,8a,9,10,11,12,12a,12b,13,14b-icosahydropicene-2-carboxylic acid region of the compound produced hydrogen bonding with residues (Tyr110 and Arg150) and van der Waal contacts with residues (Asn49,Thr50, Als51, Phe54, Thr55, Tyr112, Pro118, Pro152, and Ala157) of β1, β2, β4, β5, β6, and β7 of the β-sheet core of the protein. Bhowmik et al. reported strong binding of Rutin, Doxycycline, Caffeic acid, Ferulic acid, Simeprevir, and Grazoprevir with several functional residues of the SARS-CoV-2 nucleocapsid protein reported in this study [71].

2.2. MD Simulation Analysis

In computer-aided drug design, MD simulations are essential in providing detailed biomolecule dynamical structural information and surface wealth of protein–ligand interactions, energetic data that are foremost to understanding the structural–functionality relationship of target protein principle in ligand recognition/interactions [37,72,73]. This set of information has tremendous applications in guiding novel drug design, thereby making MD simulation a successful tool in the modern drug discovery framework.

2.2.1. Root-Mean-Square Deviation (RMSD) Analysis

MD simulation of 50 ns was performed for each receptor with bound glycyrrhizin to elucidate the compound binding stability and extract receptors/compound structural information that is key in the binding that may be altered to iMprove binding conformation and, ultimately, compound affinity for the target biomolecules. First, RMSD of receptors in each complex was estimated as carbon alpha deviations by superimposing 50,000 snapshots over the initial reference structure versus time (Figure 7A). RMSDs of all three complexes were found: Mpro (maximum, 3.14 Å; mean, 1.97 Å), PLpro (maximum, 2.59 Å; mean, 1.64 Å), and Nucleocapsid (maximum, 2.34 Å; mean, 1.32 Å). All of the receptors are relatively stable in terms of 3D structure, and no flexibility in secondary structures was noticed. As a consequence, glycyrrhizin binding pose was not altered, thus reflecting strong and stable complex formation.

2.2.2. Glycyrrhizin Conformation Stability

In addition, the MD simulation trajectories were examined to disclose information about the glycyrrhizin conformation stability with the receptors (Figure 7B). The glycyrrhizin RMSD with the receptors is Mpro (maximum, 2.56 Å; mean, 0.93 Å), PLpro (maximum, 2.14 Å; mean, 0.96 Å), and Nucleocapsid (maximum, 4.20 Å; mean, 3.48 Å). The molecules disclosed high stable, except for some deviations in the glycyrrhizin binding mode with the Nucleocapsid protein; therefore, the end MD simulation snapshot over the initial was superimposed to understand the compound dynamics. The (2S,3S,4S,5R,6R)-6-(((2S,3R,4S,5S,6S)-6-carboxy-2,4,5-trihydroxytetrahydro-2H-pyran-3-yl)oxy)-3,4,5-trihydroxytetrahydro-2H-pyran-2-carboxylic acid fragment of the glycyrrhizin is flexible in an attempt to establish a more stable conformation. This moiety left its original site of interaction and moved more towards the β-core sheet for binding (Figure 8).

2.2.3. Root-Mean-Square Fluctuation (RMSF) Analysis

The residual flexibility and stability of the receptors in the presence of glycyrrhizin were further elucidated (Figure 7C). Mean RMSF for Mpro is 1.4 Å, PLpro is 1.57 Å, and Nucleocapsid is 1.9 Å. These values suggest good agreement on intermolecular stability.

2.2.4. Radius of Gyration (Rg) Analysis

Additionally, Rg analysis was performed to evaluate protein compactness and structural equilibrium over the simulation time (Figure 7D). The Rg of the systems follows: Mpro–glycyrrhizin (45.62 Å and 42.28 Å), PLpro–glycyrrhizin (50.29 Å and 46.23 Å), and Nucleocapsid–glycyrrhizin (35.71 Å and 30.70 Å). All three systems are quite stable and remain compact.

2.3. MMGB/PBSA Analysis

To get a deeper insight into the compounds binding potential with the SARS-CoV-2 enzymes used, binding free energies were estimated using MMGBSA and MMPBSA techniques. Additionally, per residue decomposition assay was accomplished to highlight residues that contribute majorly to the compound’s stability at the docked position and, ultimately, to the strong intermolecular interactions. To this objective, 100 frames were picked at time intervals of 50 ps from the simulation trajectories, discarding the water molecules and counterions. Detailed binding energies of the complexes are listed in Table 1 All of the binding interactions are energetically favorable, resulting in the formation of stable complexes. In all of the complexes, gas-phase energy dominates the system energy with significant contribution from van der Waals compared to electrostatic energy’s minor role. The polar solvation energy is illustrated to play a nonfavorable part in binding, whereas the nonpolar energy seems to be vital in complex equilibration. The MMGBSA net binding-energy-ranked stability of the complexes follows: PLpro–glycyrrhizin > Spike–glycyrrhizin > Nucleocapsid–glycyrrhizin > Mpro–glycyrrhizin. The MMPBSA ranking follows: PLpro–glycyrrhizin > Spike–glycyrrhizin > Mpro–glycyrrhizin > N–glycyrrhizin.

2.4. Per-Residue Decomposition

The atomic-level contribution of each residue from the enzymes to the compound binding was elucidated further. Those with an average binding energy of <1 kcal/mol were categorized as hotspot residues because of their significant overall complex stability contribution [74,75]. In the case of Mpro–glycyrrhizin interaction, Asn238 and Asp289 are vital in holding the compound at the docked site. Phe69, Asn128, Gln174, and Asp179 residues are critical in bridging PLpro enzyme with glycyrrhizin compound. The primary hotspot residues in Nucleocapsid–glycyrrhizin complex are Thr92, Arg94, Tyr110, and Arg150. It was further noticed that the van der Waals energy, as noted earlier, dominates the overall binding interaction energy. Hotspot residues of each receptor that are in direct contact and key in the stabilization of glycyrrhizin are presented in Table 2.

2.5. WaterSwap Binding Energy

WaterSwap uses an explicit solvation system that considers interaction details of protein–water, protein–water–ligand, and ligand–water. Such information is not provided in the MMGB/PBSA; therefore, it is not reliable for predicting the role of water molecules in biomolecule–ligand interactions [76]. Specifically, this holds great importance in an instance where the ligand is bridged to the receptor through water molecules. The WaterSwap method has been successfully applied to various biological systems and proved critical in determining absolute binding free energy. For each complex, the WaterSwap energies converged significantly after running 1000 frames. All the values also concluded good stability of intermolecular docked conformation. WaterSwap energies for each complex are shown in Table 3.

3. Materials and Methods

3.1. Target Proteins Preparation

The anti-SARS-CoV-2 targets (Mpro PDB code: 7BQY, PLpro PDB code: 6XAA, and Nucleocapsid PDB code: 6M3M) were retrieved and prepared using the AMBER18 program [77]. Ff14SB force field [78] was used for amino acid parameterization. To add complementary hydrogen atoms missed by the crystallography, the tleap module of AmberTools18 was employed. Energy minimization of the targeted proteins was done first for 1000 steepest descent steps, and then by 500 conjugate gradient steps, allowing the step size to be 0.02 Å. Charge addition was done through the Gasteiger method.

3.2. Compound Preparation

The MPD3 phytochemical database (https://www.bioinformation.info/), in addition to reported natural antiviral compounds, were used in this study to filter molecules that show best binding affinity to the selected SARS-CoV-2 multiple targets. The library containing ~5000 natural compounds was imported to PyRx 0.8 software [79], where they were minimized for optimal energy and followed by conversion to pdbqt format for use in virtual screening against the mentioned targets.

3.3. Structure-Based Virtual Screening

Virtual screening of the compounds against of the targets used was done using the AutoDock Vina in PyRx [80] on Windows 10-supported Dell system (processor: Intel(R) Core(TM) i7-8550U CPU @ 1.80 GHz with a 64-bit operating system, ×64-based processor, a memory of 8.00 GB). First, the docking protocol was validated by docking cocrystallized ligands to the protein keeping the docking parameters default except for the sphere around the binding site, which was set to 15 Å. Validation was also done by comparing the best-ranked compounds conformation relative to the crystallized ligand by root-mean-square deviation (RMSD) [81]. Docking of the compound to the targets was accomplished by using the same set of parameters described for the validation procedure and run in triplicates to absolute consistency of the results. The docked solutions were clustered, considering an RMSD value of 1 Å. The binding mode of compounds with the lowest binding energy in kcal/mol was refined in MD simulations.

3.4. MD Simulations

MD simulations of the docked solutions were performed using AMBER18 [77]. Each top complex was explicitly solvated with water molecules, and then to get a neutral system, counter ions were added. Afterward, using the TIP3P solvent model, a water box of thickness 12 Å was created to surround the complex [82]. Simulation of the complex was done through periodic boundary conditions where electrostatic interactions were modeled with the particle–mesh Ewald procedure [74]. In the process, a threshold value of 8 Å was defined for nonbounded interactions. Water molecules were minimized for 500 cycles, followed by complete system minimization for 1000 rounds. Then, each system temperature was gradually scaled to 300 K. Equilibration of the systems was achieved under the NPT ensemble for 100 ps. This involves equilibration of both counter ions and water molecules while considering restraint on solutes in the first phase for 50 ps; subsequent protein side chains were relaxed. MD simulation of 50 ns was performed at 300 K and 1 atm for two fs under the NPT ensemble. Hydrogen and covalent bonds were constrained using the SHAKE algorithm [83], whereas system temperature was controlled through Langevin dynamics [84]. The initial structure was used as a reference, and CPPTRAJ [85] of AMBER was run to generate a root-mean-square deviation (RMSD) plot to check the system MD simulation convergence [81]. Ligand structural flexibilities were calculated by ligand RMSD. Furthermore, hydrogen bond analysis was performed to investigate hydrogen bonds formed between the compounds and amino acids present within the docked site vicinity.

3.5. MMGB/PBSA Analysis

The binding free energy (ΔG binding) of the complexes was estimated using the AMBER18 MM/PBSA method [42,86]. One hundred snapshots were considered from simulation trajectories at a regular time interval to calculate the free energy difference.
ΔGbinding = Gcomplex − (Gprotein + Gligand)
ΔG = ΔGgas + ΔGsolv − TΔS
ΔGgas = Δele + ΔGvdw
ΔGsolv = ΔGGB + ΔGSA
ΔGSA = γ × SASA × b
In these equations, Gcomplex is delta free energy of the complex, Gprotein is delta free energy of the protein, and Gligand is delta free energy of the ligand; ΔGgas represents gas-phase energy and can be split into delta electrostatic (ΔEele), and delta van der Waals (ΔEvdw) energy; and the ΔGsolv term stands for solvation free energy, which comprises polar (ΔGGB) and nonpolar (ΔGSA) energy. In the ΔGGB, the εw value is set to 80, and εp is selected as 1.0. Linear combinations of the pairwise overlap method are used to estimate the solvent-accessible surface area (SASA).

3.6. WaterSwap Analysis

WaterSwap [76,87] was additionally done over the last 10 ns of MD simulation for a total of default 1000 iterations, keeping the sample size of Monte Carlo simulation to 1.6 × 109. The absolute binding energy of each complex was estimated using three useful algorithms: thermodynamics integration, free energy perturbation, and Bennett’s. The energy value <1 kcal/mol represents a good convergence of the system [75].

4. Conclusions

In this study, we found glycyrrhizin as the most significant natural compound that can act as a double-edged sword and inhibit multiple proteins of SARS-CoV-2. This compound has a high binding affinity for all of the SARS-CoV-2 receptors used in this study and had a stable binding mode in the MD simulation time. The compound revealed important interactions with all receptors, and thus requires further consideration in future anti-SARS-CoV-2 therapeutic studies. Glycyrrhizin has been previously documented to have therapeutic applications against SARS-CoV, chronic hepatitis C, and HIV-1 [88]. The molecule is clinically useful and had few toxic reactions. One way to overcome toxicity is by allowing low concentration of the drug in the cells (<100 µg/mL) [89]. Glycyrrhizin has been reported to inhibit viral penetration and effective both during the viral infection and postinfection [90]. It was previously demonstrated that the glycyrrhizin binds with good affinity to the human ACE2 and interacts with Asp30, Gln288, Arg393, and Arg559 residues, hence also underlines its potential to target the SARS-CoV-2 Spike protein RBD attachment to the human ACE2 receptor [90]. It also was shown that glycyrrhizin can be employed in synergism along with other plant-based molecules to treat SARS-CoVs [91]. From a pharmacological perspective, the glycyrrhizin prevents the production of intracellular reactive oxygen species, activates interferon production, downregulates proinflammatory cytokines, lowers airway exudate production, and inhibits thrombin [45,92]. The compound was also computationally characterized previously to bind with good affinity to SARS-CoV-2 main protease [93]. Therefore, additional structural modification to lower the side effects and enhance the clinical efficacy of this compound is of high interest to treat SARS-related infections.

Author Contributions

Conceptualization, G.L.; data curation, Z.T.M., A.R.H., and H.M.H.A.-H.; funding acquisition, G.L.; investigation, Z.T.M.; project administration, G.L.; software, S.A.; supervision, G.L.; validation, A.R.H., H.M.H.A.-H., and S.A.; visualization, Z.T.M.; writing—original draft, Z.T.M.; writing—review and editing, A.R.H., H.M.H.A.-H., S.A., and G.L. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Natural Science Foundation of China (grant numbers: 31770333, 31370329, and 11631012).

Data Availability Statement

The data presented in this study are available within the article.

Acknowledgments

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

Conflicts of Interest

The authors declare that the research was conducted in the absence of any commercial or finan-cial relationships that could be construed as a potential conflict of interest.

Sample Availability

Not available.

References

  1. Yang, Y.; Peng, F.; Wang, R.; Guan, K.; Jiang, T.; Xu, G.; Sun, J.; Chang, C. The Deadly Coronaviruses: The 2003 SARS Pandemic and the 2020 Novel Coronavirus Epidemic in China. J. Autoimmun. 2020, 109, 102434. [Google Scholar] [CrossRef] [PubMed]
  2. Wu, F.; Zhao, S.; Yu, B.; Chen, Y.-M.; Wang, W.; Song, Z.-G.; Hu, Y.; Tao, Z.-W.; Tian, J.-H.; Pei, Y.-Y.; et al. A New Coronavirus Associated with Human Respiratory Disease in China. Nature 2020, 579, 265–269. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  3. Hui, D.S.; Azhar, E.I.; Madani, T.A.; Ntoumi, F.; Kock, R.; Dar, O.; Ippolito, G.; Mchugh, T.D.; Memish, Z.A.; Drosten, C.; et al. The Continuing 2019-nCoV Epidemic Threat of Novel Coronaviruses to Global Health—The latest 2019 Novel Coronavirus Outbreak in Wuhan, China. Int. J. Infect. Dis. 2020, 91, 264–266. [Google Scholar] [CrossRef] [Green Version]
  4. Ye, Z.-W.; Yuan, S.; Yuen, K.-S.; Fung, S.-Y.; Chan, C.-P.; Jin, D.-Y. Zoonotic Origins of Human Coronaviruses. Int J. Biol. Sci. 2020, 16, 1686–1697. [Google Scholar] [CrossRef] [Green Version]
  5. Tahir ul Qamar, M.; Saleem, S.; Ashfaq, U.A.; Bari, A.; Anwar, F.; Alqahtani, S. Epitope-Based Peptide Vaccine Design and Target Site Depiction against Middle East Respiratory Syndrome Coronavirus: An Immune-Informatics Study. J. Transl. Med. 2019, 17, 362. [Google Scholar] [CrossRef]
  6. Kim, D.; Lee, J.-Y.; Yang, J.-S.; Kim, J.W.; Kim, V.N.; Chang, H. The Architecture of SARS-CoV-2 Transcriptome. Cell 2020, 181, 914–921.e10. [Google Scholar] [CrossRef]
  7. Zhu, N.; Zhang, D.; Wang, W.; Li, X.; Yang, B.; Song, J.; Zhao, X.; Huang, B.; Shi, W.; Lu, R.; et al. A Novel Coronavirus from Patients with Pneumonia in China, 2019. N. Engl. J. Med. 2020, 382, 727–733. [Google Scholar] [CrossRef]
  8. Tahir ul Qamar, M.; Alqahtani, S.M.; Alamri, M.A.; Chen, L.-L. Structural Basis of SARS-CoV-2 3CLpro and Anti-COVID-19 Drug Discovery from Medicinal Plants. J. Pharm. Anal. 2020, 10, 313–319. [Google Scholar] [CrossRef] [PubMed]
  9. Zhou, P.; Yang, X.-L.; Wang, X.-G.; Hu, B.; Zhang, L.; Zhang, W.; Si, H.-R.; Zhu, Y.; Li, B.; Huang, C.-L.; et al. A Pneumonia Outbreak Associated with a New Coronavirus of Probable Bat Origin. Nature 2020, 579, 270–273. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  10. Zhang, Y.-Z.; Holmes, E.C. A Genomic Perspective on the Origin and Emergence of SARS-CoV-2. Cell 2020, 181, 223–227. [Google Scholar] [CrossRef]
  11. Fahmi, M.; Kubota, Y.; Ito, M. Nonstructural Proteins NS7b and NS8 are Likely to be Phylogenetically Associated with Evolution of 2019-nCoV. Infect. Genet. Evol. 2020, 81, 104272. [Google Scholar] [CrossRef]
  12. Tahir ul Qamar, M.; Rehman, A.; Tusleem, K.; Ashfaq, U.A.; Qasim, M.; Zhu, X.; Fatima, I.; Shahid, F.; Chen, L.-L. Designing of a Next Generation Multiepitope Based Vaccine (MEV) against SARS-COV-2: Immunoinformatics and in Silico Approaches. PLoS ONE 2020, 15, e0244176. [Google Scholar] [CrossRef]
  13. Lenzen, M.; Li, M.; Malik, A.; Pomponi, F.; Sun, Y.-Y.; Wiedmann, T.; Faturay, F.; Fry, J.; Gallego, B.; Geschke, A.; et al. Global Socio-Economic Losses and Environmental Gains from the Coronavirus Pandemic. PLoS ONE 2020, 15, e0235654. [Google Scholar] [CrossRef]
  14. Boopathi, S.; Poma, A.B.; Kolandaivel, P. Novel 2019 Coronavirus Structure, Mechanism of Action, Antiviral drug Promises and Rule Out against Its Treatment. J. Biomol. Struct. Dyn. 2020, 1–14. [Google Scholar] [CrossRef] [Green Version]
  15. Tahir ul Qamar, M.; Shahid, F.; Aslam, S.; Ashfaq, U.A.; Aslam, S.; Fatima, I.; Fareed, M.M.; Zohaib, A.; Chen, L.-L. Reverse Vaccinology Assisted Designing of Multiepitope-Based Subunit Vaccine Against SARS-CoV-2. Infect. Dis. Poverty 2020, 9, 1–14. [Google Scholar] [CrossRef] [PubMed]
  16. Yi, C.; Sun, X.; Ye, J.; Ding, L.; Liu, M.; Yang, Z.; Lu, X.; Zhang, Y.; Ma, L.; Gu, W.; et al. Key Residues of the Receptor Binding Motif in the Spike Protein of SARS-CoV-2 That Interact with ACE2 and Neutralizing Antibodies. Cell. Mol. Immunol. 2020, 17, 621–630. [Google Scholar] [CrossRef]
  17. Chen, R.; Fu, J.; Hu, J.; Li, C.; Zhao, Y.; Qu, H.; Wen, X.; Cao, S.; Wen, Y.; Wu, R.; et al. Identification of the Immunodominant Neutralizing Regions in the Spike Glycoprotein of Porcine Deltacorona-Virus. Virus Res. 2020, 276, 197834. [Google Scholar] [CrossRef] [PubMed]
  18. Gui, M.; Song, W.; Zhou, H.; Xu, J.; Chen, S.; Xiang, Y.; Wang, X. Cryo-Electron Microscopy Structures of the SARS-CoV Spike Glycoprotein Reveal a Prerequisite Conformational State for Receptor Binding. Cell Res. 2017, 27, 119–129. [Google Scholar] [CrossRef]
  19. Jin, Z.; Zhao, Y.; Sun, Y.; Zhang, B.; Wang, H.; Wu, Y.; Zhu, Y.; Zhu, C.; Hu, T.; Du, X.; et al. Structural Basis for the Inhibition of SARS-CoV-2 Main Protease by Antineoplastic Drug Carmofur. Nat. Struct. Mol. Biol. 2020, 27, 529–532. [Google Scholar] [CrossRef] [PubMed]
  20. Kang, S.; Yang, M.; Hong, Z.; Zhang, L.; Huang, Z.; Chen, X.; He, S.; Zhou, Z.; Zhou, Z.; Chen, Q.; et al. Crystal Structure of SARS-CoV-2 Nucleocapsid Protein RNA Binding Domain Reveals Potential Unique Drug Tar-Geting Sites. Acta Pharm. Sin. B 2020, 19, 1228–1238. [Google Scholar] [CrossRef]
  21. Elfiky, A.A. SARS-CoV-2 RNA Dependent RNA Polymerase (RdRp) Targeting: An in Silico Perspective. J. Biomol. Struct. Dyn. 2020, 1–15. [Google Scholar] [CrossRef]
  22. Shin, D.; Mukherjee, R.; Grewe, D.; Bojkova, D.; Baek, K.; Bhattacharya, A.; Schulz, L.; Widera, M.; Mehdipour, A.R.; Tascher, G.; et al. Papain-Like Protease Regulates SARS-CoV-2 Viral Spread and Innate Immunity. Nature 2020, 587, 657–662. [Google Scholar] [CrossRef] [PubMed]
  23. Gyebi, G.A.; Ogunro, O.B.; Adegunloye, A.P.; Ogunyemi, O.M.; Afolabi, S.O. Potential Inhibitors of Coronavirus 3-Chymotrypsin-Like Protease (3CLpro): An in Silico Screening of Alkaloids and Terpenoids from African Medicinal Plants. J. Biomol. Struct. Dyn. 2020, 1–19. [Google Scholar] [CrossRef] [PubMed]
  24. Tomar, P.P.S.; Arkin, I.T. SARS-CoV-2 E Protein Is a Potential Ion Channel That Can Be Inhibited by Gliclazide and Memantine. Biochem. Biophys. Res. Commun. 2020, 530, 10–14. [Google Scholar] [CrossRef] [PubMed]
  25. Yan, R.; Zhang, Y.; Li, Y.; Xia, L.; Guo, Y.; Zhou, Q. Structural Basis for the Recognition of SARS-CoV-2 by full-Length Human ACE2. Science 2020, 367, 1444–1448. [Google Scholar] [CrossRef] [Green Version]
  26. Hirano, T.; Murakami, M. COVID-19: A New Virus, but a Familiar Receptor and Cytokine Release Syndrome. Immunity 2020, 52, 731–733. [Google Scholar] [CrossRef]
  27. Du, L.; Zhao, G.; Lin, Y.; Chan, C.; He, Y.; Jiang, S.; Wu, C.; Jin, D.-Y.; Yuen, K.-Y.; Zhou, Y.; et al. Priming with rAAV Encoding RBD of SARS-CoV S Protein and Boosting with RBD-Specific Peptides for T Cell Epitopes Elevated Humoral and Cellular Immune Responses Against SARS-CoV Infection. Vaccine 2008, 26, 1644–1651. [Google Scholar] [CrossRef] [PubMed]
  28. Zhang, L.; Lin, D.; Sun, X.; Curth, U.; Drosten, C.; Sauerhering, L.; Becker, S.; Rox, K.; Hilgenfeld, R. Crystal Structure of SARS-CoV-2 Main Protease Provides a Basis for Design of Improved $α$-Ketoamide Inhibitors. Science 2020, 368, 409–412. [Google Scholar] [CrossRef] [Green Version]
  29. Kumar, Y.; Singh, H.; Patel, C.N. In Silico Prediction of Potential Inhibitors for the Main Protease of SARS-CoV-2 Using Molecular Docking and Dynamics Simulation Based Drug-Repurposing. J. Infect. Public Health 2020, 13, 1210–1223. [Google Scholar] [CrossRef]
  30. Alamri, M.A.; Tahir ul Qamar, M.; Mirza, M.U.; Bhadane, R.; Alqahtani, S.M.; Muneer, I.; Froeyen, M.; Salo-Ahen, O.M.H. Pharmacoinformatics and Molecular Dynamics Simulation Studies Reveal Potential Covalent and FDA-Approved Inhibitors of SARS-CoV-2 Main Protease 3CLpro. J. Biomol. Struct. Dyn. 2020, 1–13. [Google Scholar] [CrossRef]
  31. Dai, W.; Zhang, B.; Jiang, X.-M.; Su, H.; Li, J.; Zhao, Y.; Xie, X.; Jin, Z.; Peng, J.; Liu, F.; et al. Structure-Based Design of Antiviral Drug Candidates Targeting the SARS-CoV-2 Main Protease. Science 2020, 368, 1331–1335. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  32. Joshi, R.S.; Jagdale, S.S.; Bansode, S.B.; Shankar, S.S.; Tellis, M.B.; Pandya, V.K.; Chugh, A.; Giri, A.P.; Kulkarni, M.J. Discovery of Potential Multi-Target-Directed Ligands by Targeting Host-Specific SARS-CoV-2 Structurally Conserved Main Proteases. J. Biomol. Struct. Dyn. 2020, 1–16. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  33. Havranek, B.; Islam, S.M. An in Silico Approach for Identification of Novel Inhibitors as Potential Therapeutics Targeting COVID-19 Main Protease. J. Biomol. Struct. Dyn. 2020, 1–12. [Google Scholar] [CrossRef] [PubMed]
  34. Alamri, M.A.; Tahir ul Qamar, M.; Mirza, M.U.; Alqahtani, S.M.; Froeyen, M.; Chen, L.-L. Discovery of Human Coronaviruses Pan-Papain-like Protease Inhibitors Using Computational Approaches. J. Pharm. Anal. 2020, 10, 546–559. [Google Scholar] [CrossRef] [PubMed]
  35. Riaz, M.; Ashfaq, U.A.; Qasim, M.; Yasmeen, E.; Tahir ul Qamar, M.; Anwar, F. Screening of Medicinal Plant Phytochemicals as Natural Antagonists of p53-MDM2 Interaction to Reactivate p53 Functioning. Anticancer Drugs 2017, 28, 1032–1038. [Google Scholar] [CrossRef] [PubMed]
  36. Rehan Khalid, R.; Tahir ul Qamar, M.; Maryam, A.; Ashique, A.; Anwar, F.H.; Geesi, M.; Siddiqi, A.R. Comparative Studies of the Dynamics Effects of BAY60-2770 and BAY58-2667 Binding with Human and Bacterial H-NOX Domains. Molecules 2018, 23, 2141. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  37. Durdagi, S.; Tahir ul Qamar, M.; Salmas, R.E.; Tariq, Q.; Anwar, F.; Ashfaq, U.A. Investigating the Molecular Mechanism of Staphylococcal DNA Gyrase Inhibitors: A Combined Ligand-Based and Structure-Based Resources Pipeline. J. Mol. Graph. Model. 2018, 85, 122–129. [Google Scholar] [CrossRef]
  38. Muneer, I.; Tusleem, K.; Abdul Rauf, S.; Hussain, H.M.J.; Siddiqi, A.R. Others Discovery of Selective Inhibitors for Cyclic AMP Response Element-Binding Protein: A Combined Ligand and Structure-Based Resources Pipeline. Anti-Cancer Drugs 2019, 30, 363–373. [Google Scholar]
  39. Yu, W.; MacKerell, A.D. Computer-Aided Drug Design Methods. In Antibiotics; Springer: Berlin/Heidelberg, Germany, 2017; pp. 85–106. [Google Scholar]
  40. Tahir ul Qamar, M.; Ashfaq, U.A.; Tusleem, K.; Mumtaz, A.; Tariq, Q.; Goheer, A.; Ahmed, B. In-Silico Identification and Evaluation of Plant Flavonoids as Dengue NS2B/NS3 Protease Inhibitors Using Molecular Docking and Simulation Approach. Pak. J. Pharm. Sci. 2017, 30, 2119–2137. [Google Scholar]
  41. Durrant, J.D.; McCammon, J.A. Molecular Dynamics Simulations and Drug Discovery. BMC Biol. 2011, 9, 71. [Google Scholar] [CrossRef] [Green Version]
  42. Genheden, S.; Ryde, U. The MM/PBSA and MM/GBSA Methods to Estimate Ligand-Binding Affinities. Expert Opin. Drug Discov. 2015, 10, 449–461. [Google Scholar] [CrossRef] [PubMed]
  43. Mumtaz, A.; Ashfaq, U.A.; Qamar, U.M.T.; Anwar, F.; Gulzar, F.; Ali, M.A.; Saari, N.; Pervez, M.T. MPD3: A Useful Medicinal Plants Database for Drug Designing. Nat. Prod. Res. 2017, 31, 1228–1236. [Google Scholar] [CrossRef] [PubMed]
  44. Santos, I.d.A.; Grosche, V.R.; Bergamini, F.R.G.; Sabino-Silva, R.; Jardim, A.C.G. Antivirals against Coronaviruses: Candidate Drugs for SARS-CoV-2 Treatment? Front. Microbiol. 2020, 11, 1818. [Google Scholar] [CrossRef] [PubMed]
  45. LuoLiu, P.; Li, J. Pharmacologic Perspective: Glycyrrhizin May Be an Efficacious Therapeutic Agent for COVID-19. Int. J. Antimicrob. Agents 2020, 105995. [Google Scholar]
  46. Elshabrawy, H.A. SARS-CoV-2: An Update on Potential Antivirals in Light of SARS-CoV Antiviral Drug Discoveries. Vaccines 2020, 8, 335. [Google Scholar] [CrossRef]
  47. Kato, F.; Matsuyama, S.; Kawase, M.; Hishiki, T.; Katoh, H.; Takeda, M. Antiviral Activities of Mycophenolic Acid and IMD-0354 Against SARS-CoV-2. Microbiol. Immunol. 2020, 64, 635–639. [Google Scholar] [CrossRef]
  48. Case, D.; Ben-Shalom, I.; Brozell, S.; Cerutti, D.; Cheatham III, T.; Cruzeiro, V.; Darden, T.; Duke, R.; Ghoreishi, D.; Gilson, M.; et al. AMBER 18; University of California: San Francisco, CA, USA, 2018. [Google Scholar]
  49. Maier, J.A.; Martinez, C.; Kasavajhala, K.; Wickstrom, L.; Hauser, K.E.; Simmerling, C. ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters From ff99SB. J. Chem. Theory Comput. 2015, 11, 3696–3713. [Google Scholar] [CrossRef] [Green Version]
  50. Dallakyan, S.; Olson, A.J. Small-Molecule Library Screening by Docking with PyRx. In Chemical Biology; Springer: Berlin/Heidelberg, Germany, 2015; pp. 243–250. [Google Scholar]
  51. Trott, O.; Olson, A.J. Auto Dock Vina: Improving the Speed and Accuracy of Docking with a New Scoring Function, Efficient Optimization, and Multi-Threading. J. Comput. Chem. 2010, 31, 455–461. [Google Scholar]
  52. Maiorov, V.N.; Crippen, G.M. Significance of Root-Mean-Square Deviation in Comparing Three-Dimensional Structures of Globular Proteins. J. Mol. Biol. 1994, 235, 625–634. [Google Scholar] [CrossRef] [Green Version]
  53. Andleeb, S.; Imtiaz-Ud-Din; Rauf, M.K.; Azam, S.S.; Badshah, A.; Sadaf, H.; Raheel, A.; Tahir, M.N.; Raza, S. A One-Pot Multicomponent Facile Synthesis of Dihydropyrimidin-2(1: H)-Thione Derivatives Using Triphenyl-Germane as a Catalyst and Its Binding Pattern Validation. RSC Adv. 2016, 6, 79651–79661. [Google Scholar] [CrossRef]
  54. Abro, A.; Azam, S.S. Binding Free Energy Based Analysis of Arsenic (+3 Oxidation State) Methyltransferase with S-Adenosylmethionine. J. Mol. Liq. 2016, 220, 375–382. [Google Scholar] [CrossRef]
  55. Kräutler, V.; Van Gunsteren, W.F.; Hünenberger, P.H. A fast SHAKE Algorithm to Solve Distance Constraint Equations for Small Molecules in Molecular Dynamics Simulations. J. Comput. Chem. 2001, 22, 501–508. [Google Scholar] [CrossRef]
  56. Izaguirre, J.A.; Catarello, D.P.; Wozniak, J.M.; Skeel, R.D. Langevin Stabilization of Molecular Dynamics. J. Chem. Phys. 2001, 114, 2090–2098. [Google Scholar] [CrossRef]
  57. Roe, D.R.; Cheatham III, T.E. PTRAJ and CPPTRAJ: Software for Processing and Analysis of Molecular Dynamics Trajectory Data. J. Chem. Theory Comput. 2013, 9, 3084–3095. [Google Scholar] [CrossRef] [PubMed]
  58. Miller, B.R.; McGee, T.D.; Swails, J.M.; Homeyer, N.; Gohlke, H.; Roitberg, A.E. MMPBSA.py: An Efficient Program for End-State Free Energy Calculations. J. Chem. Theory Comput. 2012, 8, 3314–3321. [Google Scholar] [CrossRef]
  59. Woods, C.J.; Malaisree, M.; Michel, J.; Long, B.; McIntosh-Smith, S.; Mulholland, A.J. Rapid Decomposition and Visualisation of Protein-Ligand Binding Free Energies by Residue and by Water. Faraday Discuss. 2014, 169, 477–499. [Google Scholar] [CrossRef] [Green Version]
  60. Woods, C.J.; Malaisree, M.; Hannongbua, S.; Mulholland, A.J. A Water-Swap Reaction Coordinate for the Calculation of Absolute Protein-Ligand Binding Free Energies. J. Chem. Phys. 2011, 134. [Google Scholar] [CrossRef] [Green Version]
  61. Kiani, Y.S.; Ranaghan, K.E.; Jabeen, I.; Mulholland, A.J. Molecular Dynamics Simulation Framework to Probe the Binding Hypothesis of CYP3A4 Inhibitors. Int. J. Mol. Sci. 2019, 20, 4468. [Google Scholar] [CrossRef] [Green Version]
  62. Ahmed, B.; Ashfaq, U.A.; Tahir ul Qamar, M.; Ahmad, M. Anticancer Potential of Phytochemicals against Breast Cancer: Molecular Docking and Simulation Approach. Bangladesh J. Pharmacol. 2014, 9, 545–550. [Google Scholar] [CrossRef] [Green Version]
  63. Tahir ul Qamar, M.; Mumtaz, A.; Ashfaq, U.A.; Adeel, M.M.; Fatima, T. Potential of Plant Alkaloids as Dengue ns3 Protease Inhibitors: Molecular Docking and Simulation Approach. Bangladesh J. Pharmacol. 2014, 9, 262–267. [Google Scholar] [CrossRef] [Green Version]
  64. Ashfaq, U.A.; Jalil, A.; Tahir ul Qamar, M. Antiviral Phytochemicals Identification from Azadirachta Indica Leaves against HCV NS3 Protease: An in Sili-Co Approach. Nat. Prod. Res. 2016, 30, 1866–1869. [Google Scholar] [CrossRef] [PubMed]
  65. Tahir ul Qamar, M.; Kiran, S.; Ashfaq, U.A.; Javed, M.R.; Anwar, F.; Ali, M.A.; Gilani, A. ul H. Discovery of Novel Dengue NS2B/NS3 Protease Inhibitors Using Pharmacophore Modeling and Molecular Docking Based Virtual Screening of the Zinc Database. Int. J. Pharmacol. 2016, 12, 621–632. [Google Scholar] [CrossRef]
  66. Morris, G.M.; Lim-Wilby, M. Molecular Docking. In Molecular Modeling of Proteins; Springer: Berlin/Heidelberg, Germany, 2008; pp. 365–382. [Google Scholar]
  67. Muhseen, Z.T.; Hameed, A.R.; Al-Hasani, H.M.H.; Tahir ul Qamar, M.; Li, G. Promising Terpenes as SARS-CoV-2 Spike Receptor-Binding Domain (RBD) Attachment Inhibitors to the Human ACE2 Receptor: Integrated Computational Approach. J. Mol. Liq. 2020, 320, 114493. [Google Scholar] [CrossRef] [PubMed]
  68. Ton, A.-T.; Gentile, F.; Hsing, M.; Ban, F.; Cherkasov, A. Rapid Identification of Potential Inhibitors of SARS-CoV-2 Main Protease by Deep Docking of 1.3 Billion Compounds. Mol. Inform. 2020, 39, 2000028. [Google Scholar] [CrossRef] [Green Version]
  69. Jin, Z.; Du, X.; Xu, Y.; Deng, Y.; Liu, M.; Zhao, Y.; Zhang, B.; Li, X.; Zhang, L.; Peng, C.; et al. Structure of M Pro from SARS-CoV-2 and Discovery of its Inhibitors. Nature 2020, 582, 289–293. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  70. Xue, X.; Yu, H.; Yang, H.; Xue, F.; Wu, Z.; Shen, W.; Li, J.; Zhou, Z.; Ding, Y.; Zhao, Q.; et al. Structures of Two Coronavirus Main Proteases: Implications for Substrate Binding and Antiviral Drug Design. J. Virol. 2008, 82, 2515–2527. [Google Scholar] [CrossRef] [Green Version]
  71. Kneller, D.W.; Phillips, G.; O’Neill, H.M.; Jedrzejczak, R.; Stols, L.; Langan, P.; Joachimiak, A.; Coates, L.; Kovalevsky, A. Structural Plasticity of the SARS-CoV-2 3CL Mpro Active Site Cavity Revealed by Room Temperature X-ray Crystallography. Nat. Commun. 2020, 11, 3202. [Google Scholar] [CrossRef]
  72. Sacco, M.D.; Ma, C.; Lagarias, P.; Gao, A.; Townsend, J.A.; Meng, X.; Dube, P.; Zhang, X.; Hu, Y.; Kitamura, N.; et al. Structure and Inhibition of the SARS-CoV-2 Main Protease Reveal Strategy for Developing Dual Inhibitors against Mpro and Cathepsin L. Sci. Adv. 2020, 6, eabe0751. [Google Scholar] [CrossRef]
  73. Alberto, J.-A.; Ribas-Aparicio, R.M.; Ozores, A.G.; Vega, C.J.A. Virtual Screening of Approved Drugs as Potential SARS-CoV-2 Main Protease Inhibitors. Comput. Biol. Chem. 2020, 88, 107325. [Google Scholar] [CrossRef]
  74. Ngo, S.T.; Pham, Q.A.N.; Le, T.L.; Pham, D.-H.; Vu, V.V. Computational Determination of Potential Inhibitors of SARS-CoV-2 Main Protease. J. Chem. Inf. Model. 2020, 60, 5771–5780. [Google Scholar] [CrossRef]
  75. Wang, J. Fast Identification of Possible Drug Treatment of Coronavirus Disease-19 (COVID-19) through Computational Drug Repurposing Study. J. Chem. Inf. Model. 2020, 60, 3277–3286. [Google Scholar] [CrossRef] [PubMed]
  76. Gao, X.; Qin, B.; Chen, P.; Zhu, K.; Hou, P.; Wojdyla, J.A.; Wang, M.; Cui, S. Crystal structure of SARS-CoV-2 papain-like protease. Acta Pharm. Sin. B 2020, 11, 237–245. [Google Scholar] [CrossRef] [PubMed]
  77. Fu, Z.; Huang, B.; Tang, J.; Liu, S.; Liu, M.; Ye, Y.; Liu, Z.; Xiong, Y.; Cao, D.; Li, J.; et al. Structural Basis for the Inhibition of the Papain-Like Protease of SARS-CoV-2 by Small Molecules. Biorxiv 2020. [Google Scholar] [CrossRef]
  78. Petushkova, A.I.; Zamyatnin, A. A Papain-Like Proteases as Coronaviral Drug Targets: Current Inhibitors, Opportunities, and Limitations. Pharmaceuticals 2020, 13, 277. [Google Scholar] [CrossRef] [PubMed]
  79. Park, J.-Y.; Kim, J.H.; Kim, Y.M.; Jeong, H.J.; Kim, D.W.; Park, K.H.; Kwon, H.-J.; Park, S.-J.; Lee, W.S.; Ryu, Y.B. Tanshinones as Selective and Slow-Binding Inhibitors for SARS-CoV Cysteine Proteases. Bioorg. Med. Chem. 2012, 20, 5928–5935. [Google Scholar] [CrossRef] [PubMed]
  80. Chu, H.-F.; Chen, C.-C.; Moses, D.C.; Chen, Y.-H.; Lin, C.-H.; Tsai, Y.-C.; Chou, C.-Y. Porcine Epidemic Diarrhea Virus Papain-Like Protease 2 Can Be Noncompetitively Inhibited by 6-Thioguanine. Antiviral Res. 2018, 158, 199–205. [Google Scholar] [CrossRef]
  81. Kouznetsova, V.L.; Zhang, A.; Tatineni, M.; Miller, M.A.; Tsigelny, I.F. Potential COVID-19 Papain-like Protease PLpro Inhibitors: Repurposing FDA-Approved Drugs. PeerJ 2020, 8, e9965. [Google Scholar] [CrossRef]
  82. Amin, S.A.; Ghosh, K.; Gayen, S.; Jha, T. Chemical-Informatics Approach to COVID-19 Drug Discovery: Monte Carlo based QSAR, vIrtual Screening and Molecular Docking Study of Some in-House Molecules as Papain-Like Protease (PLpro) Inhibitors. J. Biomol. Struct. Dyn. 2020, 1–10. [Google Scholar] [CrossRef]
  83. Mirza, M.U.; Ahmad, S.; Abdullah, I.; Froeyen, M. Identification of Novel Human USP2 Inhibitor and its Putative Role in Treatment of COVID-19 by Inhibiting SARS-CoV-2 Papain-Like (PLpro) Protease. Comput. Biol. Chem. 2020, 89, 107376. [Google Scholar] [CrossRef]
  84. Bhati, S. Structure-Based Drug Designing of Naphthalene Based SARS-CoV PLpro Inhibitors for the Treatment of COVID-19. Heliyon 2020, 6, e05558. [Google Scholar] [CrossRef]
  85. Bhowmik, D.; Nandi, R.; Jagadeesan, R.; Kumar, N.; Prakash, A.; Kumar, D. Identification of Potential Inhibitors against SARS-CoV-2 by Targeting Proteins Responsible for Envelope for-Mation and Virion Assembly Using Docking Based Virtual Screening, and Pharmacokinetics Approaches. Infect. Genet. Evol. 2020, 84, 104451. [Google Scholar] [CrossRef] [PubMed]
  86. Tahir ul Qamar, M.; Maryam, A.; Muneer, I.; Xing, F.; Ashfaq, U.A.; Khan, F.A.; Anwar, F.; Geesi, M.H.; Khalid, R.R.; Rauf, S.A.; et al. Computational Screening of Medicinal Plant Phytochemicals to Discover Potent Pan-Serotype Inhibitors against Dengue Virus. Sci. Rep. 2019, 9, 1–16. [Google Scholar] [CrossRef] [PubMed]
  87. Karplus, M.; McCammon, J.A. Molecular Dynamics Simulations of Biomolecules. Nat. Struct. Mol. Biol. 2002, 9, 646. [Google Scholar] [CrossRef] [PubMed]
  88. Cinatl, J.; Morgenstern, B.; Bauer, G.; Chandra, P.; Rabenau, H.; Doerr, H.W. Glycyrrhizin, an Active Component of Liquorice Roots, and Replication of SARS-Associated Coronavirus. Lancet 2003, 361, 2045–2046. [Google Scholar] [CrossRef] [Green Version]
  89. Ashfaq, U.A.; Masoud, M.S.; Nawaz, Z.; Riazuddin, S. Glycyrrhizin as Antiviral Agent against Hepatitis C Virus. J. Transl. Med. 2011, 9, 1–7. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  90. Chen, H.; Du, Q. Potential Natural Compounds for Preventing SARS-CoV-2 (2019-nCoV) Infection. Preprints 2020. [Google Scholar] [CrossRef]
  91. Prasad, A.; Muthamilarasan, M.; Prasad, M. Synergistic Antiviral Effects against SARS-CoV-2 by Plant-Based Molecules. Plant Cell Rep. 2020, 39, 1109–1114. [Google Scholar] [CrossRef]
  92. Chrzanowski, J.; Chrzanowska, A.; Graboń, W. Glycyrrhizin: An Old Weapon against a Novel Coronavirus. Phyther. Res. 2020. [Google Scholar] [CrossRef]
  93. Pham, M.Q.; Vu, K.B.; Pham, T.N.H.; Tran, L.H.; Tung, N.T.; Vu, V.V.; Nguyen, T.H.; Ngo, S.T. Rapid Prediction of Possible Inhibitors for SARS-CoV-2 Main Protease Using Docking and FPL Simulations. RSC Adv. 2020, 10, 31991–31996. [Google Scholar] [CrossRef]
Figure 1. Schematic presentation of the methodology used in this current study.
Figure 1. Schematic presentation of the methodology used in this current study.
Molecules 26 00674 g001
Figure 2. AutoDock binding affinity score of the compounds to the SARS-CoV-2 enzymes.
Figure 2. AutoDock binding affinity score of the compounds to the SARS-CoV-2 enzymes.
Molecules 26 00674 g002
Figure 3. Two-dimensional presentation of high affinity binders to the SARS-CoV-2 proteins.
Figure 3. Two-dimensional presentation of high affinity binders to the SARS-CoV-2 proteins.
Molecules 26 00674 g003
Figure 4. Binding pose of glycyrrhizin (in yellow stick) at the substrate binding pocket of Mpro (in the green surface). A 2D illustration of the glycyrrhizin chemical interactions at the docked site is also provided.
Figure 4. Binding pose of glycyrrhizin (in yellow stick) at the substrate binding pocket of Mpro (in the green surface). A 2D illustration of the glycyrrhizin chemical interactions at the docked site is also provided.
Molecules 26 00674 g004
Figure 5. Binding pose of glycyrrhizin (in yellow stick) at the junction binding pocket of PLpro (in blue surface). A 2D illustration of the glycyrrhizin chemical interactions at the docked site is also provided.
Figure 5. Binding pose of glycyrrhizin (in yellow stick) at the junction binding pocket of PLpro (in blue surface). A 2D illustration of the glycyrrhizin chemical interactions at the docked site is also provided.
Molecules 26 00674 g005
Figure 6. Binding pose of glycyrrhizin (in yellow stick) at the junction binding pocket of PLpro (in deep sky blue surface). A 2D illustration of the glycyrrhizin chemical interactions at the docked site is also provided.
Figure 6. Binding pose of glycyrrhizin (in yellow stick) at the junction binding pocket of PLpro (in deep sky blue surface). A 2D illustration of the glycyrrhizin chemical interactions at the docked site is also provided.
Molecules 26 00674 g006
Figure 7. Statistical parameters calculated based on simulation trajectories. Receptor RMSD plots (A), glycyrrhizin RMSD plots (B), receptor RMSF (C), and receptor Rg (D).
Figure 7. Statistical parameters calculated based on simulation trajectories. Receptor RMSD plots (A), glycyrrhizin RMSD plots (B), receptor RMSF (C), and receptor Rg (D).
Molecules 26 00674 g007
Figure 8. Binding mode dynamics of glycyrrhizin in the MD simulation at the initial time (in tan stick) versus at the end time (in plum stick).
Figure 8. Binding mode dynamics of glycyrrhizin in the MD simulation at the initial time (in tan stick) versus at the end time (in plum stick).
Molecules 26 00674 g008
Table 1. Binding free energy components of SARS-CoV-2 enzyme complexes with glycyrrhizin. The energy values are provided in units of kcal/mol.
Table 1. Binding free energy components of SARS-CoV-2 enzyme complexes with glycyrrhizin. The energy values are provided in units of kcal/mol.
MethodEnergy ComponentMpro–GlycyrrhizinPLpro–GlycyrrhizinNucleocapsid–Glycyrrhizin
MMGBSAVan der Waals Energy−36.50−61.10−37.97
Electrostatic Energy−13.92−8.538.75
Polar Solvation Energy30.1926.793.15
Nonpolar Solvation Energy−4.19−5.85−3.97
Gas Phase Energy−50.42−69.63−29.22
Solvation Energy25.9920.93−0.82
Total Binding Energy−24.42−48.69−30.05
MMPBSAVan der Waals Energy−36.50−61.10−37.97
Electrostatic Energy−13.92−8.538.75
Polar Solvation Energy42.5635.776.65
Nonpolar Solvation Energy−2.94−4.31−3.38
Gas Phase Energy−50.42−69.63−29.22
Solvation Energy39.6231.463.27
Total Binding Energy−10.80−38.17−25.95
Table 2. Hotspot residues identified that played a significant role in interaction with the glycyrrhizin.
Table 2. Hotspot residues identified that played a significant role in interaction with the glycyrrhizin.
ComplexResiduesMMGBSAMMPBSA
Mpro–GlycyrrhizinLys137−1.74−1.51
Asp197−1.76−0.45
Thr198−1.50−1.76
Thr199−1.18−2.84
Tyr237−1.46−1.89
Asn238−2.98−3.45
Tyr239−1.54−1.48
Leu271−1.24−1.69
Leu272−1.65−3.48
Gln273−1.14−1.24
Asn274−1.56−1.42
Met276–1.73−1.98
Ser284−1.89−1.51
Leu286−1.98−1.47
Leu287−1.48−2.84
Glu288−1.44−1.84
Asp289−3.74−3.54
Glu290−1.88−5.45
PLpro–GlycyrrhizinPhe69−2.54−3.54
His73−2.11−2.45
Thr74−1.82−1.45
Asp76−1.99−1.68
Phe79−1.47−1.46
Arg82−1.82−1.12
Asn128–4.41−1.39
Tyr154−1.61−1.48
Asn156−1.61−5.24
Phe173−1.11−1.58
Gln174−5.48−3.61
His175−1.94−1.48
Ala176−1.69−1.12
Asn177−1.81−1.62
Leu178−1.64−1.11
Asp179−2.47−1.83
Val202−1.43−1.19
Nucleocapsid–GlycyrrhizinAsn49−1.99−2.54
Thr50−1.81−1.42
Ala51−1.25−2.45
Phe54−1.66−3.15
Thr55−1.65−1.12
Arg89−2.74−1.84
Thr92−2.45−3.65
Arg94−4.66−2.48
Tyr112−1.45−1.24
Tyr110−3.74−3.51
Pro118−1.89−1.48
Arg150−2.78−3.58
Table 3. WaterSwap absolute binding energy estimation for all four complexes.
Table 3. WaterSwap absolute binding energy estimation for all four complexes.
Algorithm Mpro–GlycyrrhizinPLpro–GlycyrrhizinNucleocapsid–Glycyrrhizin
Bennett’s−22.39−25.84−22.34
Free energy perturbation−22.48−25.94−23.83
Thermodynamic integration −22.47−24.61−23.45
Mean −22.44−25.46−23.30
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Muhseen, Z.T.; Hameed, A.R.; Al-Hasani, H.M.H.; Ahmad, S.; Li, G. Computational Determination of Potential Multiprotein Targeting Natural Compounds for Rational Drug Design Against SARS-COV-2. Molecules 2021, 26, 674. https://doi.org/10.3390/molecules26030674

AMA Style

Muhseen ZT, Hameed AR, Al-Hasani HMH, Ahmad S, Li G. Computational Determination of Potential Multiprotein Targeting Natural Compounds for Rational Drug Design Against SARS-COV-2. Molecules. 2021; 26(3):674. https://doi.org/10.3390/molecules26030674

Chicago/Turabian Style

Muhseen, Ziyad Tariq, Alaa R. Hameed, Halah M. H. Al-Hasani, Sajjad Ahmad, and Guanglin Li. 2021. "Computational Determination of Potential Multiprotein Targeting Natural Compounds for Rational Drug Design Against SARS-COV-2" Molecules 26, no. 3: 674. https://doi.org/10.3390/molecules26030674

Article Metrics

Back to TopTop