Modulating functional allostery of the host-cell receptor protein hACE2 to inhibit viral entry of SARS-CoV-2

Abstract

The emergence of new SARS-CoV-2 Omicron sub-variants with faster transmission has necessitated accelerated scientific efforts to confront a possible health emergency. Conventional anti-CoV strategies targeting viral proteins often fail due to frequent mutations. Thus, targeting the conserved host receptor angiotensin-converting enzyme 2 (hACE2), which mediates viral entry via interaction with the spike protein's receptor-binding domain (RBD), presents a rational therapeutic alternative. This, however, requires identification of non-orthosteric hACE2 sites and suitable modulators that retain hACE2's physiological function. Using blind docking and unbiased molecular dynamics (MD) simulations, we identify a novel allosteric site on hACE2, distant from its peptidase domain. Simulations show that an allosteric modulator can disrupt hACE2–RBD interaction by perturbing the spike RBD while stabilizing hACE2's binding to its natural substrate, angiotensin II (AngII). Pharmacophore modeling and high-throughput virtual screening (HTVS) of large databases yield more effective modulators. These allosteric binders downregulate hACE2–RBD interaction across three SARS-CoV-2 variants of concern (Beta, Delta, and Omicron). Dynamic residue network analysis reveals the shortest suboptimal pathway through which the allosteric signal is transmitted to the RBD. We believe that the identified site and mechanistic insight offer a promising basis for developing variant-agnostic SARS-CoV-2 therapeutics.

---

1. Introduction

Despite the World Health Organization (WHO)'s official declaration of ending the COVID-19 pandemic,^1,2 scientific research continues with the hope of discovering effective treatment options for possible future emergencies that may be more challenging.³ Severe acute respiratory syndrome-coronavirus-2 (SARS-CoV-2) has posed a serious threat to humanity since its onset around late 2019 in a Chinese province and has resulted in more than 775 million confirmed cases and 7 million deaths worldwide to date.⁴ Promulgated vaccination programs have brought the infection rate under control,⁵ but complete protection against the virus remains questionable.⁶ Moreover, future mutants of the virus or forthcoming emergence of new microbes might lead to even more drastic immune invasion and fatalities.⁷ These serve as motivation to carefully extend related investigations, if not with the earlier urgency. Currently, some new variants are reported to cause a rise in cases, among which JN.1 is estimated to spread most rampantly.⁸

Interaction of the α1 helix of the N-terminal peptidase domain (PD) of hACE2 and the RBD of the SARS-CoV-2 spike protein is crucial for viral entry into the host cell.^9,10 Hence, targeting the host protein may be an unusual but effective strategy to prevent viral entry and transmission.¹¹ Although most drugs developed/discovered focus on inhibiting viral replication by targeting the key protease enzymes of the virus,^12,13 works on interrupting the hACE2–RBD interaction to prevent viral entry are relatively scarce.¹⁴ However, with subsequent mutations in the viral genome, including the spike protein, a specific drug cannot be effective against all variants, necessitating repeated optimization and verification of drugs.¹⁵ In this context, targeting the conserved host receptor protein hACE2 could be a promising and unique approach.¹⁶

hACE2 (EC 3.4.17.23) is a hydrolase classified as a mono-carboxypeptidase, contrary to di-carboxypeptidase ACE, which is found in both membrane-bound (mACE2) and soluble (sACE2) forms in cells.^17,18 The primary enzymatic function of hACE2 is to aid in the positive regulation of endogenous cardiovascular activities as a part of the renin-angiotensin-aldosterone system (RAAS) by hydrolyzing bioactive vasoconstrictor octapeptide angiotensin II (AngII) into angiotensin (1–7).^17,19 Angiotensin (1–7) serves as a vasodilator with anti-proliferative and anti-apoptotic activities;²⁰ it indirectly contributes to vasoprotection by binding to a 7TM GPCR, a Mas receptor that induces a cascade of reactions, activating protein kinase B and nitric oxide production. This is contrary to the AngII-AT1/2R axis, known to cause hypertension; the failure of hACE2-mediated cleavage of AngII is regarded as the possible cause of hypertension in SARS-CoV-2-infected individuals.^21,22 Hence, there is a need to develop therapeutic strategies that modulate hACE2 to weaken its interaction with the viral spike RBD while preserving its native functions, since its active site cannot be targeted due to its critical physiological role.

The phenomenon of synergistic changes in conformational dynamics of two different binding sites on a protein is called allostery, enabled by an effector that binds to a protein region other than its catalytic substrate binding (or orthosteric) site.^23,24 An orthosteric inhibitor binds directly to the active site of the enzyme or receptor, the same site where the natural substrate or ligand would bind. An allosteric inhibitor, by contrast, binds to a distinct site on the target that is different from the active site, called the allosteric site. An allosteric regulator or modulator binding induces a conformational change in the protein's structure, which can either enhance or suppress its activity. Positive allosteric modulators (PAMs) increase the enzyme's activity by stabilizing a conformation that makes the active site more favorable for substrate binding or catalysis. Negative allosteric modulators (NAMs) impair the enzyme's catalytic activity by inducing conformational changes that reduce the active site's affinity for the substrate.^24,25

While counteracting viral pathogens, targeting host receptors is a novel but challenging therapeutic strategy. A great deal of investigation is indispensable for discovering, characterizing, and assessing cryptic non-orthosteric sites to reduce adverse effects on the native function of the host protein. Allosteric inhibition also requires deeper investigation to assess its impact on the molecular mechanism and its effectiveness against pathogenic mutant variants, which remain poorly understood to date. Hence, we employ an in silico approach to identify a novel allosteric site and map how its modulation propagates to disrupt the host-virus interaction. Using microsecond-scale unbiased MD simulations, we demonstrate that binding a potential agonist to this designated pocket not only attenuates the hACE2–RBD assembly but also potentiates the binding of the hACE2 active site with its natural substrate AngII. Fostered by pharmacophore modeling, a few putative allosteric modulators are also identified as effective against the major SARS-CoV-2 variants of concern. Finally, dynamic network analysis of Cα atoms discerns a pathway linking the allosteric site to the RBD interface, offering mechanistic insights into functional allostery.

2. Methods

2.1. Protein and ligand preparation

Different crystal structures of hACE2 bound to the SARS-CoV-2 spike RBD of wild-type (PDB ID: 6M0J),⁹ Beta (PDB ID: 8DLN),²⁶ Delta (PDB ID: 7TEW),²⁷ and Omicron (PDB ID: 7T9L)²⁸ variants were retrieved from RCSB²⁹ for in silico studies. The co-crystallized ligands and water molecules were removed from the crystal structures, and spike RBD and hACE2 structures were saved as separate PDB files. The reference molecule SB27012 was drawn using ChemDraw and was optimized utilizing density functional theory (DFT) at the B3LYP/6-311G (d,p) level using Gaussian 16³⁰ with the GaussView 6.0³¹ interface.

2.2. Molecular docking

The docking protocol adopted herein is the same as that used in our previous work.³² Given the widespread use and reasonable performance of AutoDock Vina in generating plausible binding poses, we employed it for our docking studies.^33,34 The optimized reference molecule SB27012 was first docked into the hACE2 protein after assigning the Gasteiger charges with Vina, using the Lamarckian genetic algorithm (GA) and grid-based energy estimation incorporated in PyRx 0.8.³⁵ The wild-type spike RBD was then docked to this protein–ligand complex using HADDOCK2.4.³⁶

2.3. High-throughput virtual screening

We developed a structure-based pharmacophore by inspecting critical interactions between the reference molecule and hACE2, which were retained throughout the simulation. It was subjected to high-throughput virtual screening (HTVS) along with the hACE2 protein on the Pharmit server³⁷ against 10 available databases containing approximately 0.35 billion small molecules or ligands. HTVS is a popular method for rapidly identifying small molecules that might exhibit properties similar to those of the ref. 38. Hit molecules with a docking score cutoff of −9.5 kcal mol⁻¹ were chosen for further investigation.

2.4. Determination of pharmacokinetic properties

The hit molecules were further filtered out on the basis of their pharmacokinetic and drug metabolism properties, which were computed using two web servers: SwissADME³⁹ and admetSAR 2.0.⁴⁰ Lipinski's rule of five⁴¹ was taken as the primary consideration, along with absorption, distribution, metabolism, excretion, and toxicity (ADMET) profiles, to pick out molecules with better druggability potential, which were then optimized for further use. The missing hydrogens were added to the small molecules and optimized following the same method described previously.

2.5. MD simulations

MD simulations were carried out using GROMACS (Groningen Machine for Chemical Simulations) v2022.⁴² ACPYPE (AnteChamber PYthon Parser interfacE)^43,44 was used to generate the topology and parameter files of the small molecule SB27012, and those for the proteins were generated using the AMBER ff99SB force field.⁴⁵ Following this, the proteins and ligands were complexed into single entities, and the unit cell was solvated using the TIP3P model⁴⁶ into a cube of 10 Å diameter. To neutralize the system, 0.15 mol L⁻¹ of counterions (Na⁺ and Cl⁻) were added to it. Here, we used the AMBER 99SB force field for ions. Energy minimization of the complex was performed by employing the steepest descent integrator⁴⁷ for 5000 steps with a force convergence of <1000 kcal mol⁻¹ nm⁻¹.

Afterward, we equilibrated the complexes for 5 ns employing canonical isothermal–isovolumetric (NVT) and isothermal–isobaric (NPT) ensembles,⁴⁸ during which Berendsen temperature⁴⁹ and Parrinello–Rahman pressure⁵⁰ controllers coupled with the systems were used to maintain a temperature of 300 K and a pressure of 1 bar, respectively. The particle mesh Ewald (PME) algorithm,⁵¹ with a Fourier grid spacing of 0.12 nm, was utilized to compute the long-range Coulomb interaction. Further, Lennard-Jones (LJ) potential⁵² with a cutoff distance of 1 nm was used to model the short-range van der Waals interactions, and all bond lengths were constrained by the linear constraint solver (LINCS) method.⁵³ Finally, 1000 ns production MD simulations were run under the microcanonical ensemble by relaxing the coupling to the thermostats. A time step of 2 fs was used to save the coordinates. The same protocol was used to run MD simulations for all the complexes thereafter. The trajectories were analysed for various post-MD parameters as described in (Section S1.1).

2.6. Binding free energy calculations

Binding free energy is an important metric for estimating the affinity between two molecular partners in their bound state. In computational studies, this is often approximated using endpoint methods such as MM/PBSA (molecular mechanics/Poisson–Boltzmann surface area)⁵⁴ or MMGBSA (molecular mechanics/generalized born surface area).⁵⁵ In this study, we employed the MM/PBSA approach using the g_mmpbsa⁵⁶ tool to evaluate binding energetics. Importantly, we computed binding free energies for two phenomena: (i) small molecule binding with hACE2 and (ii) spike RBD binding with hACE2. For each system, binding free energies were calculated over the first and last 10 ns of the 1000 ns MD trajectory to compare the initial and final interaction strengths. The formulae employed for the estimation of MM/PBSA binding free energy are described in Section S1.2.

2.7. Pull-force calculation

Steered MD (SMD) simulations⁵⁷ were conducted to estimate the mechanical unbinding force between hACE2 and the spike RBD. The starting structure for SMD was selected from the final equilibrated phase of the 1000 ns MD simulation of the hACE2–RBD complex. This frame was chosen based on its structural stability (RMSD plateau) and the preservation of key intermolecular interactions. Prior to the pulling simulation, the system was re-solvated in a cubic water box and neutralized with 0.15 M counterions to ensure compatibility with the pulling protocol. The system then underwent brief energy minimization using the steepest descent algorithm, followed by a short 100 ps NPT equilibration to remove any steric clashes or pressure artifacts introduced during re-solvation. This setup ensures a stable baseline before applying the external pulling force, as is standard in SMD workflows. In the simulation, the RBD was pulled along the z-axis, whereas the hACE2 remained immobile.

2.8. Dynamic network analysis

Dynamic network analysis is an efficient method to study the intramolecular dynamic nature of proteins. It is mainly used to identify nodes (residues) that might be critical for regulating the behavior of another node (residue) of the protein placed distantly from the former.⁵⁸ The network analysis was done using NAPS (network analysis of protein structures).^59,60 The dynamical coordinates and the average structure of the last 10 ns trajectory were used as inputs. The shortest suboptimal path between the chosen residues was determined by visualizing nodes. While the complex network consisting of several C_α nodes and edges is colored in blue, the suboptimal path is shown in maroon. The node centrality calculations were also carried out as a part of the network analysis by the same tool, from which four types of centralities, as mentioned in (Section S1.3), have been considered for the work.

3. Results and discussion

3.1. Discovery of a novel allosteric pocket on hACE2 targeted by SB27012

SB27012 is considered a reference molecule based on an experimental report⁶¹ about its druggable efficacy against the hACE2–spike RBD binding. To elucidate binding modes, we performed blind docking against both the hACE2–spike RBD complex and apo-hACE2 protein, which revealed that the molecule binds to a novel non-orthosteric site of hACE2 (Fig. 1A), far from its peptidase domain (PD). The binding score of this particular mode is −9.5 kcal mol⁻¹. This pocket comprises residues like Leu95, Gln98, Ala99, Gln102, Asn103, Tyr196, Tyr202, Leu392, Lys562, and Trp203 and has an average volume of 640 ų, as determined through DoGSiteScorer.^62,63

Fig. 1 (A) The binding pose of SB27012 in the hACE2–spike RBD complex and a magnified representation of the designated allosteric pocket, (B) binding poses of three SB27012 molecules at 0 and 1000 ns as derived from an unbiased MD simulation, and (C) their RDF and distance analyses with Respect to the designated residues.

To validate the accuracy of the binding mode, determined through global docking, we conducted unbiased MD simulations to observe the dynamic behavior of molecules in a solvent environment. For this purpose, we introduced 10 molecules of SB27012 into a simulation box containing the hACE2–RBD complex and TIP3P water molecules using the gmx insert-molecules tool in GROMACS, which randomly places the ligands in the solvent region while avoiding steric clashes. During the simulation, two molecules, labeled A and B, gradually approach the hACE2 receptor for interaction. Interestingly, another molecule labeled as C moves to the binding interface of the two proteins. The position of these three molecules before and after the simulation is illustrated in Fig. 1B. To assess the propensity of ligand molecules near a binding site, we calculated the radial distribution function (RDF) and the minimum distance of molecules A and B from Tyr196, a residue in the discovered allosteric pocket. Our analysis revealed that molecule A maintains a stable position within the defined pocket with a high-intensity probability peak. Conversely, molecule B exhibits a lower intensity peak and significant deviation from its initial position. We calculated RDF and the minimum distance of molecule C from Glu493, an interface residue of the spike RBD. Fig. 1C shows that molecule C is unstable in its position. This aligns with the experimental data⁶¹ that reported no interaction between SB27012 and the complex at the binding interface.

3.2. Counter-directional repulsion induced by SB27012 between hACE2 and the spike RBD

To evaluate the impact of SB27012 binding to the hACE2 allosteric site on the interaction of hACE2 with the spike RBD and AngII, we examined the hACE2–AngII–spike RBD complex with and without SB27012 through a 1000 ns simulation. Analysis of the superimposed structures at 0 ns and 1000 ns revealed that SB27012 induces a repulsive effect, with hACE2 and the spike RBD moving away from each other (counter-directional), resulting in an increased distance across the protein–protein interface (Fig. S1). Additionally, SB27012 exhibits a significant translational shift from its initial docked pose during the simulation, relocating to a distinct position within the same pocket as indicated in Fig. S1. Alterations in conformational dynamics of the complex in the aqueous environment were quantified using post-molecular dynamics (MD) parameters. Specifically, a radius of gyration (R_g) and root-mean-square-deviation (RMSD) correlation map was generated for the spike RBD of the two systems, and structures from randomly selected frames within their energy-minima basins were superimposed (Fig. 2A). The analysis showed that in the presence of SB27012, the spike RBD exhibits a higher propensity for expansion compared to the complex without SB27012; five helices, H1, H3, H4, H6, and H7, were found to have unfolded into flexible loops. Furthermore, the average R_g and RMSD values of the spike RBD of both systems, along with their respective maximum and minimum propensities (Table 1), indicate reduced compactness and increased atomistic deviation in the SB27012-bound complex, respectively. Additionally, individual R_g and RMSD plots of the spike RBD demonstrate relatively greater fluctuations in the complex with SB27012 (Fig. S2A).

Fig. 2 (A) R_g–RMSD correlation plot presented as a free energy landscape and a superimposed structure of the spike RBD from the two complexes with and without SB27012 at 1000 ns, (B) average SASA and RMSF bar plots of the spike RBD, (C) RDF plots of the hACE2 α1 helix with respect to the two interface β sheets of the spike RBD, (D) the MMPBSA binding free energy bar plot between hACE2 and the spike RBD, (E) pulling force profiles to dissociate the spike RBD from hACE2, and (F) cluster analysis diagrams for both complexes, where each complex is an average structure over every 200 ns taken successively in the order of dark red, cyan, goldenrod, dark green, and magenta.

Table 1 The average, maximum, and minimum scores of RMSD and R_g calculated for the spike RBD from both complexes

Parameters	RMSD (nm)			Rg (nm)
Complex	Average	Highest	Lowest	Average	Highest	Lowest
hACE2–AngII–RBD	0.19942	0.2568	1.1 × 10^−6	2.48279	2.44318	2.58371
hACE2–AngII–SB27012–RBD	0.22559	0.29442	9 × 10^−7	2.49212	2.44477	2.57596

---

Further analysis of the dynamic perturbation of the spike RBD in the complexes was conducted by examining its solvent accessibility and residual behavior. This was illustrated through averaged bar plots depicting solvent-accessible surface area (SASA) and root-mean-square fluctuation (RMSF), as shown in Fig. 2B, respectively. A larger solvent-accessible surface area of the spike RBD in the presence of SB27012 (436.92 nm²) compared to the complex without SB27012 (430.86 nm²) suggests reduced interaction between hACE2 and the spike RBD. Moreover, a higher RMSF value observed for the SB27012-bound complex (0.1255 nm) relative to the complex without SB27012 (0.1234 nm) aligns with previous observations, further supporting the pronounced conformational variations induced by SB27012.

The observed greater perturbations prominent in residual behavior in the former suggest a pronounced influence induced by the binding of the small molecule to the allosteric pocket of hACE2. Specifically, the distal region exhibits higher fluctuations than the proximal interacting surface. Individual plots of SASA and RMSF for the spike protein in both complexes are provided in Fig. S2B. Principal component analysis (PCA) of the spike RBD shows larger eigenvalues and broader subspace coverage in the allosteric ligand-bound complex, indicating enhanced conformational flexibility (Fig. S3A). A PCA-derived eigenvalue quantifies the variance of atomic motion captured along a specific principal component, reflecting the magnitude of collective fluctuations during simulation.

The observed conformational shift of hACE2 and the spike RBD was quantified in terms of distance through RDF calculations of the hACE2 α1 helix relative to the two interface β sheets of the spike RBD (Fig. 2C). The RDF profile shows that the SB27012-bound complex exhibits a peak at 1.594 nm, while the complex without SB27012 shows a peak at 1.276 nm, indicating a greater average separation and reduced proximity between hACE2 and the spike RBD in the presence of the modulator. The MMPBSA binding free energy between hACE2 and the spike RBD was computed using the last 10 ns of the trajectories to quantify the observed decrease in affinity between the two proteins, as inferred from the analyses above. Notably, there is a significant increase in the binding energy between the two complexes, with the SB27012-bound complex exhibiting a value of −334.83 kJ mol⁻¹, while the other complex displays −432.82.16 kJ mol⁻¹. Evidently, low van der Waals energy for the latter (−66.25 kJ mol⁻¹) contributes to the decreased energy, while there are no substantial differences in other contributing energies between the two complexes. The MMPBSA calculation is crucial for assessing the stability of a system. In this case, it indicates that the binding of SB27012 in the allosteric pocket leads to instability of the complex (Fig. 2D). These calculations unanimously demonstrate substantial movements of one or both proteins upon SB27012 binding, lowering their probability of being proximal compared to the absence of the allosteric inhibitor. The force required to pull apart two bonded molecules from each other, i.e., pull force, was employed to validate MMPBSA and RDF results by assessing the interaction between hACE2 and the spike RBD. Interestingly, the pulling simulation reveals that the presence of SB27012 lowers the dissociation force from 772.76 kJ mol⁻¹ (in its absence) to 692.84 kJ mol⁻¹ (Fig. 2E), indicating that binding of SB27012 to the allosteric pocket destabilizes the hACE2–spike RBD complex.

Furthermore, Fig. 2F overlays structures averaged over every 200 ns of the simulation trajectory in different clusters for both systems. Again, a notable trend emerges in the system with SB27012, where the spike RBD displays a greater tendency to dissociate from hACE2 compared to the system without SB27012. This observation aligns with the conclusions drawn from the binding energy calculations. Analysis of the mean smallest distance (MSD) with the residue index of the complex (Fig. S3B) confirms a reduction in the number of contacts between hACE2 and the spike RBD. In addition, we used PDBePISA,⁶⁴ a computational tool for analyzing macromolecular interfaces, which helps estimate properties such as interface area and quantify and interpret changes in protein–protein interactions (PPI) during simulations. The analysis of the last frame of the MD trajectory revealed a reduced interface area in the SB27012-bound complex (814.1 Ų) compared to the unbound system (966.5 Ų), further supporting weakened PPI. Altogether, these findings demonstrate that even in the presence of AngII, SB27012 effectively prevents the interaction between hACE2 and the spike RBD upon binding to the non-orthosteric pocket.

3.3. SB27012 enables better hACE2–AngII interaction

To determine whether SB27012-induced dissociation of hACE2 and the spike RBD maintains hACE2's natural protease activity, we analyzed the structural changes of AngII over time during the MD simulation. Comparing the initial (0 ns) and final (1000 ns) structures reveals a significant alteration in the AngII orientation in the active site pocket of hACE2, presumably triggered by SB27012 binding to the allosteric site (Fig. 3A). However, our focus remains on determining whether this alteration favors the interaction between AngII, the designated natural substrate, and hACE2. To investigate this, we studied the relative conformational modulations through various parameters. Due to the small size of AngII and the comparatively ample hollow space of hACE2 where it binds, there remains great flexibility for the peptide to move, which is why the RMSD profile of AngII in both complexes is seen to be fluctuating throughout the simulation (Fig. S4A). Nevertheless, the RMSD of AngII in both complexes fluctuates throughout the simulation; however, the SB27012-bound system shows a slightly lower average RMSD (0.35 ± 0.05 nm) than the unbound system (0.37 ± 0.04 nm), consistent with a modest stabilization effect.

Fig. 3 (A) Shift in the position of the AngII substrate in the active site pocket of hACE2 in both the complexes, (B) RMSF, (C) SASA of AngII octapeptide residues, (D) MMPBSA binding energy bar plot, and (E) hydrogen bond plot between AngII and hACE2.

RMSF analysis identifies Asp1 and Arg2 as the most fluctuating residues in both complexes. In contrast, others exhibit relatively stable residual fluctuations, though the overall RMSF of the octapeptide is greater in the case of the complex lacking SB27012 (Fig. 3B). SASA calculation for the peptide suggests enhanced compatibility of AngII within the catalytic pocket. Despite a smooth profile, the average SASA value of AngII in the complexes differs by 0.4 nm², with the value for the SB27012-bound complex being lower at 13.87 nm² (Fig. 3C). R_g calculation (Fig. S4B) also aligns with the RMSD and SASA analyses. Lesser solvent accessibility of the substrate in the complex with the allosteric inhibitor implies that the small molecule enables better enzyme-substrate interaction.

The MMPBSA calculation (Fig. 3D) demonstrates that SB27012-less and SB27012-bound complexes furnish −21.02 kJ mol⁻¹ and −38.9 kJ mol⁻¹ binding energies, respectively, which signify robust interaction of the substrate AngII with the enzyme catalytic pocket of hACE2. The low van der Waals and electrostatic energies are countered by very high polar solvation energy in both cases; hence, the latter is inferred to govern the binding of AngII to hACE2. Again, to support these deductions, we performed hydrogen bond analysis between hACE2 and AngII. Surprisingly, we observed a sizable difference in the number of hydrogen bonds between the two complexes (Fig. 3E). The SB27012-bound complex retains considerably more intermolecular H-bonds throughout the simulation. This suggests enhanced substrate binding stability, consistent with the more favorable MMPBSA binding energy and reduced conformational deviation observed. This supports our hypothesis that allosteric modulation by SB27012 preserves or improves hACE2's natural function while disrupting spike RBD binding.

3.4. Pharmacophore modeling and HTVS to yield more potent inhibitors

Through an extensive analysis across multiple time frames spanning 1000 ns, we discovered persistent interactions between SB27012 and the hACE2 allosteric pocket. Leveraging these findings, we constructed a structure-based pharmacophore (Fig. 4) and employed it subsequently for high-throughput virtual screening of vast compound libraries. Over 0.35 billion compounds were screened, identifying five molecules with docking scores better than that of SB27012, which is −9.5 kcal mol⁻¹. The databases from which the compounds were retrieved, along with their identifiers, are given in Table S1. The identified hit compounds contain triazole or tetrazole functional groups, which are commonly used in drug design due to their favorable hydrogen bonding capacity, metabolic stability, and bio-isosteric properties. However, these heterocycles may raise concerns during synthesis, particularly regarding safety and scalability under certain conditions. Therefore, while their presence supports potential efficacy, their incorporation should be carefully evaluated during subsequent lead optimization and experimental development stages.

Fig. 4 Structure-based pharmacophore development from persistent interactions and high-throughput virtual screening leading to the identification of five hit molecules. Their docking scores are given in parentheses.

Pharmacokinetic and ADMET analysis (shown in Table S2) helps identify candidates with favorable drug-likeness, reduced side effects, and optimal therapeutic potential, accelerating the development of safer and more effective treatments.⁶⁵ The identified hit compounds, labeled hit1–5, demonstrate lower carcinogenicity and higher human oral bioavailability than the reference SB27012, as shown in the table. Additionally, they exhibit reduced permeability across the blood–brain barrier (BBB) and possess favorable human intestinal absorption properties. Clearly, all hit molecules adhere to Lipinski's rule of 5 and exhibit acceptable lipophilicity, whereas SB27012 violates two parameters with its excessive rotatable bonds and molecular weight of more than 500 Da. These findings suggest that the hit molecules may serve as more powerful allosteric inhibitors of SARS-CoV-2.

Deliberating the HTVS results, ADMET properties, and docking scores of small molecules from databases, we selected five hit molecules for further investigation. Once these molecules were docked into the allosteric pocket of hACE2, the resulting complexes were docked to the spike RBD and subjected to MD simulations spanning 1000 ns. Detailed 2D interaction analyses for both pre (0 ns) and post (1000 ns) simulation trajectories are presented in Fig. S5. Subsequently, we analyzed various MD parameters across the trajectories of five complexes and generated comparative plots including the SB27012-bound complex. First, we assessed the compatibility and stability of the molecules in the hACE2 allosteric pocket by analyzing the C_α atoms. The hACE2–SB27012 complex yields an average RMSD score of 0.25 nm (Fig. 5A). Among the other complexes with hit molecules, those with hit2, hit4, and hit5 exhibit lower RMSD values of 0.19 nm, 0.23 nm, and 0.21 nm, respectively. Conversely, the complexes with hit1 and hit3 show higher RMSD values compared to the hACE2–SB27012 complex. The time-averaged RMSD plots depict relatively smoother profiles for the hACE2-hit2 and hACE2-hit4 complexes (Fig. S6A).

Fig. 5 (A) The average RMSD bar plots of Cα atoms for hACE2 complexed with the small molecules, (B) plots of eigenvalues versus the corresponding eigenvector indices derived from the hACE2 Cα covariance matrix during the MD simulation, (C) RDF plots of hACE2 Gln42 with respect to the RBD Gly446 residue, (D) the MMPBSA binding free energies between hACE2-small molecules and the spike RBD depicted as bar plots, and (E) mean smallest distance analysis between the residues of the complexes, with the highlighted regions indicating those affected by ligand binding.

The principal component (PC) analysis was conducted on hACE2 C_α atoms to evaluate conformational changes induced by small molecule binding. An essential subspace analysis was performed on the first 15 modes, encompassing >95% of the protein's variance, revealing an exponentially decaying curve of eigenvalues against eigenvectors. While plotting eigenvalues against eigenvectors (Fig. 5B), the hACE2–SB27012 complex renders the highest value at the beginning, whereas hACE2-hit2 shows the lowest, followed by hACE2-hit4/5, hACE2-hit3, and hACE2-hit1. Protein conformation subspaces delineated by the projection of the first two PCs (Fig. S6B) are congruent with the trend of eigenvalues for the complexes. Given that all hit complexes display lower subspace coverage than hACE2–SB27012, the notably lower subspace covered by hit2 and hit4 or hit5 complexes suggests reduced conformational movements and enhanced stability. RMSF and R_g profiles (Fig. S6C and D) further support the RMSD and PCA results by reproducing almost the same trend.

It is pertinent to demonstrate the extent to which the hit molecules perturb the spike RBD, disrupting hACE2–RBD interaction, compared to the reference molecule. The superimposed docking poses of hACE2–spike RBD complexes of hit/reference molecules are shown in Fig. S7. The perturbation in the spike protein was quantified by estimating MD parameters for structural and position changes. Interestingly, the average RMSD, RMSF, R_g, and SASA values (Fig. S8A and D) are all in line with the previous inferences, i.e., the hit2 and hit4 complexes exhibit lower values for these parameters than the SB27012 complex. The rest of the hit molecule complexes exhibit higher values, indicating that they are relatively less stable in the allosteric pocket of hACE2 and cannot exert a significant effect to prevent hACE2–spike RBD interactions. The average values of these parameters are provided in Table S3.

RDF analysis focuses on two key residues, Gln42 of hACE2 and Gly446 of the spike RBD. These residues are consistently observed to engage in a predominant interaction, specifically a hydrogen bond, across all docked protein complexes (Fig. 5C). The analysis reveals minimal displacement in the complexes involving hit1 and hit3, with peak values falling within the permissible range of 0.5 nm, indicative of sustained atomic interactions. On the contrary, other systems, including SB27012, exhibit peaks surpassing this threshold, suggesting no significant interaction between the residues. Notably, the peak displacements are most pronounced in the hit2 and hit4 complexes. In all systems, however, the hydrogen bond initially present between the two residues breaks since the threshold distance for the H-bond is 0.35 nm.⁶⁶ This may result from allosterically induced conformational flexibility, as supported by elevated RMSF and PCA values near the interface, which can perturb the optimal donor–acceptor geometry required for stable hydrogen bonding.

The MMPBSA calculation helps capture any change in the affinity between hACE2 and RBD upon introducing small molecules, as depicted in Fig. 5D and Table S4. The analysis reveals that among the tested compounds, only hit2 (−215.104 kJ mol⁻¹) and hit4 (−153.766 kJ mol⁻¹), when bound to the allosteric site of hACE2, demonstrate superior inhibitory effects on the binding between hACE2 and the spike RBD compared to SB27012 (−351.835 kJ mol⁻¹). Notably, electrostatic energy emerges as the primary contributor to the overall binding free energy, with negligible contributions from van der Waals and SASA energies exhibited by the proteins. The hit5 complex exhibits the lowest binding free energy, representative of minimal alteration in hACE2–spike RBD interaction, attributed markedly to its highly negative electrostatic energy (−578.735 kJ mol⁻¹).

To provide a clear perspective, we computed the MMPBSA binding energy between hACE2 and the ligands again. Hit2 yields the lowest value (−75.23 kJ mol⁻¹), followed by SB27012 (−62.05 kJ mol⁻¹) and hit4 (−61.34 kJ mol⁻¹) with comparable values (Fig. S9A and Table S5). Notably, the van der Waals and polar solvation energy components emerge as the primary contributors to binding. Conversely, hit5 displays the weakest affinity for hACE2 (−23.34 kJ mol⁻¹). The assessment of ligand RMSD (Fig. S9B) corroborates the MMPBSA findings, indicating that hit5 exhibits the highest average deviation, which increases toward the end of the simulation. Additionally, we provide superimposed ligand images to depict positional changes in Fig. S9C. Based on these analyses, it is inferred that hit2 and hit4, among the five hit molecules assessed, hold potential as more efficacious allosteric inhibitors of the hACE2–spike RBD interaction than SB27012.

Moving ahead, PCA of the whole spike RBD complexed with SB27012, hit2, and hit5 reveals larger eigenvalues and 2-D-projected subspaces (Fig. S10). Thus, it is perceived that the binding of these compounds in the allosteric pocket of hACE2 is responsible for the highest conformational movement of the spike protein. Similarly, other hit molecules bound to the hACE2 lead to large conformational movements in the spike RBD. Additionally, we provide contact maps for the reference system (Fig. 5E) and the systems containing the hit molecules (Fig. S11), highlighting the major contact regions affected. The residue contact analysis carried out to examine the changes in contact demonstrates that for the hit2 and hit4 complexes, there is a reduction in the number of contacts not only between hACE2 and the spike RBD but also within the spike RBD.

3.5. The effectiveness of SB27012 against SARS-CoV-2 VOCs

SARS-CoV-2 has evolved through mutations in the spike protein into several variants of concern (VOCs). Since SB27012 is an allosteric inhibitor of hACE2 that is unlikely to get mutated, it is expected that they ideally should work against all VOCs. As a proof of concept, we considered three systems with the spike RBD of three lethal VOCs. These variants, Beta, Delta, and early Omicron, were selected based on their structural availability and global prevalence at the time of model development. At that moment, the structural data for newer sublineages such as BA.1 and BA.2 were either unavailable or incomplete, which limited their inclusion in our simulations. We conducted 500 ns MD simulations for all VOCs complexed with hACE2–SB27012 and performed post-MD analysis.

It is worth noting from the superimposed structures that the Beta spike RBD has the highest deviation (Fig. 6A). Despite the predominantly smooth profiles exhibited by all complexes, the RMSD of the Beta variant displays the greatest deviation, indicating heightened structural perturbation, followed by Omicron, Delta, and wild-type variants (Fig. 6B). Similarly, the RMSF analysis reveals that the Omicron variant exhibits the most pronounced fluctuations, followed by the Beta, wild-type, and Delta variants (Fig. S12A). Furthermore, the calculation of SASA indicates that all variants, including the wild-type, have a comparable SASA (Fig. S12B). Regarding compactness, R_g calculations identify that the wild-type and Delta variants' RBD has the most compact structure, followed by the Beta and Omicron variants (Fig. S12C). Despite apparent discrepancies across analyses, a comprehensive evaluation shows that among the mutants, the spike RBD of the Beta variant experiences the most significant perturbation upon SB27012 binding, followed by the Omicron and Delta variants. The degrees of perturbation are higher than that of the wild-type. The PCA of the spike RBD of these complexes also shows similar results; i.e., the eigenvalues and the corresponding essential subspaces (Fig. 6C and Fig. S12D) capture larger conformational movements for the complexes of Omicron and Beta variants. Distinctly, the Delta spike RBD shows the lowest eigenvalue and subspace in the 2D projection plots.

Fig. 6 (A) Superimposed images of the mutated spike RBDs from the VOC complexes [wild-type: tan, Beta: light cyan, Delta: aero green, and Omicron: pink], (B) RMSD profiles of spike RBDs in the complexes, (C) the 2D projection plots of spike RBDs' first two principal eigenvectors, (D) the percentage of decrease in the binding energy of hACE2 and the spike RBD in VOCs shown as bar plots, and (E) RDF plots of the hACE2 α1 helix with respect to the two interfacial β sheets of RBDs.

To explore how affinity changes over the simulation period, we compared the MMPBSA binding energy between hACE2 and the spike RBD for the final 10 ns of the trajectories. Upon comparison with the initial 10 ns of their respective complexes (data not provided), we observed a higher percentage of decrease in binding affinity between the hACE2 and wild-type RBD, followed by the Delta, Omicron, and Beta variants (Fig. 6D), which corroborates the experimental findings. Thus, the RBD of the variants, especially Omicron, exhibits much higher binding affinity for the host hACE2 than the wild-type RBD, which is also supported by previous computational reports.^67–69 Based on the exhaustive analysis performed, the hit molecules identified to inhibit the wild-type spike RBD-hACE2 interaction may also be more effective against the spike variants compared to SB27012. The energy profiles are provided in Table S6.

To verify the binding energy results, we again performed the RDF analysis of the hACE2 α1 helix with respect to two interfacial β sheets of the RBD. A small probability peak at 0.17 nm appears for all variants except Beta, suggesting a close proximity between the spike RBD and hACE2 at some point during the simulations for all except Beta. Again, a second more intense probability peak occurs between 1.3 and 1.6 nm. This indicates the residence of hACE2 and RBD at a distance of 1.3–1.6 nm away from each other (Fig. 6E). Notably, the ascending loop in the plot shows a significant shift for all variants. Additionally, the SB27012 binding is not much affected by simulation as there is no apparent difference among the protein–ligand binding energies across all the complexes (Fig. S12E and Table S7).

Given that the allosteric and orthosteric sites of the receptor protein are conserved motifs, the identified hits (hit2 and hit4) should, in principle, be unanimously effective across all variants of concern. To test this hypothesis, we examined six systems of hACE2-hit2/4 complexes, each docked separately to the spike RBD of three variants, i.e., Beta, Delta, and Omicron. Each system was subjected to 500 ns of MD simulation; post-MD parameters were evaluated compared to those of the SB27012-bound complexes. The analysis (Fig. S13) shows that both hit molecules exhibit similar efficacy to SB27012 in perturbing the spike RBD. While both hits impact the Delta variant's RBD comparably to SB27012, significant conformational deviations were observed in the Beta and Omicron RBDs. MMPBSA binding free energy calculations between hACE2 and the RBD across all complexes also revealed a comparable difference for the Delta RBD. However, for both the Beta and Omicron variants, the hit molecules significantly reduce the binding affinity (Fig. S14). These findings support our hypothesis that hit2 and hit4 are more effective allosteric modulators, disrupting the hACE2–RBD interaction in the variants of concern.

3.6. Dynamic network analysis unraveling a path of allosteric modulation

To find out the intramolecular signaling pathway, we carried out dynamic network analysis (DNA), a method to study and analyze the dynamic behavior of a complex system over time, where the system is represented by a network of nodes connected by edges.⁷⁰ DNA representation of protein generally consists of the C_α atoms of the residues as their nodes, and the shortest optimal path between two different residues is considered to be the means of modulation.⁷¹ We have previously surmised Gln42–Gly446 as a vital interaction for hACE2–spike RBD binding. In order to discover an optimal pathway through which the SB27012 binding at the non-orthosteric pocket of hACE2 allosterically modulates its spike RBD binding at the α1 helix interface, a dynamic network map was generated and analyzed. In the network of hACE2 shown in Fig. 7A, the nodes are shown in blue except for the suboptimal route, which is shown in maroon color. The identified suboptimal or shortest path consists of 6 residues connecting Lys562 and Gln42. While Lys562 is a residue from the allosteric pocket that constantly interacts with not just SB27012 but all the hit molecules, Gln42 of the α1 helix interface interacts with Gly446 of the spike RBD interacting surface in the hACE2–spike RBD complex. The suboptimal path residues are shown in Fig. 7B as Lys562–Leu392–Phe390–Glu37–Asp38–Gln42.

Fig. 7 The dynamic network analysis of the hACE2–SB27012–spike RBD (wild-type) complex: (A) the shortest suboptimal path marked in maroon from Gln42 to Lys562 via 4 other residues, (B) a clear visualization of the suboptimal-path residues and the hACE2–spike RBD interface, (C) centrality analysis of the suboptimal-path residues, and (D) the contact map analysis of hACE2 showing contacts between the suboptimal path residues.

There could be multiple other paths of the same length, where only two out of six residues get replaced with other residues, such as Asp38 (Arg39, Leu40, and Tyr41) and Leu392 (Arg393), to generate alternative paths. This necessitated performing node centrality analysis to study the relative influence and prominence of the residues involved in the shortest path, illustrated in Fig. 7C. The analysis proved conclusive in confirming the most suitable path already proposed. The betweenness centrality plot suggests that the alternative residues are less likely to mediate protein conformational changes and regulate any form of allostery in comparison to the original residues. Lys562 displays the highest degree of closeness centrality, suggesting that it may play a critical role in determining the compactness of the protein. The alternative residues exhibit approximately the same level of closeness centrality as their counterparts. The degree centrality plot indicates the residues that are important for protein folding and facilitating interaction with other biomolecules. Among the alternative residues, Arg39 has the highest degree of centrality. Ultimately, eigenvector centrality suggests that Lys562 may serve as a key regulator of hACE2 dynamics as it is connected to more residues than other members along the path. In contrast, the alternative residues show an incomparable centrality. A contact map analysis can confirm the interaction between these residues (Fig. 7D). Considering all centrality-based outcomes, we hypothesize that all mentioned residues are competent to be part of the pathway that can produce a significant allosteric effect from the designated non-orthosteric site to the hACE2 peptidase domain that binds to the spike RBD.

Through the multiple sequence alignment (MSA) of the spike RBD of various SARS-CoV-2 variants as well as the hACE2 receptor of humans and closely related species, it was confirmed that hACE2 Gln42 is a conserved residue across the mentioned species, whereas the spike Gly446 undergoes a point mutation only in Omicron as Q446S (Fig. S15). These mutations cause a stronger binding between hACE2 and the spike RBD for Omicron. The residues constituting the suboptimal path of hACE2 for allosteric modulation remain also conserved.

Previous efforts have aimed to develop allosteric inhibitors of hACE2 that interfere with its interaction with the spike RBD.⁷² Our comprehensive in silico study identifies a potential allosteric pocket on hACE2 that may serve as a focal point for future investigations and the design of more effective modulators. Notably, SB27012 has also been found to support hACE2's interaction with its natural substrate, angiotensin II. In comparison, we identify two promising understudied compounds, hit2 and hit4, which demonstrate superior potential as allosteric modulators. SB27012, meanwhile, shows reduced efficacy against the three SARS-CoV-2 variants of concern. Importantly, while the proposed modulators weaken hACE2–RBD binding through allosteric structural changes, they do not disrupt AngII binding at the catalytic site, suggesting a dual therapeutic advantage. Additionally, we uncover a potential pathway of allosteric propagation linking the pocket to the spike-binding interface. This structural insight offers a promising strategy to counter evolving SARS-CoV-2 variants through variant-agnostic allosteric inhibition.

4. Conclusion

The emergence of new SARS-CoV-2 Omicron sub-variants strongly underscores the need for continued investigation on SARS-CoV-2. Deviating far from the conventional target-specific drug discovery, this research uniquely highlights the modulation of a host receptor protein to intercept pathogen entry into hosts. In our investigation, we identified a potential allosteric site on the human ACE2 receptor and elucidated the mechanism by which allostery is propagated from its non-orthosteric site to the spike RBD-interacting surface. The binding of the agonistic modulators induces counter-directional repulsion of hACE2 and the spike RBD, leading to large-scale perturbation of the RBD dynamics to the extent that it cannot form a stable complex any longer. However, these adversities for the spike RBD turn out to be beneficial for AngII, the natural substrate of hACE2 involved in many crucial functions. Surprisingly, the allosteric inhibition of spike RBD binding aids in better hACE2–AngII interaction. Eventually, the study led us to propose two molecules that can more effectively bind to the same allosteric site and produce stronger inhibitory effects on the hACE2–RBD interaction. We also showed that these molecules prove effective in disrupting complex formation for the RBDs of the three major variants of concern. Finally, we delineated a pathway comprising six residues through which allosteric signals can be selectively transduced to suppress the interaction with virus particles. Last but not least, we encourage rigorous experimental exploration of mechanistic insights that might help develop even more effective therapeutic options in the future.

Author contributions

Pratyush Pani: conceptualization, data curation, investigation, formal analysis, writing – original draft; Saroj Kumar Panda: data curation, formal analysis, writing – original draft; Malay Kumar Rana: supervision, writing – review & editing.

Conflicts of interest

The authors report no conflicts of interest.

Data availability

All the computational software and web servers used for this work have been appropriately cited in the manuscript with their version numbers. Gaussian16 and GaussView6 used (https://gaussian.com/) here have a commercial license. PyRx v0.8 was obtained from SOURCEFORGE (https://sourceforge.net/projects/pyrx/files/latest/download). The source code of GROMACS v2022 can be found at https://doi.org/10.5281/zenodo.6103835. The g_mmpbsa package was obtained using https://rashmikumari.github.io/g_mmpbsa/. Additional data have been provided as SI.

The SI contains further computational details, additional figures and tables, and relevant references. See DOI: https://doi.org/10.1039/d5cp01740h.

Acknowledgements

The authors acknowledge IISER Berhampur for computational support. PP thanks the University Grants Commission (UGC), New Delhi, India, for a Research Fellowship to conduct his doctoral study. SKP is also grateful for the Institute PhD Fellowship provided by IISER Berhampur to carry out this work. The authors thank Dr. Parth Sarthi Sen Gupta, DY Patil International University, for the valuable discussion and help at the beginning of this work.

References

1. WHO, Statement on the 1st meeting of the IHR Emergency Committee on the 2014 Ebola outbreak in West Africa.
2. T. Burki, WHO ends the COVID-19 public health emergency, Lancet Respir. Med., 2023, 11, 588.
3. A. Behl, A. Nair, S. Mohagaonkar, P. Yadav, K. Gambhir, N. Tyagi, R. K. Sharma, B. S. Butola and N. Sharma, Threat, challenges, and preparedness for future pandemics: A descriptive review of phylogenetic analysis based predictions, Infect., Genet. Evol., 2022, 98, 105217.
4. W. H. Organization, WHO COVID-19 dashboard, https://data.who.int/dashboards/covid19/cases?n=c.
5. O. J. Watson, G. Barnsley, J. Toor, A. B. Hogan, P. Winskill and A. C. Ghani, Global impact of the first year of COVID-19 vaccination: a mathematical modelling study, Lancet Infect. Dis., 2022, 22, 1293–1302.
6. K. I. Notarte, J. A. Catahay, J. V. Velasco, A. Pastrana, A. T. Ver, F. C. Pangilinan, P. J. Peligro, M. Casimiro, J. J. Guerrero, M. M. L. Gellaco, G. Lippi, B. M. Henry and C. Fernández-de-las-Peñas, Impact of COVID-19 vaccination on the risk of developing long-COVID and on existing long-COVID symptoms: A systematic review, eClinicalMedicine, 2022, 53, 101624.
7. A. A. Dawood, Mutated COVID-19 may foretell a great risk for mankind in the future, New Microbes New Infect., 2020, 35, 100673.
8. A. K. Chakraborty, Higher Omicron JN.1 Coronavirus Transmission due to Unique 17MPLF Spike Insertion compensating 24LPP, 69HV, 145Y, 211N and 483V deletions in the spike, Res. Sq., 2024 DOI:10.21203/rs.3.rs-3830998/v1.
9. J. Lan, J. Ge, J. Yu, S. Shan, H. Zhou, S. Fan, Q. Zhang, X. Shi, Q. Wang, L. Zhang and X. Wang, Structure of the SARS-CoV-2 spike receptor-binding domain bound to the ACE2 receptor, Nature, 2020, 581, 215–220.
10. Q. Wang, Y. Zhang, L. Wu, S. Niu, C. Song, Z. Zhang, G. Lu, C. Qiao, Y. Hu, K.-Y. Yuen, Q. Wang, H. Zhou, J. Yan and J. Qi, Structural and Functional Basis of SARS-CoV-2 Entry by Using Human ACE2, Cell, 2020, 181, 894–904.e9.
11. G. P. Pattnaik and H. Chakraborty, Entry Inhibitors: Efficient Means to Block Viral Infection, J. Membr. Biol., 2020, 253, 425–444.
12. S. A. Amin, S. Banerjee, K. Ghosh, S. Gayen and T. Jha, Protease targeted COVID-19 drug discovery and its challenges: Insight into viral main protease (Mpro) and papain-like protease (PLpro) inhibitors, Bioorg. Med. Chem., 2021, 29, 115860.
13. R. Cannalire, C. Cerchia, A. R. Beccari, F. S. Di Leva and V. Summa, Targeting SARS-CoV-2 Proteases and Polymerase for COVID-19 Treatment: State of the Art and Future Opportunities, J. Med. Chem., 2022, 65, 2716–2746.
14. M. Hoffmann, H. Kleine-Weber, S. Schroeder, N. Krüger, T. Herrler, S. Erichsen, T. S. Schiergens, G. Herrler, N.-H. Wu, A. Nitsche, M. A. Müller, C. Drosten and S. Pöhlmann, SARS-CoV-2 Cell Entry Depends on ACE2 and TMPRSS2 and Is Blocked by a Clinically Proven Protease Inhibitor, Cell, 2020, 181, 271–280.e8.
15. C. Muus, M. D. Luecken, G. Eraslan, L. Sikkema, A. Waghray, G. Heimberg, Y. Kobayashi, E. D. Vaishnav, A. Subramanian, C. Smillie, K. A. Jagadeesh, E. T. Duong, E. Fiskin, E. Torlai Triglia, M. Ansari, P. Cai, B. Lin, J. Buchanan, S. Chen, J. Shu, A. L. Haber, H. Chung, D. T. Montoro, T. Adams, H. Aliee, S. J. Allon, Z. Andrusivova, I. Angelidis, O. Ashenberg, K. Bassler, C. Bécavin, I. Benhar, J. Bergenstråhle, L. Bergenstråhle, L. Bolt, E. Braun, L. T. Bui, S. Callori, M. Chaffin, E. Chichelnitskiy, J. Chiou, T. M. Conlon, M. S. Cuoco, A. S. E. Cuomo, M. Deprez, G. Duclos, D. Fine, D. S. Fischer, S. Ghazanfar, A. Gillich, B. Giotti, J. Gould, M. Guo, A. J. Gutierrez, A. C. Habermann, T. Harvey, P. He, X. Hou, L. Hu, Y. Hu, A. Jaiswal, L. Ji, P. Jiang, T. S. Kapellos, C. S. Kuo, L. Larsson, M. A. Leney-Greene, K. Lim, M. Litviňuková, L. S. Ludwig, S. Lukassen, W. Luo, H. Maatz, E. Madissoon, L. Mamanova, K. Manakongtreecheep, S. Leroy, C. H. Mayr, I. M. Mbano, A. M. McAdams, A. N. Nabhan, S. K. Nyquist, L. Penland, O. B. Poirion, S. Poli, C. Qi, R. Queen, D. Reichart, I. Rosas, J. C. Schupp, C. V. Shea, X. Shi, R. Sinha, R. V. Sit, K. Slowikowski, M. Slyper, N. P. Smith, A. Sountoulidis, M. Strunz, T. B. Sullivan, D. Sun, C. Talavera-López, P. Tan, J. Tantivit, K. J. Travaglini, N. R. Tucker, K. A. Vernon, M. H. Wadsworth, J. Waldman, X. Wang, K. Xu, W. Yan, W. Zhao, C. G. K. Ziegler, G. H. Deutsch, J. Dutra, K. J. Gaulton, J. Holden-Wiltse, H. L. Huyck, T. J. Mariani, R. S. Misra, C. Poole, S. Preissl, G. S. Pryhuber, L. Rogers, X. Sun, A. Wang, J. A. Whitsett, Y. Xu, J. Alladina, N. E. Banovich, P. Barbry, J. E. Beane, R. P. Bhattacharyya, K. E. Black, A. Brazma, J. D. Campbell, J. L. Cho, J. Collin, C. Conrad, K. de Jong, T. Desai, D. Z. Ding, O. Eickelberg, R. Eils, P. T. Ellinor, A. Faiz, C. S. Falk, M. Farzan, A. Gellman, G. Getz, I. A. Glass, A. Greka, M. Haniffa, L. P. Hariri, M. W. Hennon, P. Horvath, N. Hübner, D. T. Hung, H. L. Huyck, W. J. Janssen, D. Juric, N. Kaminski, M. Koenigshoff, G. H. Koppelman, M. A. Krasnow, J. A. Kropski, M. Kuhnemund, R. Lafyatis, M. Lako, E. S. Lander, H. Lee, M. E. Lenburg, C.-H. Marquette, R. J. Metzger, S. Linnarsson, G. Liu, Y. M. D. Lo, J. Lundeberg, J. C. Marioni, S. A. Mazzilli, B. D. Medoff, K. B. Meyer, Z. Miao, A. V. Misharin, M. C. Nawijn, M. Z. Nikolić, M. Noseda, J. Ordovas-Montanes, G. Y. Oudit, D. Pe’er, J. E. Powell, S. R. Quake, J. Rajagopal, P. R. Tata, E. L. Rawlins, A. Regev, M. E. Reid, P. A. Reyfman, K. M. Rieger-Christ, M. Rojas, O. Rozenblatt-Rosen, K. Saeb-Parsy, C. Samakovlis, J. R. Sanes, H. B. Schiller, J. L. Schultze, R. F. Schwarz, A. V. Segre, M. A. Seibold, C. E. Seidman, J. G. Seidman, A. K. Shalek, D. P. Shepherd, R. Sinha, J. R. Spence, A. Spira, X. Sun, E. Sundström, S. A. Teichmann, F. J. Theis, A. M. Tsankov, L. Vallier, M. van den Berge, T. A. Van Zyl, A.-C. Villani, A. Weins, R. J. Xavier, A. Ö. Yildirim, L.-E. Zaragosi, D. Zerti, H. Zhang, K. Zhang and X. Zhang, N. L. C. The and N. The Human Cell Atlas Lung Biological, Single-cell meta-analysis of SARS-CoV-2 entry genes across tissues and demographics, Nat. Med., 2021, 27, 546–559.
16. W. Li, M. J. Moore, N. Vasilieva, J. Sui, S. K. Wong, M. A. Berne, M. Somasundaran, J. L. Sullivan, K. Luzuriaga, T. C. Greenough, H. Choe and M. Farzan, Angiotensin-converting enzyme 2 is a functional receptor for the SARS coronavirus, Nature, 2003, 426, 450–454.
17. M. Donoghue, F. Hsieh, E. Baronas, K. Godbout, M. Gosselin, N. Stagliano, M. Donovan, B. Woolf, K. Robison, R. Jeyaseelan, R. E. Breitbart and S. Acton, A Novel Angiotensin-Converting Enzyme–Related Carboxypeptidase (ACE2) Converts Angiotensin I to Angiotensin 1-9, Circ. Res., 2000, 87, e1–e9.
18. L. T. Skeggs Jr., J. R. Kahn and N. P. Shumway, THE PREPARATION AND FUNCTION OF THE HYPERTENSIN-CONVERTING ENZYME, J. Exp. Med., 1956, 103, 295–299.
19. S. Keidar, M. Kaplan and A. Gamliel-Lazarovich, ACE2 of the heart: From angiotensin I to angiotensin (1–7), Cardiovasc. Res., 2007, 73, 463–469.
20. F. Jiang, J. Yang, Y. Zhang, M. Dong, S. Wang, Q. Zhang, F. F. Liu, K. Zhang and C. Zhang, Angiotensin-converting enzyme 2 and angiotensin 1–7: novel therapeutic targets, Nat. Rev. Cardiol., 2014, 11, 413–426.
21. M. M. Gironacci, H. P. Adamo, G. Corradi, R. A. Santos, P. Ortiz and O. A. Carretero, Angiotensin (1-7) Induces Mas Receptor Internalization, Hypertension, 2011, 58, 176–181.
22. Y. Hu, L. Liu and X. Lu, Regulation of Angiotensin-Converting Enzyme 2: A Potential Target to Prevent COVID-19?, Front. Endocrinol., 2021, 12 DOI:10.3389/fendo.2021.725967.
23. S. Xiao, H. Tian and P. Tao, PASSer2.0: Accurate Prediction of Protein Allosteric Sites Through Automated Machine Learning, Front. Mol. Biosci., 2022, 9 DOI:10.3389/fmolb.2022.879251.
24. W. Huang, S. Lu, Z. Huang, X. Liu, L. Mou, Y. Luo, Y. Zhao, Y. Liu, Z. Chen, T. Hou and J. Zhang, Allosite: a method for predicting allosteric sites, Bioinformatics, 2013, 29, 2357–2359.
25. A. F. Abdel-Magid, Allosteric Modulators: An Emerging Concept in Drug Discovery, ACS Med. Chem. Lett., 2015, 6(2), 104–107.
26. D. Mannar, J. W. Saville, Z. Sun, X. Zhu, M. M. Marti, S. S. Srivastava, A. M. Berezuk, S. Zhou, K. S. Tuttle, M. D. Sobolewski, A. Kim, B. R. Treat, P. M. Da Silva Castanha, J. L. Jacobs, S. M. Barratt-Boyes, J. W. Mellors, D. S. Dimitrov, W. Li and S. Subramaniam, SARS-CoV-2 variants of concern: spike protein mutational analysis and epitope for broad neutralization, Nat. Commun., 2022, 13, 4696.
27. J. W. Saville, D. Mannar, X. Zhu, S. S. Srivastava, A. M. Berezuk, J. P. Demers, S. Zhou, K. S. Tuttle, I. Sekirov, A. Kim, W. Li, D. S. Dimitrov and S. Subramaniam, Structural and biochemical rationale for enhanced spike protein fitness in delta and kappa SARS-CoV-2 variants, Nat. Commun., 2022, 13, 742.
28. D. Mannar, J. W. Saville, X. Zhu, S. S. Srivastava, A. M. Berezuk, K. S. Tuttle, A. C. Marquez, I. Sekirov and S. Subramaniam, SARS-CoV-2 Omicron variant: Antibody evasion and cryo-EM structure of spike protein–ACE2 complex, Science, 2022, 375(6582), 760–764.
29. H. Berman, K. Henrick and H. Nakamura, Announcing the worldwide Protein Data Bank, Nat. Struct. Mol. Biol., 2003, 10, 980.
30. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson and H. Nakatsuji, Gaussian 16, Gaussian, Inc., Wallingford CT.
31. R. Dennington, T. A. Keith and J. M. Millam, GaussView 6, 2016.
32. P. S. Sen Gupta, H. R. Bhat, S. Biswal and M. K. Rana, Computer-aided discovery of bis-indole derivatives as multi-target drugs against cancer and bacterial infections: DFT, docking, virtual screening, and molecular dynamics studies, J. Mol. Liq., 2020, 320, 114375.
33. O. Trott and A. J. Olson, AutoDock Vina: Improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading, J. Comput. Chem., 2010, 31, 455–461.
34. N. T. Nguyen, T. H. Nguyen, T. N. H. Pham, N. T. Huy, M. Van Bay, M. Q. Pham, P. C. Nam, V. V. Vu and S. T. Ngo, Autodock Vina Adopts More Accurate Binding Poses but Autodock4 Forms Better Binding Affinity, J. Chem. Inf. Model., 2020, 60, 204–211.
35. S. Dallakyan and A. J. Olson, in Chemical Biology: Methods and Protocols, ed. J. E. Hempel, C. H. Williams and C. C. Hong, Springer, New York, NY, 2015, pp. 243–250.
36. C. Dominguez, R. Boelens and A. M. J. J. Bonvin, HADDOCK:  A Protein−Protein Docking Approach Based on Biochemical or Biophysical Information, J. Am. Chem. Soc., 2003, 125, 1731–1737.
37. J. Sunseri and D. R. Koes, Pharmit: interactive exploration of chemical space, Nucleic Acids Res., 2016, 44, W442–W448.
38. T. I. Oprea and H. Matter, Integrating virtual screening in lead discovery, Curr. Opin. Chem. Biol., 2004, 8, 349–358.
39. A. Daina, O. Michielin and V. Zoete, SwissADME: a free web tool to evaluate pharmacokinetics, drug-likeness and medicinal chemistry friendliness of small molecules, Sci. Rep., 2017, 7, 42717.
40. H. Yang, C. Lou, L. Sun, J. Li, Y. Cai, Z. Wang, W. Li, G. Liu and Y. Tang, admetSAR 2.0: web-service for prediction and optimization of chemical ADMET properties, Bioinformatics, 2019, 35, 1067–1069.
41. C. A. Lipinski, F. Lombardo, B. W. Dominy and P. J. Feeney, Experimental and computational approaches to estimate solubility and permeability in drug discovery and development settings1PII of original article: S0169-409X(96)00423-1. The article was originally published in Advanced Drug Delivery Reviews 23 (1997) 3, Adv. Drug Deliv. Rev., 2001, 46, 3–26.
42. M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess and E. Lindahl, GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers, SoftwareX, 2015, 1–2, 19–25.
43. A. W. Sousa da Silva and W. F. Vranken, ACPYPE - AnteChamber PYthon Parser interfacE, BMC Res. Notes, 2012, 5, 367.
44. J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman and D. A. Case, Development and testing of a general amber force field, J. Comput. Chem., 2004, 25, 1157–1174.
45. V. Hornak, R. Abel, A. Okur, B. Strockbine, A. Roitberg and C. Simmerling, Comparison of multiple Amber force fields and development of improved protein backbone parameters, Proteins:Struct., Funct., Bioinf., 2006, 65, 712–725.
46. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, Comparison of simple potential functions for simulating liquid water, J. Chem. Phys., 1983, 79, 926–935.
47. E. J. Haug, J. S. Arora and K. Matsui, A steepest-descent method for optimization of mechanical systems, J. Optim. Theory Appl., 1976, 19, 401–424.
48. D. Frenkel and B. Smit, Understanding molecular simulation: from algorithms to applications, 2nd edn, 1996, vol. 50.
49. H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola and J. R. Haak, Molecular dynamics with coupling to an external bath, J. Chem. Phys., 1984, 81, 3684–3690.
50. M. Parrinello and A. Rahman, Polymorphic transitions in single crystals: A new molecular dynamics method, J. Appl. Phys., 1981, 52, 7182–7190.
51. U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee and L. G. Pedersen, A smooth particle mesh Ewald method, J. Chem. Phys., 1995, 103, 8577–8593.
52. J. E. Lennard-Jones, Cohesion, Proc. Phys. Soc., 1931, 43, 461.
53. B. Hess, H. Bekker, H. J. C. Berendsen and J. G. E. M. Fraaije, LINCS: A linear constraint solver for molecular simulations, J. Comput. Chem., 1997, 18, 1463–1472.
54. J. Srinivasan, T. E. Cheatham, P. Cieplak, P. A. Kollman and D. A. Case, Continuum Solvent Studies of the Stability of DNA, RNA, and Phosphoramidate−DNA Helices, J. Am. Chem. Soc., 1998, 120, 9401–9409.
55. S. Genheden and U. Ryde, The MM/PBSA and MM/GBSA methods to estimate ligand-binding affinities, Expert Opin. Drug Discovery, 2015, 10, 449–461.
56. R. Kumari, R. Kumar and A. Lynn, g_mmpbsa—A GROMACS Tool for High-Throughput MM-PBSA Calculations, J. Chem. Inf. Model., 2014, 54, 1951–1962.
57. V. Govind Kumar, A. Polasa, S. Agrawal, T. K. S. Kumar and M. Moradi, Binding affinity estimation from restrained umbrella sampling simulations, Nat. Comput. Sci., 2023, 3, 59–70.
58. G. M. Süel, S. W. Lockless, M. A. Wall and R. Ranganathan, Evolutionarily conserved networks of residues mediate allosteric communication in proteins, Nat. Struct. Biol., 2003, 10, 59–69.
59. B. Chakrabarty and N. Parekh, NAPS: Network Analysis of Protein Structures, Nucleic Acids Res., 2016, 44, W375–W382.
60. B. Chakrabarty, V. Naganathan, K. Garg, Y. Agarwal and N. Parekh, NAPS update: network analysis of molecular dynamics data and protein–nucleic acid complexes, Nucleic Acids Res., 2019, 47, W462–W470.
61. Y.-H. Shin, K. Jeong, J. Lee, H. J. Lee, J. Yim, J. Kim, S. Kim and S. B. Park, Inhibition of ACE2-Spike Interaction by an ACE2 Binder Suppresses SARS-CoV-2 Entry, Angew. Chem., Int. Ed., 2022, 61, e202115695.
62. A. Volkamer, D. Kuhn, T. Grombacher, F. Rippmann and M. Rarey, Combining Global and Local Measures for Structure-Based Druggability Predictions, J. Chem. Inf. Model., 2012, 52, 360–372.
63. A. Volkamer, A. Griewel, T. Grombacher and M. Rarey, Analyzing the Topology of Active Sites: On the Prediction of Pockets and Subpockets, J. Chem. Inf. Model., 2010, 50, 2041–2052.
64. E. Krissinel and K. Henrick, Inference of Macromolecular Assemblies from Crystalline State, J. Mol. Biol., 2007, 372(3), 774–797.
65. M. P. Doogue and T. M. Polasek, The ABCD of clinical pharmacokinetics, Ther. Adv. Drug Saf., 2013, 4, 5–7.
66. K. P. Tan, K. Singh, A. Hazra and M. S. Madhusudhan, Peptide bond planarity constrains hydrogen bond geometry and influences secondary structure conformations, Curr. Res. Struct. Biol., 2021, 3, 1–8.
67. C. H. S. da Costa, C. A. B. de Freitas, C. N. Alves and J. Lameira, Assessment of mutations on RBD in the Spike protein of SARS-CoV-2 Alpha, Delta and Omicron variants, Sci. Rep., 2022, 12, 8540.
68. X. Zhang, B. Hong, P. Wei, P. Pei, H. Xu, L. Chen, Y. Tong, J. Chen, S.-Z. Luo, H. Fan and C. He, Pathogen-host adhesion between SARS-CoV-2 spike proteins from different variants and human ACE2 studied at single-molecule and single-cell levels, Emerg. Microbes Infect., 2022, 11, 2658–2669.
69. C. S. Lupala, Y. Ye, H. Chen, X.-D. Su and H. Liu, Mutations on RBD of SARS-CoV-2 Omicron variant result in stronger binding to human ACE2 receptor, Biochem. Biophys. Res. Commun., 2022, 590, 34–41.
70. A. Majdandzic, B. Podobnik, S. V. Buldyrev, D. Y. Kenett, S. Havlin and H. Eugene Stanley, Spontaneous recovery in dynamical networks, Nat. Phys., 2014, 10, 34–38.
71. J. C. F. Ng and F. Fraternali, Understanding the structural details of APOBEC3-DNA interactions using graph-based representations, Curr. Res. Struct. Biol., 2020, 2, 130–143.
72. J. E. Hochuli, S. Jain, C. Melo-Filho, Z. L. Sessions, T. Bobrowski, J. Choe, J. Zheng, R. Eastman, D. C. Talley, G. Rai, A. Simeonov, A. Tropsha, E. N. Muratov, B. Baljinnyam and A. V. Zakharov, Allosteric Binders of ACE2 Are Promising Anti-SARS-CoV-2 Agents, ACS Pharmacol. Transl. Sci., 2022, 5, 468–478.

---

Footnote

† Authors contributed equally.

---