key: cord-0767590-9u88qusm
authors: Barros, Emilia P.; Casalino, Lorenzo; Gaieb, Zied; Dommer, Abigail C.; Wang, Yuzhang; Fallon, Lucy; Raguette, Lauren; Belfon, Kellon; Simmerling, Carlos; Amaro, Rommie E.
title: The flexibility of ACE2 in the context of SARS-CoV-2 infection
date: 2020-11-13
journal: Biophys J
DOI: 10.1016/j.bpj.2020.10.036
sha: 7126d80ee882752303ac2fb184136becadb859ea
doc_id: 767590
cord_uid: 9u88qusm

The COVID-19 pandemic has swept over the world in the past months, causing significant loss of life and consequences to human health. Although numerous drug and vaccine development efforts are underway, there are many outstanding questions on the mechanism of SARS-CoV-2 viral association to angiotensin-converting enzyme 2 (ACE2), its main host receptor, and host cell entry. Structural and biophysical studies indicate some degree of flexibility in the viral extracellular spike glycoprotein and at the receptor binding domain-receptor interface, suggesting a role in infection. Here, we perform explicitly solvated all-atom molecular dynamics simulations of the glycosylated, full-length membrane-bound ACE2 receptor, in both an apo and spike receptor binding domain (RBD) bound state, in order to probe the intrinsic dynamics of the ACE2 receptor in the context of the cell surface. A large degree of fluctuation in the full length structure is observed, indicating hinge bending motions at the linker region connecting the head to the transmembrane helix, while still not disrupting the ACE2 homodimer or ACE2-RBD interfaces. This flexibility translates into an ensemble of ACE2 homodimer conformations that could sterically accommodate binding of the spike trimer to more than one ACE2 homodimer, and suggests a mechanical contribution of the host receptor towards the large spike conformational changes required for cell fusion. This work presents further structural and functional insights into the role of ACE2 in viral infection that can potentially be exploited for the rational design of effective SARS-CoV-2 therapeutics.

Angiotensin-converting enzyme 2 (ACE2) acts as the extracellular receptor for the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) (1-3), the virus responsible for the COVID-19 pandemic that has catastrophically affected the world since its first identification in December 2019 (4-7). ACE2 is a membrane protein found in lungs, kidneys, heart and intestine cells (8, 9) that plays a physiological role in cardiovascular regulation via the cleaving of intermediates in the maturation process of angiotensin, a peptide hormone involved in vasoconstriction control (10) (11) (12) (13) (14) . ACE2 is a homodimer with a large claw-like extracellular head domain, a small transmembrane domain and a short intracellular segment (8). The head can be further subdivided into the catalytic zinc-binding peptidase domain (PD, residues 19 to 615) (15) , and the smaller neck domain (residues 616 to 726), which is where the majority of the homodimer interactions seems to lie (16) . The neck domain is further connected to the singlehelix transmembrane (TM) domain by a long linker (Figure 1a ). ACE2 can also function as a membrane-trafficking chaperone for B 0 AT1, an amino acid transporter (17) , and it was in fact only in complex with this partner that the single TM helix of ACE2 could be resolved (16) .

SARS-CoV (18) (responsible for 8096 cases worldwide in 2002 (19) ) and now the closely related SARS-CoV-2 hijack ACE2 as the host cell receptor to its large extracellular spike (S) glycoprotein (1, 20) . The spike's receptor binding domain (RBD) in the "up" conformation binds to ACE2's PD with high affinity (21) , and the resolved ACE2-RBD complex consists of a dimer of heterodimers, with each monomer of the ACE2 homodimer interacting with one RBD (thus forming a heterodimer, Figure 1a ). Cellular recognition and binding to ACE2's peptidase domain via the RBD is proposed to initiate a series of complex conformational transitions in the S homotrimeric protein, leading to the shedding of its S1 subunit and fusion to the host cell membrane driven by the S2 subunit (22) (23) (24) (25) (26) , ultimately resulting in the infection of the host cell.

Downregulation of ACE2 and accumulation of angiotensin II due to spike binding is also associated with acute respiratory distress syndrome (ARDS) and acute lung failure (27-31), contributing to SARS-associated symptoms. As such, the S glycoprotein and ACE2-S complex are considered key targets for drug and antibody development efforts aiming to curtail the virus' remarkable transmissibility and negative effect on human health (1, (32) (33) (34) (35) , including exploiting the ACE2-S high affinity with recombinant soluble ACE2-antibody constructs (35) (36) (37) (38) .

Experimental and biophysical studies of the SARS-CoV-2 RBD and soluble ACE2's PD complex have suggested structural factors likely responsible for the higher affinity and infectivity of SARS-CoV-2 compared to SARS-CoV (39) (40) (41) , and revealed significant dynamics at the RBD-PD interface in the form of rocking motions between the two subunits (41, 42) . At the individual PD level, opening and closing of the active site cleft has been observed in X-ray structures of the extracellular region (43) . Additionally, the recent cryoEM structures of fulllength ACE2 indicate that the homodimer can adopt a less populated open conformation, as defined by the distance between the head domains, in addition to the closed conformation shown in Figure 1a CryoEM and computational studies of the full-length spike glycoprotein have recently suggested a significant degree of flexibility of the spike's stalk and at the ACE2-RBD interface (44, 45), evidencing the need to study these macromolecular complexes in the context of the cell surface.

Here, we perform all-atom molecular dynamics (MD) simulations of the full length, membraneembedded and glycosylated ACE2 homodimer both in the apo state and in complex with RBD to study the dynamics and molecular origins of the ACE2-S flexibility on the host receptor side that evade experimental characterization (46) . Seven complex N-glycans and one O-glycan in ACE2

were modeled according to glycoanalytic data (47) (48) (49) , as well as glycan N343 in RBD for the RBD-bound simulations (Figure 1b and Supplementary Table 1 ). The B 0 AT1 transporter solved in the cryoEM structure was not included in our simulations in order to computationally probe the intrinsic dynamics of ACE2. B 0 AT1 is mainly expressed in kidneys and intestines (50) , whereas ACE2 can also be found in lungs and heart tissues, supporting the likelihood that ACE2 can be found un-complexed with B 0 AT1 upon cellular recognition and binding to S.

Distributions and average values were calculated considering the total sampling in the three replicas for each system. Except when the dimer structure or inter-dimer contacts were analyzed, the monomer sampling was accumulated, such that the calculation correspond to 6 µs of combined sampling (3 replicas x 2 monomers/replica). Reported uncertainties of average values correspond to the standard deviation across the replicas.

To quantify the range of motion of ACE2 in the simulations, several angle and distance metrics were developed. Calculation was performed using MDTraj (64) with visualization through VMD (65) . The 6M17 cryoEM structure was used as the reference structure.

Head tilt angle relative to the transmembrane domain was calculated by first aligning the dimer's coordinates to the reference cryoEM TM domains, and angle calculated between the centers of mass of the reference's dimer peptidase domains (residues 18 to 600), reference's TM helices (residues 747 to 774), and monomer's PD at each frame in the simulation. Helix tilt angle was computed as the angle between a vector defining the membrane's normal and a vector connecting residues 741 and 765 at the extremities of the helix.

Revolution angle was calculated between the center of mass of the monomer's PD in the reference conformation, the center of mass of the reference's TM domain, and the center of mass of the monomer's PD at each frame in the simulation following alignment of the monomer's TM helices. Buckling angle was calculated using the xy projections of the center of mass of monomer's A PD at frame f, the center of mass of the reference dimer's PDs, and center of mass of monomer's B PD at frame f, after alignment of the trajectories to the reference dimer neck domains (residues 617 to 726).

Distance between the monomer's head domains in the homodimer was calculated by determining the distance between each monomer's PD center of mass. Distance between the head domain and membrane corresponds to the minimum distance between the PD's heavy atoms and membrane's phosphorous atoms at each frame of the simulation. Distance between TM helices was calculated based on the distance between their centers of mass.

J o u r n a l P r e -p r o o f Principal Component Analysis. PCA was performed using the Scikit-learn library in python 3.6 (66). The monomers in each of the ACE2 homodimer trajectories were saved separately and all aligned to the backbone of the transmembrane helix, residues 747 to 774. The simulation coordinates of the apo and RBD-bound systems were concatenated and used to fit the transformation function, so that both systems were transformed in the same PC space.

RMSF and secondary structure analysis. Root-mean-square fluctuation (RMSF) of Cα carbons was calculated using CPPTRAJ (67) . Calculations were performed for the head and TM domains separately, following alignment of the backbone atoms of the respective domain.

Secondary structure calculation was performed with MDTraj's compute_dssp function (64) using the simplified 3-category assignment scheme. The proportion of frames in which each residue's secondary structure was assigned as a coil, helix or strand was computed and compared to the assignment obtained for the starting cryoEM structure.

Fraction of native contacts and glycan contacts. Fraction of native contacts was calculated according to Mehdipour and Hummer (68) . The 6M17 cryoEM structure was used as the reference structure for identification of native contacts. The ACE2-RBD interface was subdivided into three interacting regions according to the interacting residues pairs listed on Supplementary Table 2. A systematic characterization of contacts established by each glycan in the system was performed using MDTraj (64), using a cutoff of 3.5 Å between the heavy atoms.

The spike model was obtained from our previous simulations (69) . The simulations included only the solvated spike. All atoms except for the spike protein and glycans were removed, along with the lower part of the stalk region of each protomer (residues 1165-1273). Since the cryoEM model was missing density for portions of the RBD, we replaced the RBD coordinates (residues 355 to 494 for closed RBD, 339 to 523 for open RBD) of each protomer in both models with the RBD coordinates from the crystal structure of the RBD bound to ACE2 (PDB 6M0J (70)). The RBD structure from 6M0J was aligned with the backbone heavy atoms (alpha-carbon, carbonyl-carbon, and nitrogen) of each RBD in the initial spike model. We J o u r n a l P r e -p r o o f then grafted the RBD coordinates onto the spike at the hinge region , which resolved missing loops as well as introducing a disulfide bond in the RBD. The remaining disulfides not resolved in the cryoEM structures were assessed based on distance criteria and sequence conservation.

The system was built using the ff14SBonlysc (71) and GLYCAM (72) force fields for the protein and glycan atoms, respectively. These were explicitly solvated in OPC3 water (73) with a 200 mM NaCl buffer (74) . The RBD-up and -down systems both consisted of 1,298,646 atoms, and were simulated on Frontera at TACC, and SDCC at BNL using the pmemd.CUDA module of Amber20. The spike systems were equilibrated using a 10-step protocol. First, the water molecules were minimized for 1,000 steps using steepest descent, and then for an additional 9,000 steps with conjugate gradient while the rest of the system was positionally restrained with 1 kcal/( mol * Å 2 ) restraints. The systems were then heated to 310 K at constant volume, again with all atoms except hydrogens and waters restrained with 100 kcal/(mol * Å 2 ) positional restraints. The box size and density were then equilibrated over 1 ns with constant pressure, with the same positional restraints as the previous step. The restraints were then lowered to 10 kcal/(mol* Å) for an additional 1 ns of equilibration, before a second minimization. This minimization consisted of 10,000 steps of conjugate gradient with positional restraints now only applied to backbone atoms (alpha-carbon, carbonyl-carbon, and nitrogen), using a force constant of 10 kcal/(mol * Å 2 ). The next three steps of equilibration were MD for 1 ns each at constant NPT with positional restraints on protein backbone atoms at 10, 1, and 0.1 kcal/(mol * Å 2 ), respectively. This was followed by a final 1 ns of unrestrained MD at constant NPT before beginning production.

To generate a 3-up model of the prefusion spike protein, steered molecular dynamics (SMD) was used as implemented in the Amber NFE toolkit (75) . The initial structure was an allclosed model from the equilibration described above. To generate a structure used as reference 

Additional supporting research data of our models to enable other groups to use and explore this dynamic system in atomic detail (77) may be accessed through the NSF MolSSI and BioExcel COVID-19 Molecular Structure and Therapeutics site at https://covid.molssi.org/. Figure 5a) , which will be discussed in more detail in the following sections. Despite the global motions and backbone fluctuations, the protein secondary structure is overall stable, retaining the motifs identified in the cryoEM structure ( Supplementary Figures 6 and 7) .

The simulations indicate that the conformational variability of ACE2 occurs not only due to flexibility at the linker connecting the transmembrane and head domains, but also due to motions of the transmembrane helix in the membrane. In contrast to other multimeric transmembrane domains such as the coiled coil trimer of S (21, 26, 69) , each ACE2 monomer is anchored to the membrane by a single helix, which does not interact with that of the opposing monomer ( Figure 4a ), but rather explores a range of tilt angles relative to the membrane's normal ( Figure 4b ). This tilt motion, however, does not affect the overall integrity of the transmembrane helix, as indicated by RMSF and secondary structure characterization of the TM-aligned conformations (Supplementary Figures 5b, 6c and 7c) . (41, 42) .

In addition to protein interactions, five glycans in ACE2 are in close proximity to the RBD and have been suggested to play a role in S binding. In agreement with other studies (48, 68) , N90 and to a lesser extent N322 of ACE2 establish contacts with RBD. Besides these glycans, we also find that N53 can form a large number of contacts with the RBD residues, while the RBD glycan N343 makes very few contacts with ACE2's head protein residues (Figure 5c) or glycans (Supplementary Figure 8) .

Besides the ACE2-RBD heterodimer interface, we considered the interactions within ACE2 that could contribute to maintaining such a stable homodimer head interface despite the pronounced flexibility of the protein. Experimental structures and simulations suggest that the majority of the protein contacts in the homodimer are located in the neck domain, with only two other interactions, in the form of hydrogen bonds, are observed in the larger peptidase domain (16, 68) . In agreement with these observations, we find that the dimer interface is mainly held together in the simulations by residues at the neck (Figure 6a for RBD-bound simulations and Supplementary Figure 9a for apo) . Computation of the glycan-protein and glycan-glycan contacts enrich the characterization of the inter-and intra-monomer interactions and suggest that, while the eight ACE2 glycans form several contacts with protein residues located within the same monomer (Figure 6b) , protein-glycan interactions with the opposite monomer are limited to N53 and N690 (Figure 6c) . Additionally, N53 is the only glycan that can be found to form intermonomer glycan-glycan contacts, established between the equivalent N53 copies in each of the J o u r n a l P r e -p r o o f monomers (Figure 6d and Supplementary Figure 10) . In this sense, it is interesting that computational predictions indicate that disruption of the N53 glycosylation motif due to mutation at T55, leading to removal of the glycan, can have a destabilizing effect on ACE2 stability (78) .

Similarly to the lack of RBD effect on the ACE2 dimer flexibility, we find that the homodimer contact distributions are comparable between apo and RBD-bound states of ACE2 (Supplementary Figure 9) .

Our systematic analysis of all glycan interactions thus indicate that N53 is the only glycan involved in both heterodimer (ACE2-RBD, Figure 5c ) and homodimer (intra ACE2 dimer, Figure 6d Thus, N53 can still be found highly solvent exposed in the apo state, suggesting optimal conformations for contact with an RBD partner and a role in S binding and infection.

All-atom simulations of apo and RBD-bound, full-length, membrane embedded ACE2 In accordance with our findings, a high deformation propensity was also observed for the TM-neck linker upon normal mode analysis of full-length ACE2 (79) . These distinct conformations were likely only observed in the simulations due to the removal of B 0 AT1 from the co-complexed structure, since they seem to bind on opposite sides of the homodimer transmembrane interface and interact with the flexible linker (16) . The presence of B 0 AT1 would thus likely reduce the observed ACE2 flexibility, restraining the sampling of highly tilted conformations. However, the expression profiles of ACE2 and B 0 AT1 suggests the likelihood of ACE2 existing in the apo state, especially in lung and heart tissues where B 0 AT1 is not expressed (50) . The observed high degree of ACE2 flexibility in the apo form can additionally explain the difficulty in the experimental characterization of the full length structure.

Remarkably, the conformations of the ACE2 peptidase and neck domains remain stable throughout the simulations, and the homodimer heads retain their relative orientation despite the dramatic global homodimer motions. In a similar fashion, the RBDs included in the holo simulations remained tightly bound to ACE2 throughout the simulations, evidencing the high affinity between them. Glycan-glycan and glycan-protein interactions suggest that the ACE2 homodimer interface is maintained not only via protein interactions at the neck domain, but also inter-monomer contacts involving N53, at the top of the PD, and N690, closer to the neck.

Interestingly, N53 also makes extensive contacts with RBD, suggesting a dual and possibly competing role between homodimer (intra-ACE2) and heterodimer (ACE2-RBD) interactions.

This dual nature may be dependent on the length of the N53 glycan, but the glycosylation heterogeneity in ACE2 in general and at this position in particular (49) Additionally, it has been found from in-situ cryoEM and molecular dynamics simulations that the spike glycoprotein also exhibits conformational plasticity, with hinge motions at three different regions of the stalk trimer (44). Large dynamical variations thus seems to be a feature of these extracellular glycoproteins. The RBD rocking motion and S conformational variability have been proposed as mechanisms for immune evasion and efficient receptor search in the host cell (44, 45) but the similar rocking motion of ACE2 we observed also suggests a mechanical aspect to ACE2-S interaction. The process of S conformational transition upon binding to the receptor and cell fusion remains elusive, but ACE2's intrinsic flexibility could promote a large swinging motion of the ACE2-S1 complex, providing a mechanical force for the approximation of the two membranes and shedding of S1 towards fusion of the S2 domains into the receptor cell ( Figure 7b ).

Finally, one can speculate that the flexibility of the host receptor might allow the accommodation of more than one ACE2 dimer bound to a single S with two or more RBDs in 

SARS-CoV-2 Cell Entry Depends on ACE2 and TMPRSS2 and Is Blocked by a Clinically Proven Protease Inhibitor

A pneumonia outbreak associated with a new coronavirus of probable bat origin

Structural and Functional Basis of SARS-CoV-2 Entry by Using Human ACE2

Clinical features of patients infected with 2019 novel coronavirus in Wuhan

A novel coronavirus outbreak of global health concern

Digestive system is a potential route of COVID-19: An analysis of single-cell coexpression pattern of key proteins in viral entry process

ACE2: From vasopeptidase to SARS virus receptor

-7)-forming activity in failing human heart ventricles: Evidence for upregulation of the angiotensin-converting enzyme homologue ACE2

ACE2: A new target for cardiovascular disease therapeutics

Angiotensin-converting enzyme 2 is an essential regulator of heart function

Angiotensin-converting enzyme 2 is a key modulator of the renin-angiotensin system in cardiovascular and renal disease

Structural biology: Structure of SARS coronavirus spike receptor-binding domain complexed with receptor. Science (80-. )

Structural basis for the recognition of the SARS-CoV-2 by full-length human ACE2. Science (80-. )

A protein complex in the brush-border membrane explains a Hartnup disorder allele

People's Republic of China

SARS and MERS: Recent insights into emerging coronaviruses

Angiotensinconverting enzyme 2 is a functional receptor for SARS coronavirus

Characterization of severe acute respiratory syndrome-associated coronavirus (SARS-CoV) spike glycoprotein-mediated viral entry

Activation of the SARS coronavirus spike protein via sequential proteolytic cleavage at two distinct sites

Proteolytic activation of the SARS-coronavirus spike protein: Cutting enzymes at the cutting edge of antiviral research

Cryo-EM structure of the SARS coronavirus spike glycoprotein in complex with its host cell receptor ACE2

Angiotensin-converting enzyme 2 protects from severe acute lung failure

Trilogy of ACE2: A peptidase in the renin-angiotensin system, a SARS receptor, and a partner for amino acid transporters

Differential Downregulation of ACE2 by the Spike Proteins of Severe Acute Respiratory Syndrome Coronavirus and Human Coronavirus NL63

Clinical and biochemical indexes from 2019-nCoV infected patients linked to viral loads and lung injury

Structure, Function, and Evolution of Coronavirus Spike Proteins

Exploitation of glycosylation in enveloped virus pathobiology

Inhibition of SARS-CoV-2 Infections in Engineered Human Tissues Using Clinical-Grade Soluble Human ACE2

Neutralization of SARS-CoV-2 spike pseudotyped virus by recombinant ACE2-Ig

Engineering human ACE2 to optimize binding to the spike protein of SARS coronavirus 2. Science (80-. )

Engineered ACE2 receptor traps potently neutralize SARS-CoV-2

Is the Rigidity of SARS-CoV

Spike Receptor-Binding Motif the Hallmark for Its Enhanced Infectivity? Insights from All-Atom Simulations

The SARS-CoV-2 exerts a distinctive strategy for interacting with the ACE2 human receptor

Enhanced receptor binding of SARS-CoV-2 through networks of hydrogen-bonding and hydrophobic interactions

Computational simulations reveal the binding dynamics between human ACE2 and the receptor binding domain of SARS-CoV-2 spike protein

Investigation of the genetic variation in ACE2 on the structural recognition by the novel coronavirus (SARS-CoV-2)

Biomolecular Simulations in the Time of COVID-19, and After

Mass spectrometry analysis of newly emerging coronavirus HCoV-19 spike S protein and human ACE2 reveals camouflaging glycans and unique post-translational modifications

Virus-receptor interactions of glycosylated SARS-CoV-2 spike and human ACE2 receptor

Comprehensive characterization of N-and O-glycosylation of SARS-CoV-2 human receptor angiotensin converting enzyme 2

Mutations in SLC6A19, encoding B0AT1, cause Hartnup disorder

I-TASSER server for protein 3D structure prediction

I-TASSER: A unified platform for automated protein structure and function prediction

Protein-ligand binding site recognition using complementary binding-specific substructure comparison and sequence profile alignment

CHARMM-GUI Glycan Modeler for modeling and simulation of carbohydrates and J o u r n a l P r e -p r o o f glycoconjugates

Glycan reader: Automated sugar identification and simulation preparation for carbohydrates and glycoproteins

CHARMM-GUI: A web-based graphical user interface for CHARMM

Membrane lipids: Where they are and how they behave

Membrane lipid composition: Effect on membrane and organelle structure, function and compartmentalization and therapeutic avenues

CHARMM36 all-atom additive protein force field: Validation based on comparison to NMR data

Comparison of simple potential functions for simulating liquid water

Scalable molecular dynamics with NAMD

Particle-mesh Ewald: An N.log(N) method for Ewald sums in large systems

Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of nalkanes

MDTraj: a modern open library for the analysis of molecular dynamics trajectories

VMD-Visual Molecular Dynamics

PTRAJ and CPPTRAJ: Software for processing and analysis of molecular dynamics trajectory data

2020. Dual nature of human ACE2 glycosylation in binding to SARS-CoV-2 spike

Beyong Shielding: The Roles of Glycans in SARS-CoV-2 Spike Protein

Structure of the SARS-CoV-2 spike receptor-binding domain bound to the ACE2 receptor

ff 14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff 99SB

GLYCAM06: A generalizable biomolecular force field. Carbohydrates

Accuracy limit of rigid 3-point water models

Systematic parameterization of monovalent ions employing the nonbonded model

Long-time-step molecular dynamics through hydrogen mass repartitioning

A Community Letter Regarding Sharing Biomolecular Simulation Data for COVID-19

ACE2 gene variants may underlie interindividual variability and susceptibility to COVID-19 in the Italian population

ACE2 polymorphisms and individual susceptibility to SARS-CoV-2 infection: insights from an study

The sequence of human ACE2 is suboptimal for binding the S spike protein of SARS coronavirus 2

Human ACE2 receptor polymorphisms predict SARS-CoV-2 susceptibility