Molecular screening of glycyrrhizin-based inhibitors against ACE2 host receptor of SARS-CoV-2.

The interaction between SARS-CoV-2 Spike protein and angiotensin-converting enzyme 2 (ACE2) is essential to viral attachment and the subsequent fusion process. Interfering with this event represents an attractive avenue for the development of therapeutics and vaccine development. Here, a hybrid approach of ligand- and structure-based virtual screening techniques were employed to disclose similar analogues of a reported antiviral phytochemical, glycyrrhizin, targeting the blockade of ACE2 interaction with the SARS-CoV-2 Spike. A ligand-based similarity search using a stringent cut-off revealed 40 FDA-approved compounds in DrugBank. These filtered hits were screened against ACE2 using a blind docking approach to determine the natural binding tendency of the compounds with ACE2. Three compounds, deslanoside, digitoxin, and digoxin, were reported to show strong binding with ACE2. These compounds bind at the H1-H2 binding pocket, in a manner similar to that of glycyrrhizin which was used as a control. To achieve consistency in the docking results, docking calculations were performed via two sets of docking software that predicted binding energy as ACE2-Deslanoside (AutoDock, -10.3 kcal/mol and DockThor, -9.53 kcal/mol), ACE2-Digitoxin (AutoDock, -10.6 kcal/mol and DockThor, -8.84 kcal/mol), and ACE2-Digoxin (AutoDock, -10.6 kcal/mol and DockThor, -8.81 kcal/mol). The docking results were validated by running molecular simulations in aqueous solution that demonstrated the stability of ACE2 with no major conformational changes in the ligand original binding mode (~ 2 Å average RMSD). Binding interactions remained quite stable with an increased potential for getting stronger as the simulation proceeded. MMGB/PBSA binding free energies were also estimated and these supported the high stability of the complexes compared to the control (~ -50 kcal/mol net MMGB/PBSA binding energy versus ~ -30 kcal/mol). Collectively, the data demonstrated that the compounds shortlisted in this study might be subjected to experimental evaluation to uncover their real blockade capacity of SARS-CoV-2 host ACE2 receptor.

Introduction

In December 2019, a novel severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has emerged in Wuhan Hubei province of China [1-5]. It has been identified as a previously unknown SARS-CoV that has crossed the species barrier to be known as SARS-CoV-2 and causes an atypical pneumonia [6]. This disease is now widespread and is known as coronavirus disease 19 or COVID-19. COVID-19 has turned into a pandemic and has spread across 210 countries across the globe, despite the efforts of health officials in China to limit its spread. It was declared as a pandemic by the WHO on 11th March 2020 [5]. The number of new cases and deaths are increasing with each passing moment. Betacoronavirus is recognized as the family and genus of SARS-CoV-2 which is known for highly pathogenic emerging viruses [7, 8]. The members of this family have been known for their lethality in both birds and mammals. Research is still being carried out on the two very famous viruses of this family that have caused hundreds of deaths over a period of 2 decades: severe acute respiratory syndrome coronavirus (SARS-CoV) and Middle East respiratory syndrome coronavirus (MERS-CoV) [9]. However, to date, scientists have been unable to find an effective treatment against this family of viruses. Currently, lots of therapeutics are under examination including: remdesivir, lopinavir with ritonavir, lopinavir with ritonavir plus interferon beta-1a, and chloroquine or hydroxychloroquine [10]. Therefore, this area is very open for research and study for the establishment of successful treatment methods to limit the spread of COVID-19.

For pursuing any study on SARS-COV-2, it is necessary to understand the structural units of a coronavirus. The basic structure of coronavirus includes an envelope [11, 12], membrane [13], nucleocapsid [14], and Spike proteins [15-18]. These are actually non-segmented positive sense RNA viruses [19]. Spike protein plays a vital role in attachment, fusion, and entry of virus into the host machinery, and hence, it is a main determinant of CoV tropism [17, 20–22]. It forms homotrimers that protrude from the surface of the virus and supports virus entry. These protrusions form spikes on the surface of virus giving it a crown-like appearance. In the initial phase of infection by the virus, angiotensin-converting enzyme 2 (ACE2) mediates the entry by fusion of viral and cellular membranes [17, 23]. This facilitates the binding of the Spike protein to the membrane receptor on the host cells. In case of coronavirus SARS-CoV-2 infection, human ACE2 is exploited for host cell entry [6, 24]. A recent study suggests another route of viral attack through CD147/S-protein [25]. Spike protein has two sub-units: S1 and S2 subunits. S1 binds to ACE2 through the receptor-binding domain (RBD); in this way, it mediates viral attachment to the cellular surface and hence entry into the host [26]. Viral entry also needs priming of Spike protein with serine protease TMPRSS2 in the host cell. This is followed by cleavage of Spike protein at sites S1/S2 and S2, which allow fusion of viral and cellular membranes. The S2 subunit drives the process of fusion [17]. Specific to SARS-CoV-2, a furin cleavage site is present at subunits S1 and S2 boundaries. The unique feature of SARS-CoV-2 from other coronaviruses is processing of the furin cleavage site during biosynthesis [27].

Human ACE2 is a type 1 membrane protein expressed in the heart, lung, intestine, and kidneys [28, 29]. It is the main target for viral entry by coronaviruses. The structure of ACE2 contains an N-terminal peptidase domain, a ~40-residue intracellular segment and a C-terminal collectrin-like domain that ends with a single transmembrane helix [30, 31]. A zinc-dependent PD domain is its single transmembrane helix which makes it daunting to establish the full length protein structure. This is the reason for limited availability for structural information regarding ACE2. Physiologically, it is a part of the renin-angiotensin system (RAS) and catalyzes cleavage of angiotensin II to angiotensin 1–7. The latter is a potent vasodilator in the cardiovascular system with defensive properties [32].

Recently, determination of a cryo-electron microscopy structure at a resolution of 3.5 Å of SARS-CoV-2 Spike trimer in prefusion conformation has been performed to facilitate countermeasure developments [21]. It was revealed that in this structure, one of the three RBDs rotated up to facilitate availability to the receptor. Biophysical and structural evidence is provided that shows binding of the SARS-CoV-2 protein ectodomain with PD of human ACE2 with an affinity ~10 to 20 fold higher than that shown in binding with SARS-CoV [21, 33]. The recognition of the viral RBD is mainly performed through polar residues of extracellular PD of ACE2. It is revealed in molecular studies that the two S protein trimers bind to an ACE2 homodimer simultaneously [33]. The residues in SARS-CoV-2 RBD critical for binding with ACE2 are also validated through X-ray crystallographic data. The majority of these residues are highly conserved and the side chain properties are similar for both SARS-CoV-2 and SARS-CoV RBDs [34]. Largely, the COVID-19 research is being focused on production of vaccines [35] and targeting of nonstructural proteins: main protease (Mpro) and RNA-dependent RNA polymerase (RdRp) [6]. However, our focus is on investigation of novel infection prevention strategies against SARS-CoV-2. In 2004, crystal structures of free and bound forms of ACE2 were first solved at the resolutions of 2.2 Å and 3.0 Å, respectively [36]. Huentelman et al. discovered a novel inhibitor against ACE2 using structure-based virtual screening approach for a chemical library of ~140,000 diverse synthetic compounds [37]. The novel compound proved effective in blocking the S protein-mediated SARS-CoV cell fusion. The findings of this study support drug repositioning as a concept. Drug repositioning is also known as drug repurposing and it describes the process of creating novel clinical opportunities for available drugs [38]. This process, in comparison to de novo drug discovery, truncates the chemical structure optimization time and reduces the toxicological testing time. Previously for the treatment of other coronavirus infections, including MERS-CoV, this strategy has also been used [39]. In the current study, employing a molecular mechanics (MM)-assisted structure-based virtual screening method, an investigation of glycyrrhizin-based inhibitors to act as potential inhibitors of ACE2 has been undertaken. Glycyrrhizin is a plant-based derivative and has been reported to show anti-HVC, anti-HIV, and anti-ACE activity [40, 41]. This compound has been demonstrated to effectively neutralize SARS-CoV-2 in vitro and bind to the ACE2 receptor [40, 42–46]. New drug molecules are uncovered in this work, which are structurally related to glycyrrhizin and have been already approved by the Food and Drug Administration (FDA). Thus, the screened drugs can be repurposed against COVID-19 pandemic disease.

Materials and methods

Figure 1 demonstrates the methodology of this work.

Fig. 1

Schematic representation of the flow of steps used in the present study

Retrieval of ACE2 receptor structure and processing

The 3D structure of the ACE2 receptor was retrieved from RCSB Protein Databank (PDB ID: 6LZG) [34]. The initial PDB file was treated in a preparatory phase where all non-relevant co-crystalized structures and ligands were discarded in UCSF Chimera, alpha version 1.15 [47]. A short energy minimization run was then performed to add missing atoms if left in the crystallization process and remove steric clashes. During energy minimization, the first hydrogen atoms were added to the protein and missing bonds were repaired. The structure was then subject to 100 steps of energy minimization via steepest descent and 10 steps of conjugate gradient algorithms while keeping the step size to 0.02 Å.

Ligand based similarity search

Glycyrrhizin, which is a phytochemical and has been extensively documented for its potential for antiviral activity as well for its ability to block SARS-CoV ACE2 receptor, was used as the parent compound in the ligand-based virtual screening [40, 41, 48]. The similarity search was done in DrugBank 5.0 [49] employing a very stringent threshold of 0.6 (60%) and FDA-approved drugs. All the structures of the potent compounds were extracted from DrugBank in MDL Mol format and converted into PDB format using Open Babel [50]. Pre-optimization was performed with Chem3D Pro software [51] employing an MM2 force field [52]. The compounds were then uploaded to PyRx 0.8 [53] where they were converted to pdbqt. Meanwhile, the ACE2 receptor was also added to the PyRx AutoDock panel in preparation for structure-based virtual screening [54].

Molecular docking

For molecular docking and scoring procedures, PyRx AutoDock 4.2.6 was used [55]. To allow free binding of the compounds with ACE2, a blind docking strategy was applied where the search space center at XYZ coordinates was set as −25.16 Å, 15.93 Å, and − 25.12 Å, respectively. The grid dimensions were as follow: X (62.34 Å), Y (69.55 Å), and Z (71.84 Å). The number of iterations set for each compound was 100. Also, the criterion for RMSD values and binding energy values in docking experiments was around 2 Å and below −7 kcal/mol, respectively. The selected compounds were subjected to PreADMET toxicity predictor for Ames mutagenicity and rat carcinogenicity evaluation [56].

Molecular dynamics simulations

Molecular dynamic studies were carried out to evaluate the binding force of compounds in the active site, as shown by the docking studies. The role of water molecules and conformational dynamics of the ACE2 were also assessed. A total of 50 ns full-scale atomic simulations were undertaken using AMBER18 software [57] with a general atom force field (GAFF) [58] for small molecules and ff14SB [59] for ACE2. The antechamber in AMBER [60] was used to add atomic charges, together with a force field for the purpose of generating parameters of compounds for non-bonded interaction with enzymes using the AM1-BCC method [61]. All the atoms and bonds were incorporated according to their topologies after analysis of the structures. Leap (tleap) was used to build topologies and initial coordinates of complexes to run the molecular dynamics simulations [62]. Simulations were performed under periodic boundary conditions in a rectangular size cell with dimensions 65 × 65 × 75 Å3 with a 10 Å non-bonded cut-off and using Particle Mesh Ewald (PME) for the long-range electrostatics [63]. Initially, ACE2-compound complexes were subjected to energy minimization using a total of 1000 minimization steps with 500 steepest descent and 500 conjugate gradient steps to 0 K. In this way, temperature was raised to 100 K in the first cycle. Later, the system was equilibrated to 300 K in the next cycle. Finally, the production run was carried out at a constant temperature and pressure for 50 ns each with an integration step size of 2 fs using SHAKE [64]. Simulations were run at a pH value of 7 and the protonation state of all residues was kept fixed. This was done according to the fact that the virus interactions with the ACE2 are activated at neutral pH [65]. System stability was examined by calculating root mean square deviation (RMSD) for both ACE2 and compounds of all four complexes using AMBER CPPPTRAJ [66]. Additionally, hydrogen bond analysis for all frames of the simulation trajectories was completed using VMD [67].

Binding free energy calculation (MMGB/PBSA)

In this study, MMPBSA implemented in MMPBSA.py script was used for the binding free energy calculations [68, 69]. For this purpose, complexes were selected on the basis of best docking scores. In the first step, multiple snapshots (100) are generated from a molecular dynamics simulation production run trajectory for the respective complex. The final step is to calculate free energy for each snapshot and the formula for estimation of free energy is stated as:

Effective binding free energy is considered an entropic term and is ignored in the binding free energy calculation. In this way, energy calculations were performed for each complex.

Results and discussion

Glycyrrhizin, a triterpene glycoside, is a promising and safe antiviral natural product extracted from the roots of Glycyrrhiza glabra [40, 41, 48, 70]. This drug has a long established history of medical use and has been in use for treating atopic dermatitis, viral hepatitis, influenza A virus, Epstein Barr virus, varicella-zoster virus, parainfluenza virus type 2, dengue virus, chikungunya virus, Semliki Forest viruses, human immunodeficiency virus (HIV), and SARS-CoV and many more make this drug as a good candidate either alone or in combination to be evaluated against SARS-CoV-2 [71]. Considering this, we performed a ligand-based similarity search that filtered 40 compounds as analogs of glycyrrhizin and used in screening against SARS-CoV-2 ACE2 receptor (Table 1).

Table 1

Glycyrrhizin based similarity screen analogs in DrugBank

Name	Accession number	Groups	Weight	Chemical formula	Score
Ciclesonide	DB01410	Approved, investigational	540.697	C32H44O7	0.741
Budesonide	DB01222	Approved	430.5339	C25H34O6	0.736
Desonide	DB01260	Approved, investigational	416.5073	C24H32O6	0.707
Steviolbioside	DB12434	Approved, investigational	642.739	C32H50O13	0.693
Amcinonide	DB00288	Approved	502.5717	C28H35FO7	0.67
Spinosad	DB08823	Approved, investigational, vet approved	1477.963	C83H132N2O20	0.674
Hydrocortisone cypionate	DB14541	Approved, investigational, vet approved	486.649	C29H42O6	0.669
Prednicarbate	DB01130	Approved, investigational	488.577	C27H36O8	0.669
Deslanoside	DB01078	Approved	943.0791	C47H74O19	0.667
Hydrocortisone valerate	DB14544	Approved, vet approved	446.5763	C26H38O6	0.667
Flurandrenolide	DB00846	Approved	436.5136	C24H33FO6	0.663
Loteprednol etabonate	DB14596	Approved	466.96	C24H31ClO7	0.663
Hydrocortisone butyrate	DB14540	Approved, vet approved	432.557	C25H36O6	0.66
Hydrocortisone probutate	DB14543	Approved, vet approved	488.613	C28H40O7	0.66
Hydrocortisone succinate	DB14545	Approved	462.539	C25H34O8	0.66
Ouabain	DB01092	Approved	584.6525	C29H44O12	0.655
Clascoterone	DB12499	Approved, investigational	402.531	C24H34O5	0.649
Digoxin	DB00390	Approved	780.9385	C41H64O14	0.648
Digitoxin	DB01396	Approved, investigational	764.9391	C41H64O13	0.648
Hydroxyprogesterone caproate	DB06789	Approved, investigational	428.6041	C27H40O4	0.647
Gestonorone caproate	DB14677	Approved	414.586	C26H38O4	0.647
Acetyldigitoxin	DB00511	Approved	806.9757	C43H66O14	0.647
Flunisolide	DB00180	Approved, investigational	434.4977	C24H31FO6	0.638
Hydrocortisone acetate	DB14539	Approved, vet approved	404.4966	C23H32O6	0.636
Ivermectin	DB00602	Approved, investigational, vet approved	1736.1589	C95H146O28	0.628
Cortisone acetate	DB01380	Approved, investigational	402.4807	C23H30O6	0.628
Prednisolone tebutate	DB14632	Approved, vet approved	458.587	C27H38O6	0.625
Methylprednisolone hemisuccinate	DB14644	Approved	474.55	C26H34O8	0.623
Methylprednisolone aceponate	DB14643	Approved, vet approved	472.578	C27H36O7	0.617
Medroxyprogesterone acetate	DB00603	Approved, investigational	386.5244	C24H34O4	0.614
Halcinonide	DB06786	Approved, investigational, withdrawn	454.96		0.611
Nandrolone decanoate	DB08804	Approved, illicit	428.6472	C28H44O3	0.607
Testosterone cypionate	DB13943	Approved	412.614	C27H40O3	0.607
Testosterone enanthate	DB13944	Approved	400.594	C26H40O3	0.607
Testosterone undecanoate	DB13946	Approved, investigational	456.711	C30H48O3	0.607
Fluocinonide	DB01047	Approved, investigational	494.5249	C26H32F2O7	0.606
Fluocinolone acetonide	DB00591	Approved, investigational, vet approved	452.4882	C24H30F2O6	0.605
Clarithromycin	DB01211	Approved	747.9534	C38H69NO13	0.602
Drospirenone	DB01395	Approved	366.4932	C24H30O3	0.6
Prednisolone acetate	DB15566	Approved, vet approved	402.4807	C23H30O6	0.6
Glycyrrhizic acid	DB13751	Approved, experimental	822.942	C42H62O16	Control

Molecular docking studies

Molecular docking is crucial in forecasting the predominant binding mode of a ligand molecule to a biological macromolecule and determining binding affinity of the complex [72, 73]. Traditionally, this process is tedious, long, and costly. The introduction of computer-aided drug design undoubtedly played a significant contribution in overcoming the many limitations of conventional drug design and rationalizing the drug discovery process [74]. Molecular docking was swiftly found to have value in the modern structure-based drug design in predicting the binding mode and intermolecular framework of chemical interactions between small drug molecules and proteins that drive inhibition of protein functionality [75]. The reported antiviral agent glycyrrhizin and a library of its closed 40 analogues were screened against the ACE2 receptor. To ensure consistency in the docking process, we initially tested the docking protocol by running a blind docking calculation of the co-crystalized ligand with ACE2 receptor. An identical binding mode was predicted, thus affirming the accuracy of the docking process. Further, consistency in the compounds’ binding mode was examined by performing docking the compounds with ACE2 via two docking platforms. It was also ensured that selected compounds strongly bind to the ACE2 residues that are in direct contact with Spike protein [34]. The compounds used in the screening process were ranked based on the binding energy, and the results were compared to opt top 3 molecules that have the least binding energy in kcal/mol. The best binding complexes are ACE2-Deslanoside (AutoDock binding energy = −10.3 kcal/mol and DockThor = −9.53 kcal/mol), ACE2-Digitoxin (AutoDock binding energy = −10.6 kcal/mol and DockThor = −8.84 kcal/mol), and ACE2-Digoxin (AutoDock binding energy = −10.6 kcal/mol and DockThor = −8.81 kcal/mol). The control ACE2-Glycyrrhizin complex has AutoDock binding energy = −10.3 kcal/mol and DockThor binding energy of −9.53 kcal/mol. A general description of the binding energy of all 41 compounds used in both of the sets of docking software is provided in Fig.2. The obtained docking result demonstrated that the minimum binding energy in AutoDock is −11.2 kcal/mol, with a maximum of −6.1 kcal/mol and mean of −8.42 kcal/mol. In case of DockThor, the minimum, maximum, and mean binding energy values are −9.64 kcal/mol, −6.8 kcal/mol, and − 7.91 kcal/mol, respectively. Moreover, for the ease of the experimentalists, the dissociation constant (Ki) value was predicted for the compounds: deslanoside (2.13 μM), digitoxin (1.48 μM), digoxin (1.17 μM), and glycyrrhizin (2.87 μM). The compounds were further screened for toxicity moieties and predicted to be negative for the Ames mutagenicity assay, and rat carcinogenicity.

Fig. 2

Distribution of binding energy of compounds in both AutoDock and DockThor

Deslanoside-ACE2 complex

Deslanoside (4-(12,14-dihydroxy-3-((4-hydroxy-5-((4-hydroxy-5-((4-hydroxy-6-methyl-5-((3,4,5-trihydroxy-6-(hydroxymethyl)tetrahydro-2H-pyran-2-yl)oxy)tetrahydro-2H-pyran-2-yl)oxy)-6-methyltetrahydro-2H-pyran-2-yl)oxy)-6-methyltetrahydro-2H-pyran-2-yl)oxy)-10,13-dimethylhexadecahydro-1H-cyclopenta[a]phenanthren-17-yl)furan-2(5H)-one) is docked along helical 1 and 2 RBD binding sites at the ACE2 receptor. The binding conformation is towards the H3 helical region where the chemical moiety 4-(3a,7-dihydroxy-7a-methyloctahydro-1H-inden-1-yl)furan-2(5H)-one interacts with the H3 helical region through balanced hydrogen and hydrophobic contacts. The 6-(hydroxymethyl)tetrahydro-2H-pyran-2,3,4,5-tetraol of the compound interacts through rich weak hydrophobic binding. The central backbone chemical moiety of 6-((4-hydroxy-2-methyl-6-((4a-methyldecahydronaphthalen-2-yl)oxy)tetrahydro-2H-pyran-3-yl)oxy)-2-methyltetrahydro-2H-pyran-3,4-diol exhibited more affinity towards the H1 helical site compared to terminal rings that favor binding H2 helix. The docked conformation of the deslanoside and ACE2 and its chemical interactions are illustrated in Fig.3.

Fig. 3

3D and 2D illustration of deslanoside binding mode and chemical interactions at the SARS-CoV-2 RBD binding site of ACE2 receptor. A Docked binding mode of Deslanoside (yellow color stick) respect to surface ACE2 (the most hydrophobic are tan whereas polar residues are colored are medium purple). B Different molecular interactions of deslanoside (shown in ball and stick) involved in the binding of deslanoside with ACE2 residues

Digitoxin-ACE2 complex

Like deslanoside, digitoxin (4-(3-((5-((5-((4,5-dihydroxy-6-methyltetrahydro-2H-pyran-2-yl)oxy)-4-hydroxy-6-methyltetrahydro-2H-pyran-2-yl)oxy)-4-hydroxy-6-methyltetrahydro-2H-pyran-2-yl)oxy)-14-hydroxy-10,13-dimethylhexadecahydro-1H-cyclopenta[a]phenanthren-17-yl)furan-2(5H)-one) is also aligned along the long cylindrical pocket of H1 and H2 binding pocket. The compound binding is dominated by hydrophobic interactions along the length of interactions with the pocket residues. The major center of the compound’s structure was elevated whereas the terminal moieties were dipped. The 4-(3a-hydroxy-7a-methyloctahydro-1H-inden-1-yl)furan-2(5H)-one moiety of the compound was configured to the H3 helix and adjusted near to the H1 helix. The opposite 6-methyltetrahydro-2H-pyran-2,4,5-triol part of the compound lay perpendicular on the cavity center to facilitate interactions with key residues. The docked conformation of digitoxin and ACE2 and its chemical interactions are illustrated in Fig.4.

Fig. 4

Binding mode and chemical interactions of SARS-CoV-2-ACE2 and digitoxin. A Docked binding mode of digitoxin (yellow color stick) and ACE2 (hydrophobicity surface). B van der Waals interactions of digitoxin (shown in ball and stick) involved in the binding with ACE2 residues

Digoxin-ACE2 complex

Among the hits, digoxin (4-(3-((5-((5-((4,5-dihydroxy-6-methyltetrahydro-2H-pyran-2-yl)oxy)-4-hydroxy-6-methyltetrahydro-2H-pyran-2-yl)oxy)-4-hydroxy-6-methyltetrahydro-2H-pyran-2-yl)oxy)-12,14-dihydroxy-10,13-dimethylhexadecahydro-1H-cyclopenta[a]phenanthren-17-yl)furan-2(5H)-one) is positioned exactly at the center of the H1-H2 pocket. The bulk of the 4-(12,14-dihydroxy-10,13-dimethylhexadecahydro-1H-cyclopenta[a]phenanthren-17-yl)furan-2(5H)-one of the compound moiety associates with the H2 helix where its furan-2(5H)-one ring is pointed to the H helix. The region 6-((4-hydroxy-6-((4-hydroxy-2-methyltetrahydro-2H-pyran-3-yl)oxy)-2-methyltetrahydro-2H-pyran-3-yl)oxy)-2-methyltetrahydro-2H-pyran-3,4-diol is sandwiched between the H1 and H2 helix. Here, the 6-methyltetrahydro-2H-pyran-2,4,5-triol ring formed strong hydrogen bonding with Gln42 and Lys68 residues of the pocket. The docked conformation of digoxin and ACE2 and its chemical interactions are illustrated in Fig.5.

Fig. 5

Binding mode and chemical interactions of SARS-CoV-2-ACE2 and digoxin. A Docked binding mode of digoxin (deep pink color stick) and ACE2 (dark cyan hydrophobicity surface). B Chemical interactions of digoxin (shown in ball and stick) involved in the binding with ACE2 residues

Glycyrrhizin-ACE2 complex

The control glycyrrhizin (6-((6-carboxy-2-((11-carboxy-4,4,6a,6b,8a,11,14b-heptamethyl-14-oxo-1,2,3,4,4a,5,6,6a,6b,7,8,8a,9,10,11,12,12a,14,14a,14b-icosahydropicen-3-yl)oxy)-4,5-dihydroxytetrahydro-2H-pyran-3-yl)oxy)-3,4,5-trihydroxytetrahydro-2H-pyran-2-carboxylic acid), as demonstrated before, showed great binding affinity for the ACE2 receptor and showed multiple strong hydrogen and hydrophobic interactions. In particular, 6-((6-carboxy-2,4,5-trihydroxytetrahydro-2H-pyran-3-yl)oxy)-3,4,5-trihydroxytetrahydro-2H-pyran-2-carboxylic acid that contributed highly to the strong interactions anchored the compound at the H1-H2 pocket base and was supported on both helices. This allowed the rest of 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 to project parallel towards the H3 helix. The docked conformation of glycyrrhizin and ACE2 and its chemical interactions are illustrated in Fig.6.

Fig. 6

Binding conformation and chemical interactions of SARS-CoV-2-ACE2 and glycyrrhizin. A Docked binding mode of Glycyrrhizin (dodger blue color stick) and ACE2 (hydrophobicity chartreuse surface). B Chemical interactions of glycyrrhizin (shown in ball and stick) involved in the binding with ACE2 residues

Molecular dynamics simulation

Molecular dynamics simulation is a generally used technique to evaluate the stability and dynamic behavior of a biomolecule or biomolecule/ligand complex under different sets of conditions [76]. An affirmation was done on the initial conformation of the hit compounds-ACE2 complexes obtained after docking in water solution for 50 ns by all atom MD simulations. Uniformity in complex dynamics was noticed as measured by RMSD (Fig.7A). The RMSD values of the ACE2 alpha-backbone carbon atoms in the complex of the four possible inhibitors were measured with respect to the original docked structure which are Deslanoside-ACE2 complex (maximum, 2.79 Å and mean 1.92 Å), Digitoxin-ACE2 complex (maximum, 2.79 Å and mean, 2.14 Å), Digoxin-ACE2 complex (maximum, 2.40 Å and mean, 1.95 Å), and Glycyrrhizin-ACE2 complex (maximum, 2.40 Å and mean 1.95 Å). As can be seen in Fig.7A, all the systems were well converged and upon visualization of the complexes in different time frames reflected stability of the attached compounds to the ACE2 despite few structural changes to get strength in interactions. The stability of complexes was further validated by performing compound RMSD calculations and these were determined as Deslanoside-ACE2 complex (maximum, 2.95 and mean, 1.28), Digitoxin-ACE2 complex (maximum, 3.89 Å and mean, 2.22 Å), Digoxin-ACE2 complex (maximum, 3.50 Å and mean, 1.99 Å), and Glycyrrhizin-ACE2 complex (maximum, 3.50 Å and mean, 1.99 Å). The ligands from all four systems remained in a stable binding and conformation state, with the exception of digitoxin which revealed binding pose variations, as depicted in Fig.7B. Hydrogen bonds played a key function in stabilizing ligand binding at the docked site. The total number of hydrogen bonds formed by ACE2 with the three candidate compounds is shown in Fig.7C. Throughout the simulation time, all four compounds formed hydrogen bonds up to 4 in number with ACE2 receptor. These bonding parameters suggested that all candidate compounds have the ability to bind the ACE2 efficiently and tightly. Thus, the relative stability of the complexes filtered in this study has been inferred in terms of dynamics.

Fig. 7

Molecular dynamics simulation trajectories analysis. A ACE receptor RMSD versus time, B ligands RMSD versus time, and C number of hydrogen bonds in each simulation frame

Binding free energy calculations

The MMGB/PBSA based binding free energy and molecular interactions between the shortlisted compounds and the ACE2 receptor complex have been calculated. Both methods have shown formation of very stable complex interactions. Table 2 illustrates the description of the obtained binding free energies. The entropy term is omitted since, and in some situations, it can not be measured due to convergence problems. As seen in the table, the term electrostatic energy, being most favorable, dominates the total energy of binding. The complex formation contributes to favorable non-polar interactions as compared to the polar part of the solution-free energy through non-favorable contributions. However, the overall van der Waals contribution is still high. The combined average values for the gas phase binding free energy using both MMGB/PBSA methods were determined to be a significant contribution towards overall complex stabilization. The 100 snapshots from the resulting trajectories were averaged for the calculation of all the energy values. All the molecular interactions were driven by polar and non-polar components in a polar solvent, where non-polar interactions and their contributions to energy were more advantageous. In the hydrophobic regions, non-polar residues show a tendency to hide and dissipate water molecules at the binding site. In order to obtain reliable results from absolute binding free energies, convergence of average values determined by MMGB/PBSA must be achieved. Energy oscillations are mainly due to limited sample size, so longer timescales are needed for snapshot extraction.

Table 2

Binding energy values for filter hits-ACE2 complex. Energy values are in kcal/mol

Energy components	Deslanoside-ACE2 complex		Digitoxin-ACE2 complex		Digoxin-ACE2 complex		Glycyrrhizin-ACE2 complex
	MMGBSA	MMPBSA	MMGBSA	MMPBSA	MMGBSA	MMPBSA	MMGBSA	MMPBSA
van der Waal energy	−77.16	−77.16	−59.81	−59.81	−60.72	−60.72	−44.21	−44.21
Electrostatic	−277.55	−277.55	−278.82	−278.82	−22.05	−22.05	−26.56	−26.56
Polar solvation energy	301.92	304.22	294.80	294.40	39.22	43.22	45.23	46.62
Non-polar solvation energy	−8.51	−6.12	−6.22	−4.65	−6.66	−4.86	−5.78	−4.59
Delta gas phase energy	−354.71	−354.71	−338.63	−338.63	−82.78	−82.78	−70.78	−70.78
Delta solvation energy	293.41	298.09	288.57	289.74	32.56	38.35	39.45	42.03
Total energy	−61.30	−56.62	−50.05	−48.88	−50.22	−44.42	−31.32	−28.75

Per residue energy decomposition

To better understand the atomistic interaction mechanism and binding pattern of the compounds with ACE2, the residue wise free energy contribution, to the overall binding free energy of the complex, was calculated using MMGB/PBSA methods. All these residues are classified hotspot amino acids which have shown a contribution of <− 1 kcal/mol and have played a significant role in stabilizing the docked complex (Table 3). The residues that together make high free energy contributions <− 1 kcal/mol in both calculations (PB/GB) are Thr9, Phe10, Lys13, Phe14, Glu17, Ala18, Asp20, Leu21, Gln24, Ala47, Lys50, Phe54, Glu57, Gln58, and Asn61. Therefore, the residues greatly add to the total free energy of the studied system. All these residues lie either within or in close proximity of the binding pocket of ACE2 and display a considerably low energy value of in the binding pocket, highlighting their importance. Aforementioned residues are involved in hydrogen and hydrophobic bond interactions with the compounds resulting in a network of chemical interactions around it. The absence of any residue among the residues listed above might ultimately result in the loss of catalytic activity. It was shown that all the hotspot residues were present in or near the protein’s active site. Additionally, the binding free energy was classified into its components such as the electrostatic, van der Waals, and solvation contribution to define the detailed contribution of each key residue. As demonstrated earlier, the total binding free energy is dominated by electrostatic term contributions, and supports the previous results of the energy calculation. The involvement of some of the binding residues in hydrogen bond interactions was not found but these residues showed favorable van der Waals and non-polar energies attributed to their van der Waals contacts with the compounds. Significant polar contributions of ≤ − 1 kcal/mol were observed by most of the active site residues while combining the results of both free energy methods. Phe10, Leu21, Phe54, and Asn61 were found to make considerable contributions in MMGB/PBASA in binding with all compounds.

Table 3

Hotspot residues identified to contribute significantly in interactions with the compound and supporting complex stability. The residues are ranked based on highest energy contribution and energy values are given in kcal/mol

Deslanoside-ACE2 complex		Digitoxin-ACE2 complex		Digoxin-ACE2 complex		Glycyrrhizin-ACE2 complex
MMGBSA	MMPBSA	MMGBSA	MMPBSA	MMGBSA	MMPBSA	MMGBSA	MMPBSA
Leu21 (−3.2)	Phe54 (−2.7)	Leu21 (−3.5)	Leu21 (−3.1)	Phe10 (−3.1)	Phe10 (2.8)	Lys50 (−2.4)	Phe54 (−2.3)
Phe54 (−3.1)	Leu21 (−2.4)	Phe10 (−3.2)	Phe10 (−2.9)	Lys13 (−3.0)	Leu21 (−2.2)	Phe54 (−2.3)	Leu21 (−1.3)
Phe10 (−2.5)	Phe10 (−2.3)	Lys13 (−2.9)	Lys13 (−2.1)	Leu21 (−2.8)	Lys13 (−2.1)	Leu21 (−2.3)	Ala47 (−1.1)
Lys50 (−2.4)	Ala47 (−1.7)	Phe54 (−2.2)	Phe54 (−1.7)	Phe54 (−2.5)	Phe54 (−2.1)	Ala47 (−1.8)	–
Gln24 (−2.0)	Asn61 (−1.5)	Glu17 (−2.0)	Glu17 (−1.5)	Lys50 (−2.3)	Glu17 (−1.2)	Asn43 (−1.2)	–
Ala47 (−1.7)	Lys13 (−1.05)	Phe14 (−1.57)	Phe14 (−1.4)	Gln58 (−1.6)	Lys50 (−1.0)	Gln24 (−1.1)	–
Asn61 (−1.5)	–	Thr9 (−1.2)	–	Glu17 (−1.1)	–	–	–
Gln58 (−1.5)	–	Gln24 (−1.2)	–	Asn61 (−1.1)	–	–	–
Lys13 (−1.4)	–	Asp20 (−1.0)	–	Thr9 (−1.0)	–	–	–
Glu57 (−1.2)	–	Asn61 (−1.0)	–	–	–	–	–
Thr9 (−1.0)	–	Ala18 (−1.0)	–	–	–	–	–

Conclusions

In this work, ligand-based and structure-based virtual screening approaches were explored for identification of FDA-approved potent chemical compounds as ACE2 inhibitors from DrugBank. Initial sorting based on ligand-based similarity screening resulted in formation of a library of 40 compounds. Further, a blind docking technique was considered that led to the selection of three bioactive molecules (deslanoside, digitoxin, and digoxin) after comparative analysis with control (glycyrrhizin). These compounds bind directly to the site where SARS-CoV-2 Spike interacts, so it is highly likely that these compounds will interfere with the Spike-ACE2 interaction thus hindering the binding and fusion event in the life cycle of the virus. Extensive dynamics evaluation of the complexes was also performed in an aqueous environment that elucidated stable conformational stability of the ACE2 and of the compounds. Binding free energy estimation showed the compounds had better binding affinity compared with the control. In the interactions, Phe10, Leu21, Phe54, and Asn61 were found to contribute significantly in stabilizing the compounds at the docked site. Hence, this investigation indicated that the selected compounds could be subjected to further investigation as a novel set of leads for the rational drug designing of ACE2 inhibitors in COVID-19 therapy.