Studying the Effects of ACE2 Mutations on the Stability, Dynamics, and Dissociation Process of SARS-CoV-2 S1/hACE2 Complexes.

A highly infectious coronavirus, SARS-CoV-2, has spread in many countries. This virus recognizes its receptor, angiotensin-converting enzyme 2 (ACE2), using the receptor binding domain of its spike protein subunit S1. Many missense mutations are reported in various human populations for the ACE2 gene. In the current study, we predict the affinity of many ACE2 variants for binding to S1 protein using different computational approaches. The dissociation process of S1 from some variants of ACE2 is studied in the current work by molecular dynamics approaches. We study the relation between structural dynamics of ACE2 in closed and open states and its affinity for S1 protein of SARS-CoV-2.

Introduction

Since the outbreak of coronavirus disease (COVID-19) in December 2019 in Wuhan, China, more than 210 countries have become involved. This highly infectious disease is caused by a coronavirus named severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) by the International Committee on Taxonomy of Viruses (ICTV).

Coronaviruses are enveloped RNA viruses. They can cause respiratory, enteric, and nervous system diseases in birds and mammals.[1] SARS-CoV-2 is a β-coronavirus. Among four genera of coronaviruses, β-coronaviruses contain species that have caused serious human respiratory tract infections: severe acute respiratory syndrome (SARS) with outbreaks in 2002–2004 and Middle East respiratory syndrome (MERS) with outbreaks from 2012, which were about 10% and 36% lethal, respectively.[2]

The main structural proteins of coronaviruses located in the viral envelope are membrane (M), envelope (E), and spike (S) proteins. The S proteins form homo-trimers on the virus surface and are necessary for entrance of the virus into the host cell. They determine the host range and tissue tropism and induce host immune response. Each S glycoprotein consists of the ectodomain, transmembrane, and intracellular regions. The ectodomain consists of S1 and S2 subunits used for attachment to receptor proteins and fusion with the membrane of host cells, respectively. The S1 subunit of the spike protein is composed of a C-terminal domain (CTD) and N-terminal domain (NTD). CTD, which is important in binding of S proteins to protein receptors, is composed of a core region and a receptor binding domain (RBD).[2]

The spike protein of SARS-CoV-2 virus attaches to human angiotensin-converting enzyme 2 (hACE2),[3−6] which is also the target of SARS coronavirus (SARS-CoV).[7−9] ACE2 is expressed in different organs of the body including lungs, intestines, heart, kidneys, and endothelium.[10−13] This protein is a zinc peptidase involved in converting angiotensin I to angiotensin 1-9,[14,15] and angiotensin II to angiotensin 1-7.[11]

The ectodomain of ACE2 is composed of a collectrin homology domain and a peptidase domain more distant from the cell membrane. The peptidase domain contains subdomains I and II that contain N-terminal and C-terminal regions of the active site cleft, respectively.[16] Virus-binding motifs of hACE2 are recognized on the outer N-terminal surface away from the catalytic cleft.[17] Similar to that observed for SARS-CoV,[17,18] the concave RBD of SARS-CoV-2 interacts with the convex N-terminal helix of hACE2 mainly through polar interactions.[19−21]

Changes in coronavirus spike proteins could alter their corresponding host.[22−26] In the same manner, the affinity between SARS-CoV spike protein and ACE2 in the first step of viral attachment was a determining factor in viral infectivity and host susceptibility and transmissions.[27−35] A few variations in critical ACE2 residues caused a lower affinity of SARS-CoV virus to mouse, rat, and Daubenton’s bat ACE2 than to human ACE2.[17,27,35,36] Variations in critical amino acids of mammalian ACE2s were predicted to decrease or increase SARS-CoV-2 recognition.[37] Therefore, it is possible that mutations in hACE2 affect the affinity and infection severity of SARS-CoV-2 in human hosts.

Iran is one of the countries infected pandemically with the virus. Iranian people show a variety of genetic backgrounds throughout the country.[38] Beside global mutation databases such as ENSEMBL,[39] Iranome is a database that contains genetic information on different ethnic groups in Iran.[40] On the basis of Iranome data, the hACE2 gene shows different mutations in different ethnic groups.

To find out whether the hACE2 missense mutations observed in human populations affect the affinity of this protein to SARS-CoV-2 protein, we created the three-dimensional (3D) structures of the mutated hACE2 proteins and computed the affinity of these proteins to the SARS-CoV-2 spike S1 protein. The interaction between S1 and its receptor in human cells is one of the interesting targets for drug design to fight the virus in the earliest stages. It is believed that blocking the virus–host interaction at the molecular level opens new windows toward therapy. In the current study, we considered the effect of the reported mutations in ACE2 on its binding affinity for the S1 protein using some computational approaches. We studied the dissociation of S1 RBD of SARS-CoV-2 from the extracellular domain of different mutants of the hACE2 protein in silico. The described dissociation process sheds light on the association of the SARS-CoV-2 S1 subunit with the ACE2 protein. We studied the dynamics and stability of some ACE2 variants using molecular dynamics approaches. The possible relation between ACE2 dynamics and its affinity for SARS-CoV-2 S1 subunit is discussed in the current study.

Methods

3D Structures

The 3D structure of RBD of SARS-CoV-2 spike subunit 1 complexed with its receptor, human ACE2, was retrieved from the Protein Data Bank (PDBID 6VW1).[19] The reported structure is chimeric. Its receptor binding motif (RBM) that interacts with ACE2 belongs to the SARS-CoV-2 virus. The 3D structures of human ACE2 mutants in association with spike RBD are generated using the FoldX suite.[41] The FoldX suite is also utilized to compute the folding stability of ACE2 along dissociation simulations. To simulate the dissociation process of S1 from wild type (WT) ACE2 in its closed state, we utilized ACE2 bonded to an ACE2-specefic inhibitor, ORE-1001 (PDBID 1R4L).[16] Then, we superimpose and align the ACE2 structure in the closed state on wild-type ACE2 in association with S1 to find the primary position of S1 in association with the closed state of ACE2.

Missense Mutations

The position of missense mutations in the hACE2 gene is derived from ENSEMBL release 99, January 2020 with accession code ENST00000427411.1.[39] This list covers the reported mutations in hACE2 for various human populations with different geographical distributions. Also, some missense mutations of ACE2 specific for Iranian ethnic groups are derived from the Iranome project public data[40] and ENSEMBL database.[39]

Prediction of ACE2 Affinity to S1 Protein Using Fast Methods

We predicted the affinity of the ACE2 mutants to the RBD of S1 by PISA,[42] FoldX,[43] Prodigy,[44] SEPAS,[45] and SAAMBE-3D.[46] We fed the noted tools with complexes whose ACE2 subunits are mutated.

PISA-Based Prediction of Complex Stability

PISA assigns the dissociation free energy (ΔGdiss) to the dimers using the following equation:[47]The solv, rot, and trans subscripts stand for solvation, rotation, and transition, respectively. Number (N) of hydrogen bonds (hb), salt bridges (sb), and disulfide bonds (db) are also considered in the PISA computation. Nonbonded interaction energy is presented by E. Δσ stands for the buried surface area during complex formation. F and C are fitting constants. The hydration energy and other contact-dependent energies are implemented in ΔGint as the binding energy of subunits. Rigid body entropy (ΔS) is considered for rotation (rot) and translation (trans) terms. We recruit the command line version of PISA modules implemented in the CCP4 suite.[48] In the current study, the input structures for PISA to predict the dissociation free energy of the complex are FoldX-generated mutant 3D structures of ACE2 in association with the RBD of the S1 protein. If the deviation of PISA-computed stability of a mutant form of ACE2/S1 complex from the stability computed for WT ACE2/S1 was more than 0.1 kcal/mol, we considered the mutation as one that changed the free energy of ACE2/S1 dissociation.

FoldX-Based Prediction of Complex Stability

FoldX computes the stability of protein complexes considering monomer and dimer stabilities. The computed complex stability denotes the intersubunit affinity. The stability of each monomer is computed by using the empirical terms:[41]This equation is a linear combination of various terms with specific coefficients. The van der Waals (VDW) term’s contribution to total energy, solvation of hydrophobic and polar groups, water bridges, hydrogen bonds, electrostatic interactions, and electrostatic interactions’ contribution to association constant are abbreviated to vdw, solvH, solvP, wb, hbond, el, and kon in subscripts of stability terms, respectively. The translational source of entropy, main-chain-, and side-chain-defined entropies are abbreviated to tr, mc, and sc in subscripts of ΔS, respectively. We utilize the stand-alone academic version of FoldX. The input structures of FoldX to predict the intersubunit affinity are FoldX-generated mutant 3D structures of ACE2 in association with the RBD of the S1 subunit. If the deviation of FoldX-computed stability of a mutant form of ACE2/S1 complex from the stability computed for WT ACE2/S1 was more than 0.1 kcal/mol, we considered the mutation as one that changed the affinity of ACE2 for the S1 protein.

Prodigy-Based Prediction of ACE2 Affinity for RBD of S1 Subunit

Prodigy utilizes inter-residue contacts between subunits and the noninterface surface to predict the binding affinity between subunits using the following equation:[44]ICs stands for inter-residue contacts, and the subscripts indicate the type of contact or type of noninteraction surface (NIS). In the current study, we used the stand-alone version of Prodigy.[49] The FoldX-generated 3D structures of ACE2 mutants in complex with S1 are the inputs of Prodigy. If the deviation of Prodigy-computed affinity of a mutant form of ACE2 to S1 from the affinity of WT ACE2 to S1 was more than 0.1 kcal/mol, we considered the mutation as one that changed the intersubunit affinity.

SAAMBE3D-Based Prediction of Complex Stability

SAAMBE was introduced to predict the effect of mutations on dimer stability. It considers MM/PBSA and knowledge-based descriptors such as amino-acid-specific dielectric constants to perform its predictions. SAAMBE-3D is the next generation of the mentioned method that uses a machine learning approach to compute affinity changes resulting from mutations using many knowledge-based properties of residues including volume, hydrophobicity, flexibility, polarity, distance-based metrics, and many other descriptors. We utilized the stand-alone version of the mentioned software.[46] The FoldX-constructed 3D structures of ACE2 mutants in complex with RBD of S1 are the inputs of SAAMBE-3D in the current study. On the basis of the software criteria, the positive binding energy of a mutant means that such a mutation weakens the binding free energy. If the amount of changes in binding energy was more than 0.1 kcal/mol, we considered the mutation as one that changed the intersubunit affinity.

SEPAS-Based Prediction of Complex Stability

The mentioned methods all require the 3D structure of the protein complex to predict affinity or changes in affinity between subunits. SEPAS utilizes the 3D structure of one monomer to perform affinity prediction using the mechanical softness of the protein binding patch (PBP). It requires the 3D structure of a PBP region on one monomer to predict the affinity between the introduced PBP and tentative partners of the PBP. It predicts the softness of PBP using the following equation:[45,50]N represents the count of residues in PBP, and Nn represents number of unique intra-PBP native contacts. The cutoff distance for defining the native contact is set at 7.5 Å. To create ensembles of near-equilibrium structures for all 240 mutants of ACE2 reported in ENSEMBL and Iranome that are considered in current study, we utilize the anisotropic network model (ANM) of the complexes.[51,52] The ANM approach provides a pool of 1000 structures for each of 240 mutants of ACE2 bonded to the S1 protein. The SEPAS-assigned affinity of ACE2 for S1 is the average of predicted affinities of 1000 samples generated by ANM for each mutant structure. Mainly, SEPAS predicts the stability in a classwise manner. It predicts four possibilities for each PBP, i.e., how much the introduced PBP belongs to high-affinity, high-medium, medium-low, and low-affinity classes. If a mutation in ACE2 leads to a 5% or more change in distribution of the computed possibilities among denoted classes with respect to WT, we count the mutation as an effective one.

Adaptive Tempering MD and MM-GBSA Approach

Free energy of ACE2 binding to RBD of SARS-CoV-2 spike protein 1 is also estimated in the current work by using MM-GBSA approach for mutations in ACE2 that are reported for Iranian ethnic groups by NAMD 2.13.[53,54] To perform MM-GBSA, we gradually heat up the complex structures to 300 K in Generalized Born implicit solvent (GBIS) after minimization steps. The ion concentration is set at 0.3 molar. The hydrophobic energy contribution from implicit solvent is considered in simulations. We use adaptive tempering molecular dynamics simulation (AT-MD) to enhance the system sampling.[55] The adaptive tempering is a single copy version of replica exchange MD. In the current work, we let the langevin thermostat use the updated temperatures from adaptive tempering for 2 × 106 steps. The resulted ensemble of structures was also used to study the properties of SARS-CoV-2 RBD/ACE2 complex for WT and ACE2 mutants. To estimate the binding energy using MM-GBSA, we utilized the single-trajectory approach[56] by considering the next equation:We used run-time-analysis feature of NAMD 2.13 to compute potential energy and its constituent terms in GBIS condition for all frames of protein complex (ΔG(ACE2/S1)), separated S1 subunit (ΔG(S1)), and separated ACE2 (ΔG(ACE2)). The run-time analysis is performed for 1000 snapshots derived from 2 × 106 steps of AT-MD simulation for WT and mutant complexes. The presented MM-GBSA-predicted affinities of ACE2 for S1 are the average values of the ensemble data.

Adaptive Biasing Force Approach

To simulate the binding of SARS-CoV-2 S1 protein to human ACE2 and to predict the binding free energy with higher accuracy, we utilized adaptive biasing force (ABF) for complexes between RBD of SARS-CoV-2 S1 protein and reported mutants of ACE2.[57,58] In ABF, we gradually dissociated monomers of the complexes along the z-axis of the complexes in a 4–18.5 ns simulation under GBIS conditions using collective variable features of NAMD 2.13. ABF is a history-dependent approach to compute free energy profiles along predefined collective variables. During simulation steps, the free energy surface is smoothed by the adaptive bias so the system behaves along the defined collective variable, the z-axis of the bonded subunits, like when it is able to perform diffusion there. In the current study, the bin size was set to 0.2 Å. We tested 0.2, 0.4, and 2 ps as the simulation time before applying the biasing force by ABF in each bin. Because the standard error of the computed forces and the computed potentials are not so small in long than in short bin-scanning and systems experience some binding/unbinding during some of simulations, here we set 0.2 ps as the simulation time before applying the biasing force by the ABF algorithm in each bine. We are interested in comparative rather than quantitative affinities. The ABF simulation time varies from 4 to 18.5 ns depending on the considered system. The total time of ABF simulations in the current study reached 248 ns. The maximum standard error (SE) of the ABF-computed free energy between two adjacent positions on the reaction coordinate (points c and b) are computed by the following equation:[59,60]The variance of the ABF force is defined by σ2; N denotes number of points that are sampled to compute PMF; and κ presents the correlation length of the calculated force. We compute the SEs for ABF simulations in two scenarios: time-based and distance-based. For computing the time-based SE, we split the ABF simulations into windows with 20 ps length. Each window has a 0.2 ps overlap with its previous window. For computing the distance-based SE, we split the distance between two subunits into windows with 0.5 Å length. This covers 2.5 bins of the original ABF simulations.

Characterizing Partially Melted Structures

The AT-MD-reported potential energies are used to measure the progression of systems along a partial unfolding path. The state that resides in the most negative potential well is considered as the most stable species. The reported potential energy and their corresponding temperature are used to define the temperature of microtransition and the energy barrier between folded and partially melted states of the protein using the following equation:[55]P stands for the normalized computed potential energy where 1 means folded state and 0 represents melted state, w adjusts the reported fraction of the melted state, T stands for the transition temperature, β denotes the reduced dimension of simulation temperature, and ΔE represents the energy difference between folded and partially melted states. We use Cα-based inter-residue distance root mean square (dRMS) deviation as a structural descriptor of the folded state.[61,62]

Miscellaneous

Some structural features like computation of dipole, salt-bridges, H-bonds, and ASA properties were computed by house-made Tcl codes run under VMD 1.9.2.[63] All MD-based computations such as pair interaction calculations, adaptive tempering, MM/GBSA, ABF, and other runtime analyses were performed by using NAMD 2.13. The utilized force fields for proteins are CHARMM 22 and 27.[64] In all simulations, the time steps are set at 2 fs. The force filed parameter for the drug molecule is utilized from SwissParam.[65] Prody is recruited to perform ANM and PRS computations.[52] The electrostatic potential surface for the complexes was computed by PMEPOT plugin of VMD 1.9.2.

Results and Discussion

There are reports indicating that SARS-CoV-2 utilizes ACE2 as its receptor for attaching to human cells.[4,7] ACE2 is expressed in a wide range of tissues such as heart, lungs, and intestines.[13] It is indicated that SARS-CoV-2 binds to its receptor via the S1 subunit of its surface spike complex.[20] The affinity of S1 for its receptor affects the virulence, its hosts, and the intensity of virus infection.[9] It is observed that in some populations the incidence of COVID-19 is not similar between members. Beside variables such as age, sex, and sanitary conditions, the effects of genetic variations are not ignorable.[66,67] The consequence of some mutations in ACE2 in its interaction with SARS-CoV or SARS-CoV-2 spikes was studied by genetic approaches.[67,68] The effects of some ACE2 mutations that reside in interface region of SARS-CoV-2 spike/ACE2 complex were studied by bioinformatics tools.[69,70] Studies suggested that ACE2 mutants possibly related to COVID-19 incidence.[71] Here, we studied the effect of ACE2 gene polymorphisms on its stability, dynamics, and binding affinity for the S1 protein of SARS-CoV-2 using multiple computational approaches. In the current study, we predict the effect of 240 widespread missense mutations in the ACE2 gene reported for different human populations and especially eight ones specific for Iranian ethnic populations on the binding affinity between ACE2 and S1 RBD of SARS-CoV-2 with different computational approaches from bioinformatic methods to a thermodynamic integration procedure.

There are many reported SNPs for the ACE2 gene in different human populations in the ENSEMBL database. We selected the missense SNPs in the human ACE2 extracellular domain. Recently, the 3D structure of RBD of the SARS-CoV-2 S1 protein in association with ACE2 was reported.[19] Using the crystallized 3D structure, we introduced the ENSEMBL-reported missense mutations in the ACE2 subunit of the assembly utilizing the FoldX suite.[41] The constructed mutated ACE2 subunits in association with S1 created the building blocks for the next steps in the current study.

Predicting the Affinity of ACE2 Mutants to SARS-CoV-2 S1 Protein Using Fast Methods

In this section, we predicted the effect of the detected mutations in ACE2 on the affinity of ACE2 variants to S1 by using different computational methods. There are many bioinformatic methods to predict the stability of the protein complex and the affinity between subunits by using various approaches. We utilized PISA,[42] FoldX,[43] Prodigy,[44] SEPAS,[45] and SAAMBE-3D[46] for predicting the intersubunit affinity of the ACE2/S1 complexes using their 3D structures.

In one category of structure-based affinity predictors, the utilized algorithm considers all parts of the interested 3D structure to predict the stability of protein complexes using the thermodynamic formulation of protein folding. PISA and FoldX are verified examples of the mentioned class of affinity predictors. In brief, PISA is an acceptable method to predict the biounit assembly of PDB-submitted structures by using an empirical energy function that considers the interaction energy and implicit representation of complex dissociation entropy to compute dissociation free energy of the complex. FoldX is a successful method in predicting the effect of mutations on protein complex stability. It performs the computations using semiempirical energy function that considers the stability of the monomers and the complex to predict complex stability.

Besides the introduced thermodynamic-based methods, some other structure-based methods use different structural aspects of the complex or interface region to predict intersubunit affinity. Prodigy, SEPAS, and SAMMB3D are examples of evaluated methods for this class of affinity/stability predictors. Prodigy is a successful method for predicting the experimentally measured binding affinities of protein subunits by counting the intersubunit contacts. SEPAS is a monomer-based predictor of binding affinity of protein binding patches to tentative partners that computes the mechanical stiffness of the proposed interface region. SAMMBE3D is a feature-based new generation of a trained version of adjusted GB/MM that is a predictor of the effect of mutations on dimer stability.

We report the predicted affinities of 240 mutated versions of ACE2 to S1 of SARS-CoV-2 using the mentioned fast structure-based computational methods in Figure . About 4.5% of the considered mutations reside in the ACE2–S1 interface region. Mutations in ACE2 that cause the ACE2/S1 complexes become more stable than wild type ACE2/S1 complex are red, whereas mutations that destabilize the ACE2/S1 complex are blue in Figure . If a method does not predict a considerable difference between the affinities of mutated and WT ACE2 for S1, the method is colored yellow for the considered mutation in Figure . The criteria for assigning different mutations outcomes are presented in methods sections –2.3.5. For three cases, three methods out of five predictors predict that the related mutations stabilize ACE2/S1 complex, and for two cases, 3/5 of the predictors expect the mutations destabilize the complex between ACE2 and S1. These methods are fast and utilize different approaches to perform the prediction.

Figure 1

Predicted affinity of the mutant ACE2 for the S1 protein is presented. The position of each ENSEMBL-reported SNPs for ACE2 protein is presented in the first row. The mutated residue and its corresponding WT residue are presented in third and second rows, respectively. Rows four to eight represent the effect of the ACE2 mutation on the stability of the ACE2/S1 complex predicted by different predictors. If a mutation is predicted as a stabilizer one, it is red, if the mutation leads to a more unstable complex than WT ACE2/S1 complex, it is blue, and if no change is predicted for stability, the color is yellow.

As the recruited fast methods in this section consider different structural aspects of protein complexes, the observed inconsistency in prediction seems natural. The utilized thermodynamic approaches are sensitive to the quality of the interested 3D structure. However, some of the other structure-based methods perform predictions by just utilizing general structural features that show less sensitivity to the global quality of the structure. The utilized methods in the presented part consider the static structure of the ACE2/S1 complex for predicting changes in affinity. Here, to alleviate this issue, we feed the SEPAS algorithm with ANM-generated ensembles of structures for predicting the affinity of each ACE2 mutant to the SARS-COV-2 S1 protein. The ANM-defined ensemble of structures lets us consider the effect of mutations on all regions of ACE2 better than using a single snapshot of the mutant structure.

Another possible source of the variation between the prediction results of thermodynamic-based and other descriptor-based affinity predictors is the interface issue. The correct 3D structure of the interface region is more important for the methods that are trained for the interface structural properties. Because the interface region between ACE2 variants and the RBD of S1 is derived by superimposing ACE2 of the WT complex and mutant ACE2, this approach possibly misses structural changes in the interface region resulting from mutations especially for mutations that do not occur at the interface region. Besides, we set restrict criteria to accept the influence of mutations on changes in the affinity between subunits. As the mentioned methods have different intrinsic errors, we assume if the difference between the predicted affinity of WT subunits and mutant ones is larger than 0.1 kcal/mol then it is a meaningful change in the affinity. If we accept a large difference in the predicted affinities between mutants and WT complexes, the amount of inconsistency will be decreased.

Despite the contrary results of the mentioned methods in predicting the effect of ACE2 mutations on their affinity for S1, the results indicate that those mutations may change the affinity of ACE2 to S1. It is a considerable result because if such predictions would be verified by other methods, it may mean the populations have different intrinsic susceptibility to COVID-19.

Adaptive Tempering Generated Ensembles, MM/GBSA Approach, and Dynamic Properties of ACE2/S1 Complexes

The utilized bioinformatic methods in previous section are fast with acceptable accuracy for prediction of the affinity between subunits but we need more accurate methods to predict the effect of ACE2 mutations on its affinity for S1 protein of SARS-CoV-2 that consider the dynamic features of protein structures. Because molecular dynamics (MD)-based approaches to compute the binding affinities are time-consuming, we decide to limit the target population for studying the consequence of ACE2 mutations on its affinity for S1 protein and computing the effects of mutations on complex dynamics by MD-based approaches.

We utilized MD-based approaches to predict the affinity between ACE2 and the S1 protein of SARS-CoV-2 for the mutations in the ACE2 gene that are reported in the Iranome database for ethnic groups of the Iranian population as SNPs.[40] We selected eight missense mutations of ACE2 from the Iranome project whose positions are available in the crystallized structure of WT ACE2. Four out of the eight SNPs are also reported in ENSEMBL for other populations (Table ).

Table 1

Properties of ACE2 Mutations in Iranian Ethnic Groups Are Presenteda

	dbSNP accession code	Iranian ethnic group	Polyphen2	Mutation Assessor Pred	MM/GBSA-based affinity (kcal/mol)
WT					–62 (9)
D225G	just in Iranome	Azeri	probably damaging	functional	–69 (8)
D494V	rs765152220	Persian Gulf Islander	probably damaging	functional	–61 (7)
F452V	just in Iranome	Turkmen	probably damaging	functional	–61 (8)
Q60R	rs759162332	Azeri	benign	nonfunctional	–63 (7)
Y199C	rs750145841	Turkmen	probably damaging	functional	–63 (7)
S331F	just in Iranome	Turkmen	probably damaging	functional	–62 (8)
T334M	just in Iranome	Azeri + Turkmen	benign	functional	–57 (8)
V485L	just in Iranome	Persian Gulf Islander	possibly damaging	functional	–63 (9)

Ethnic group, the predicted severity of the mutation, predicted functionality, and predicted affinity of the ACE2 variant to the S1 protein are presented in columns 2–5, respectively. The standard deviations of the predicted affinities are reported as parenthesized numbers in the last column. If the predicted affinity for a mutant in comparison to the WT shows a significant difference (p < 0.05), its value is in bold in the last column.

We used MM/GBSA and adaptive biasing force (ABF), a thermodynamic integration method, to compute the free energy of S1 binding to ACE2 for ACE2 mutants that are reported in Iranian ethnic groups.

The FoldX-generated 3D structures of ACE2 mutants in association with S1 were minimized and heated up gradually. To increase the sampling efficiency, we used an accelerated method, adaptive tempering (AT-MD), instead of common MD. AT-MD is a single copy version of the replica exchange method that finds structures with minimum energy faster than simple MD does.[55] After 2 × 106 steps of AT-MD simulation, we utilized a single-trajectory-based version of MM/GBSA[54,56] to make the first MD-based estimation of binding energy between varieties of ACE2 and the S1 protein of SARS-CoV-2 (Table ). As it is indicated in Table , for example, some members of the Azeri ethnic group that carry the DA225G mutation in their ACE2 have higher MM/GBSA-predicted affinity for the S1 protein of SARS-CoV-2 than members with WT ACE2 do. It may suggest that the carriers of the noted mutation may be predisposed to infection with the virus because their ACE2 protein shows high affinity for binding to the S1 protein.

The AT-MD-generated ensembles of the complexes provide us with rich pools of structures to study the structural properties of the ACE2/S1 complexes in the bonded state. AT-MD lets the ACE2/S1 complex sample the energy landscape with high efficiency. We report the AT-MD potential of the system, a measure of the ACE2 subunit structural stability in the presence of the S1 protein, as a function of the MM/GBSA-predicted affinity of ACE2 for S1 in 2D density plots (Figure S1).

We are interested to compute the structural stability of ACE2 mutants when they are in association with the RBD of the S1 protein. We propose the changes in structural stability of ACE2 will affect its binding site for S1. The structural stability of ACE2 mutants is much more informative when mutations reside in regions far from the binding site of the S1 protein.

In the utilized AT-MD, the ground state energy of the system is tuned in a history-dependent manner to pass some energy barriers ahead of the system in its energy landscape by thermostat-directed boosts. Because we do not want to unfold the subunits completely, we set the acceptable range of temperature between 300 and 320 K. In this range of temperature, all atoms RMSD of the considered structures are smaller than 5 Å. We detect faint two-state transitions in AT-MD simulations of the complexes (see Methods). The microtransition occurs in the WT complex or complexes of considered mutants within the temperature range of 308–309 K. The observed transitions are detected by considering the potential energies of complexes as a function of the simulation temperature. The transition energy between the folded state of the ACE2/S1 complexes and the locally melted state of the complexes are reported in Figure . We decomposed the amount of complex partial melting by computing the dRMS metric for ACE2 and RBD subunits during AT-MD simulation for pretransition population, i.e., populated stable species (Figure S2). To find which subunit of the considered dimer is affected more by temperature, we calculated the difference of dRMS between ACE2 and S1 in each simulation frame as ΔdRMS. Defining the ΔdRMS metric let us determine which subunit is affected more during AT-MD simulations. The 2D density plots of AT-MD potential energy as a function of ΔdRMS for WT ACE2 and mutant ACE2s in association with S1 are presented in Figure S3. We map 2D density plots presented in Figure S1 on their corresponding plots in Figure S3. It helps us to determine the ΔdRMS value that corresponds to the highly populated stable complex that has been used for determining the MM/GBSA-based affinity for the desired subunit pair. The obtained ΔdRMS values define the relative order of the ACE2 subunit structural stability when it is bonded to the S1 subunit. The summary of the procedure and the obtained relative order of ACE2 subunit structural stability are presented in Figure . These results indicate that some mutations stabilize ACE2 monomer and some decrease the structural stability of ACE2. In one subpopulation of Y199C mutant of ACE2, the structural stability of ACE2 is decreased, whereas the affinity of ACE2 for the S1 subunit is increased in that population. The ACE2 structural stability is decreased in Q60R mutation but its MM/GBSA-based affinity for S1 does not change significantly. As the studied mutations of ACE2 are far from the interface region, they modulate the affinity between ACE2 and S1 using long-range allosteric mechanisms. Although some of the mutations make ACE2 unstable in comparison to WT, they prepare the binding site region for making a stable complex with the SARS-CoV-2 S1 protein. These observations mean it is not necessary for an affinity modulator mutation to appear in the interface region of the complex for exerting its effect. In next steps, we follow how the considered mutations in ACE2 change its affinity to S1 over a distance.

Figure 2

Different states of the ACE2/S1 complex observed in AT-MD and the structural stabilities of ACE2 subunits of the stable complexes are described. The ACE2/S1 complexes show two populations in AT-MD simulations. The stable population, pretransition, transits to an unstable population after passing some barriers. The lower panel represents the ΔE for complexes with mutant ACE2 between the two faint states of complexes by fitting the computed data to the two-state model of protein transition (eq ). In the scheme that represents two-state transition, the x-axis represents 1/(kT) and the amount of melted structure is decreasing, while the metric presented in the vertical axis is close to 1. The shared highly populated states in Figures S1 and S3 show distinct amounts of ΔdRMS. We sort ACE2 subunits on the basis of the amount of ΔdRMS of their corresponding stable complexes (upper panel). The bold values of predicted affinities are meaningful at the p < 0.05 level for the t test compared to the affinity of WT assembly.

Considering variations of ACE2 and S1 protein electrostatic potential (EP) as averaged isopotential surfaces, we find that the large negative electrostatic potential of ACE2 makes the protein an interesting target for proteins with a positive accessible electrostatic potential like S1. The mutations in ACE2 affect the size and integrity of the negative EP in mutant versus WT ACE2 (Figure S4). The changes in EP for mutations at positions 225 and 452 are considerable. These mutations increase the affinity of ACE2 to S1 on the basis of MM/GBSA-based-predicted affinities (Table ).

The 3D structure of ACE2 indicates that there is a long groove in its structure. The noted groove is possibly developed for substrate binding. The active site of the enzyme also resides in one side of the groove. The accessibility of the groove to solvent is restricted by two regions assumed as “Gate”. The halves of the gate are interconnected in one side via a disordered region labeled as the “Zip” region (Figure ). We study the dynamics of ACE2 mutants and WT structures during AT-MD simulations by measuring the Euclidian distance between the halves of the gate region and the length of the zip structure. As indicated in Figure , the dynamics of the Gate and Zip regions are different between wild type ACE2 and the mutants of ACE2 in their S1-bonded states. A mutation at the position 452 causes an increment of ACE2 affinity to S1 and simultaneously closes the Gate region. Also, the Zip segment becomes shorter there. The gate is more opened in the T334M mutant of ACE2 in comparison to the WT subunit, whereas the affinity of this mutant of ACE2 for the S1 protein is decreased. It seems that the closed/open state of the ACE2 changes its affinity for the S1 protein of SARS-CoV-2. This observation and the mentioned hint may be an opportunity for drug discovery and finding a possible therapy against COVID-19.

Figure 3

Dynamics of the ACE2 groove region is measured for AT-MD-generated ensembles of ACE2/S1 complexes. Left panel depicts the schematic position of the Gate and Zip regions in the ACE2 structure. The length of the zip region and the distance between halves of the gate region are presented in middle and right panels, respectively, for different ACE2 variants in association with the S1 protein. The horizontal axes represent length in angstroms, and the vertical axes present population density.

The perturbation-response scanning (PRS) approach provides the possibility to mechanically perturb specific residues and then predict its effects on the dynamics of other regions of the protein.[72,73] Because we believe that the considered mutations of ACE2 for Iranian ethnic groups exert their effects on the binding affinity between ACE2 and S1 over a distance, we recruit the PRS approach for ACE2/S1 complexes that contain the WT or a mutated version of ACE2.

Anisotropic network model (ANM) is a coarse-grained description of a protein in which the nonbonding interactions between residues are represented by linear-elastic springs. This model lets us evaluate the effects of an external-force-driven mechanical perturbation of one node, residue, on other nodes of the network and vice versa. Here, we perturb all residues and measure the sensitivity of the mutation positions to the perturbation using the PRS method. It means how much the mutated ACE2 residues sense the overall mechanical perturbation in an ACE2/S1 complex.

The presented amounts of the residue sensitivity in vertical axis of graphs in Figure indicate that the sensitivity of S1 in S1/ACE2 complex is high in response to ACE2 perturbations. It means that S1 is affected by ACE2 dynamics more than the amount ACE2 is influenced by the S1 perturbation.

Figure 4

Average sensitivity of residues or the sensitivity of a specific residue to the perturbation of other residues is presented. In the “All” panel, the average sensitivity of all protein residues to the physical perturbation is presented. In other panels, the sensitivity of a specific residue (panel name) to the perturbation of other residues of the complex is presented. The positions of interfacial residues are presented as black bars parallel to the x-axis of the “All” panel. Here, chain A represents the ACE2 subunit, and chain E denotes the S1 subunit. The surface representation of ACE2/S1 assembly in panel “225” is colored in the following manner: the ACE2 chain is blue; the S1 subunit is red; distal residues of S1 are depicted by a gray surface; and S1 residues that reside in the interface region of S1 and ACE2 are shown as a yellow surface.

The average sensitivity profile of each residue of the complex in response to perturbation of other residues of the complex is presented in panel “All” in Figure . Peaks in Figure indicate which residues are strongly coupled with the sensor, mutated, residues. As it is presented in panel “All” of Figure , some of the S1 distal residues (resid: 334, 360–369, 388, 427–430, 518) are highly coupled with the other residues in the complex. This region of S1 is close to recently discovered allosteric modulator region of the S1 protein.[74] Residues at positions 225, 452, and 485 of ACE2 protein are highly coupled with the distal residues of the S1 protein dynamically. Therefore, the mentioned positions of ACE2 play a critical role in information transfer between ACE2 and S1 subunits of the complex. Position 60 of ACE2 shows an odd pattern; this node has no significant coupling with other residues of the complex. Among the ACE2 residue positions that we studied, corresponding to the mutated residues in Iranian ethnic groups, just positions 199 and 225, are strongly coupled with residues of the S1 subunit that are in contact with the ACE2 subunit; so, their perturbations may affect the assembly more than other positions. This analysis indicates that position 199 in ACE2 is also coupled with some residues of ACE2; therefore, a mutation in that position possibly affects the ACE2 structure and the ACE2/S1 interface simultaneously. Position 225 is just sensitive to the signal arrived from the S1 contact site with the ACE2 subunit. This observation may rationalize the higher affinity of ACE2 to S1 resulting from the D225G mutation in ACE2. In brief, perturbations in positions 225, 452, and 485 of ACE2 in complex with S1 provoke a high amount of perturbation in S1; this perturbation response is especially observed for position 225 that is coupled with the interface region of S1. Perturbations in some positions such as 199 induce dynamic changes in both partners of the complex. The PRS-predicted dynamics changes in the interface region of S1 and ACE2 may be another justification for different affinities of ACE2 mutants to S1 of SARS-CoV-2.

Dissociation Process of ACE2/S1 Complexes

The structural aspects of the complex derived from AT-MD simulations and MM/GBSA-based predicted affinity of WT and mutant ACE2s for S1 provide us with some information about the possible mechanisms behind the changes in the affinity of ACE2 to S1. We utilize another MD-based method to improve the accuracy of the predicted affinity of ACE2 for S1 protein and to simultaneously study the dynamics of the partners during the dissociation/association process. Using the adaptive biasing force (ABF) method,[57] we derive the potential of mean force (PMF) for dissociation of S1 from ACE2 in complexes that contain the WT or a mutant version of ACE2. Before performing ABF simulations, the WT and mutant ACE2 monomers in association with S1 protein pass 2 × 106 steps in AT-MD. Next, the S1 subunit is dissociated from ACE2 slowly along the z-axis of the complex (Supplementary Movie). Considering the standard errors of the ABF-computed dissociation free energies (Figure S5), we find that the ABF-computed affinities between S1 and ACE2 are increased in many of ACE2 mutations in Iranian ethnic groups in comparison to WT ACE2 (Figure ). It means that for those mutations we need higher forces to dissociate the S1 subunit from the mutant ACE2 structures than for WT ACE2. Among the considered mutations, the V485L mutant shows a lower affinity for binding to the SARS-CoV-2 S1 protein than WT shows. It may provide an intrinsic resistance against COVID-19 to its carriers.

Figure 5

ABF-computed dissociation free energies are presented. The relative free energy of dissociation of S1 from ACE2, which is computed by the ABF approach, is presented for different variants of ACE2 in Iranian ethnic groups. The computed free energy is reported as a relative value compared with the computed affinity for WT complex, 48 kcal/mol. The horizontal axis represents the relative distance between S1 and ACE2. A higher value in the x-axis denotes a higher distance between two subunits of ACE2/S1 complex along the axis interconnecting two subunits. Zero in the x-axis represents the bound state, and 1 represents the state in which S1 is dissociated from ACE2 completely.

The presented results in section suggest that the dynamics of the gate/zip regions of ACE2 possibly affects the affinity of ACE2 for S1. To study the effect of the ACE2 closed state on its affinity for S1, we perform ABF simulations also for ACE2 bonded to the ORE-1001 ligand as an investigational drug that converts the open state of ACE2 to its closed state. Studies suggested that during binding of substrate/inhibitor to the ACE2 groove it transits to the closed state.[16] Our computations indicate that the affinity between S1 protein of SARS-CoV-2 and ACE2 in the ORE-1001-induced closed state is higher than that in the open conformation (Figure ). We should note that the drug also changes the dynamics of dissociation process; therefore, we would not extrapolate the same conclusion for all closed states of ACE2.

Nowadays, there are many reports about the effects of high blood pressure (HBP) suppressor drugs, ACEi/ARBs, on the incidence of COVID-19. There are some contradictory results in such reports about the effect of such drugs on predisposition to COVID-19.[75,76] Many of the HPB suppressors affect ACE2 gene expression in a tissue-specific manner.[77−79] To the best of our knowledge, no experimental result has been published about the effect of ACE2-specific drugs on the affinity of ACE2 for S1 of SARS-CoV-2. Some studies reported the effect of binding of some molecules to ACE2 in the closed conformation on the affinity of ACE2 for the SARS spike.[80,81] They proposed such drug–ACE2 interactions may decrease ACE2 affinity for the SARS spike. They also postulated the critical role of the closed, substrate-bonded state of ACE2 in its binding to SARS or the COVID-19 virus spike protein. Recently, researchers reported that the activity of ACE2 increased up to 10 fold by binding to RBD of the SARS-CoV-2 S1 protein.[82] The authors of the mentioned paper suggested that the RBD of the SARS-CoV-2 S1 protein binds to the closed state of ACE2. We report the dynamics of Gate and Zip regions of ACE2 for some ACE2 mutations (Figure ). Our computations suggest that the closed and open states of ACE2 show a different affinity for the RBD of the SARS-CoV-2 S1 protein (Figure ). When ACE2 binds to ORE-1001, an ACE2-selective inhibitor, it transits to the closed state. In that state, ACE2 shows a higher affinity for the SARS-CoV-2 S1 protein than in the open state. These observations suggest that ACE2 possibly fluctuates between closed and open states and the RBD of the SARS-COV-2 S1 protein has a higher affinity toward the closed population of the ACE2 protein. Therefore, it may mean that if some compounds induce the closed state in ACE2, S1 may bind to ACE2 with higher affinity.

The dissociation of S1 from ACE2 in ABF simulation sets indicates that the diverse ACE2 mutants dissociate from S1 differently. In some cases, aromatic–aromatic contacts between the partners are disrupted quickly (mutations at 331, 199, 452, and 485; Figure S6). Considering the salt bridges between the subunits of ACE2/S1 complexes, we find that the bridge between residues 31 and 484 of ACE2 and S1, respectively (bridge 31–484) and the bridge between residue 329 of ACE2 and residue 439 of S1 (bridge 329–439) are two common bridges in ACE2/S1 assembly (Figure ), which were described in the crystal structure and simulations.[19,83] The second salt bridge was introduced unnaturally to stabilize the complex during crystallization steps. The structural changes resulted from the F452V and S331V mutations in the ACE2 force ACE2 to create a new short-life salt bridge between residue 23 of ACE2 and 458 of S1 (bridge 23–458). The intersubunit salt bridges and aromatic–aromatic contacts indicate that the disruption of the 329–439 bridge is concomitant with tearing aromatic contacts between S1 and ACE2. It suggests that the bridge acts as a lock for disruption of intersubunit hydrophobic contacts.

Figure 6

Intersubunit salt bridges along ABF simulations are presented. Three types of salt bridges detected between the S1 protein and different variants of ACE2. The X-axis represents the normalized simulation time along ABF simulations. Zero in the x-axis represents the bound state, and 1 represents the state in which S1 is dissociated from ACE2 completely. The vertical axis represents the total number of inter-residue contacts if the involved residues are closer than 4 Å.

In ABF-derived results, we considered the effect of eight mutations of ACE2 observed in Iranian ethnic groups on the ACE2 affinity for the S1 protein of SARS-CoV-2. We measured the distances between positions of the eight mutated residues in ACE2 structures and the interface region of ACE2 along dissociation simulations (Figure S7). These eight mutations are distributed in different regions of ACE2, and their distances from the interface region vary in different mutations. It means that not only the structure of the interface region is affected during dissociation process but also other regions of ACE2 that are far from the interface are also affected. These observations suggest that the integrity of the ACE2 structure in different mutants control the amount of forces required for dissociating S1 from ACE2. Changes in the integrity of mutants and WT ACE2 structures during the dissociation process of S1 are dynamic. To depict them, we follow all salt bridges that appear in ACE2 subunits of the complexes during the dissociation of S1 from ACE2. Then, we detect similar salt bridges, with the same donor and acceptor of the bridge, in different varieties of ACE2 along ABF simulations. We compute the similarity between ACE2 varieties on the basis of the appearance of the same salt bridges along ABF experiments (Figure ). For example, the salt bridges that appear in the Q60R ACE2 mutant are unique ones and show a low similarity to salt bridges of other ACE2 mutants in the S1/ACE2 dissociation experiment.

Figure 7

Diversity of salt bridges appearing in the ACE2 subunit during ABF simulations is compared among variants of ACE2. The unique intra-ACE2 salt bridges are detected during ABF simulations. The color represents how many of salt bridges in the ACE2 subunit are similar between different types of ACE2 along simulation time. 0–10%: light-gray, 10–20%: gray, 20–30%: dark blue, 30–40%: light blue, 40–50%: orange, 50–60%: red.

Studying the changes in VDW and electrostatic terms of interaction energy between S1 and ACE2 along the dissociation process indicates that the electrostatic portion of the interaction energy between subunits is lower for the V485L mutant than for WT (Figure S8). This is a mutation with the lowest ABF-computed affinity for S1 in the current study. On the other hand, for the most stable ACE2/S1 complex in the current study, the F452V ACE2 mutant, the VDW interaction energy diminishes quickly.

On the basis of our studies on the process of S1 dissociation from ACE2, we might divide the receptor binding motif (RBM) of the S1 protein into two distinct sections “Fist” and “Forearm”. These parts have different affinities to the ACE2. The ratio of hydrophobic (or polar) contacts between fist and ACE2 and hydrophobic (or polar) contacts between forehand and ACE2 is computed along dissociation of S1 from ACE2 (Figure ). We find that the ratios change differentially among ACE2 mutants along dissociation steps. The ratios indicate that the first segment of S1 that dissociates from ACE2 is forehand, and in final steps, when S1 protein struggles to dissociate completely from ACE2, the fist segment (residues 470–490 in the S1 protein) will be detached from ACE2. If we suppose in ideal state that the dissociation process is similar to the association process but in reverse direction, we may conclude that in the first stages of its binding, the S1 protein docks its fist section into ACE2. This step lets the RBD of S1 find the appropriate region on ACE2 for its complete binding via the forehand segment.

Figure 8

Superiority of different parts of the S1 subunit in interaction with ACE2 is traced along the dissociation simulations. The proposed segments of the RBD of the S1 protein in the current work are depicted in the right panel. The normalized number of contacts between hydrophobic residues of ACE2 and the fist region of S1 divided by the sum of hydrophobic contacts between ACE2 and fist and forearm regions of S1 is represented as blue lines. Red lines represent polar contacts between ACE2 and S1 segments. The vertical axis represents the normalized number of contacts. The X-axes represent the normalized time along dissociation experiments in ABF simulations. The contact criterion is set at <10 Å for the interatom distance.

The dynamics of the denoted segments of S1 is affected by dynamic changes of intra-ACE2 interactions which may be translated as changes in ACE2 stability along the dissociation process. Considering the changes in folding stability of ACE2 subunit during the dissociation of S1 from ACE2, we find that the structure of ACE2 becomes more stable in some mutants (Figure ). The increased stability of ACE2 in the F452V mutant even causes a resistance against dissociation of S1 from ACE2. The FoldX-defined folding stability analysis also indicates that the electrostatic portion of ACE2 stability is decreased in most ACE2 mutants. It possibly indicates that the electrostatic interactions in the considered mutants are not in the optimal state.

Figure 9

Folding stability of the ACE2 subunit of the complex is computed along the dissociation process. While the S1 subunit is slowly dissociated from the ACE2 subunit, the stability of ACE2 is computed for every frame of the ABF simulation (middle panel) by the FoldX approach. Left and right panels represent the density plots for the amount of electrostatic and vdw terms of the stability function for variants of ACE2, respectively. Error bars represent the standard error of computed values. For each variant, 100 samples are picked from AT-MDs to evaluate the stability.

Conclusion

In the current study, we report the effect of ACE2 polymorphism on its binding affinity for the S1 protein of SARS-CoV-2, especially for the reported mutations in ACE2 of Iranian ethnic groups as mutations that reside far from the S1 binding site. The dynamics and stability of the ACE2 mutants are studied in our work. We demonstrate that the considered mutations affect the affinity between ACE2 and S1 via long-range mechanisms. Many of the considered mutations in the current study enhance the affinity of ACE2 to S1 of SARS-CoV-2, and it may suggest possible intrinsic susceptibility of carriers of such mutations to COVID-19. On the other hand, we find one of ACE2 mutants has a lower affinity for the S1 protein of SARS-CoV-2 than WT does, which possibly suggests a type of immunity for its carriers. Our computations suggest that the affinity of ACE2 for S1 in its closed state may be greater than that in the open state of ACE2. We report the relation between dynamics of ACE2 structure and its affinity to the S1 protein of SARS-CoV-2. This information may be informative to design compounds that modulate ACE2 dynamics and consequently decrease its affinity for the SARS-CoV-2 S1 protein.