Evaluation of Interactions between SARS-CoV-2 RBD and Full-Length ACE2 with Coarse-Grained Molecular Dynamics Simulations.

Compared to all-atom models, coarse-grained models enable the investigation of the dynamics of simulation systems on a much larger length scale and a longer time scale, which makes them suitable for studying macromolecular systems. Hence, in this work, we performed multiple μs-scale Martini coarse-grained molecular dynamics simulations to reveal the interaction details between SARS-CoV-2 RBD and full-length human ACE2. Besides, the key coarse-grained systems were backmapped into the corresponding all-atom system for the display of structural details. Our results indicated that the plier structure in two ends of the binding interface plays a key role in the binding process of SARS-CoV-2 RBD with ACE2. Furthermore, we also found that when there is no B0AT1 in the simulation system, the N-terminus of ACE2 is more likely to approach the cell membrane, which has a strong correlation with the subsequent fusion of the virus with the cell membrane. These binding details of SARS-CoV-2 RBD and the ACE2 protease domain (PD) as well as the membrane orientation thermodynamics can promote the development of therapeutic drugs and preventive vaccines against SARS-CoV-2.

Introduction

The outbreak of the disease 2019 (COVID-19) caused by severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2) has become a worldwide public health crisis. Understanding why SARS-CoV-2 spreads so quickly and widely has a high priority for unveiling its mysterious mask and curbing its spread. The strong infectivity of SARS-CoV-2 is mainly from the higher binding affinity of the receptor binding domain (RBD) of the SARS-CoV-2 spike protein to human angiotensin-converting enzyme 2 (ACE2).[1,2] The binding process between SARS-CoV-2 and ACE2 is the first and key step for the cell entry of SARS-CoV-2.[1] In view of this, the RBD monomer,[3] tandem RBD dimer,[4] and S-trimer[5] are being developed into vaccines. Hence, understanding the interaction details between SARS-CoV-2 RBD and ACE2 will promote the development of efficient therapeutics. Full-length ACE2 consists of an N-terminal protease domain (PD) and a C-terminal collectrin-like domain (CLD), which contains a single transmembrane (TM) helix and a ferredoxin-like fold domain (neck domain) linking TM and PD.[6] In the recognition of RBD, ACE2 PD mainly engages the α1-helix with a minor contribution from the α2-helix and the linker of the β3- and β4-sheets.[6,7] To engage a host-cell receptor and evade the immune surveillance, RBD undergoes transient hinge-like conformational motions (opened or closed conformation).[7] Among the opened conformations, the ACE2 binding strength becomes stronger as the RBD angle increase.[8] SARS-CoV-2 is thought to occupy the active, open RBD state less often than in SARS-CoV to evade the immune surveillance of the human body, but it is compensated by the greater affinity of the SARS-CoV-2 RBD for ACE2.[1,9] Hence, the identification of the structural determinants of these interactions is a crucial step toward a deeper understanding and regulation of the binding process.

To date, methods, which are used for detecting the interacting pattern between ACE2 and SARS-CoV-2, are two main categories: (1) structural biology and (2) molecular dynamics (MD) simulations. Structural biology with a high resolution of the angstrom scale can get crystal structures of the SARS-CoV-2 S-ACE2 complex,[6,10,11] but this static structure cannot reveal the dynamic information of the binding complex of SARS-CoV-2 RBD and ACE2. However, the dynamic state of the RBD is essential to explain the paradox in the current research about SARS-CoV-2.[1] Molecular dynamics simulations can properly reflect dynamic features of the system by recording the trajectory of every atom. Then, theoretical insights of binding details between SARS-CoV-2 RBD and ACE2 can be provided based on the obtained trajectory. Using the all-atom model, main binding details,[12,13] binding free energy,[14] and peptide inhibitors[15] have been continuously revealed and found. However, most of the relevant dynamics and interactions within cells (typically, protein–protein docking, rearrangement upon ligand binding or after biochemical reactions, folding) occur on the time scale of microseconds or milliseconds.[16] Both the time scale and the length scale make it too expensive to study the SARS-CoV-2 RBD-ACE2 binding dynamics using the all-atom model. Hence, the study of SARS-CoV-2 RBD-ACE2 based on coarse-grained (CG) model will provide useful complementary insights to current all-atom MD simulations.

Especially, the reduction of the model complexity of macromolecular systems in the CG model is beneficial for retaining the overall features of the simulated system. Hence, in this work, we transferred the all-atom model (PDB ID: 6M17, Figure A) to the CG model (Figure B) to study the binding features on the interface of RBD and ACE2 PD in multiple μs-scale CGMD simulations (Figure C). We obtained the heat map of the average residue distances of SARS-CoV-2 RBD and ACE2 PD using the whole CGMD simulation. The distance, between residue I21-T92 (α1-helix and α2-helix), residue H345-M360 (the linker of the β3- and β4-sheets) on ACE2, and residue F400-P507 (including receptor binding motif, RBM) on SARS-CoV-2 RBD, is revealed by a two-dimensional (2D) heat map. These distances represent the contribution of the perturbed residues to the binding affinity of ACE2 PD and SARS-CoV-2 RBD. Subsequently, we analyzed the effects of RBD binding and the presence of B0AT1 on membrane orientation dynamics of ACE2’s extracellular part using the two-dimensional (2D) Gibbs free energy map.

Figure 1

Established process of the CG simulation system. (A) Initial all-atom structure (PDB ID: 6M17). (B) The corresponding CG structure with the Martini model. (C) Simulation system is composed of the complex of ACE2 with SARS-CoV-2 RBD (B0AT1: cyan, ACE2 PD: green, ACE2 neck: blue, ACE2 TM: magenta, SARS-CoV-2 RBD: red), POPC membrane, and Na+ ions (scattered purple beads).

Results and Discussion

Different Regions of ACE2 PD and SARS-CoV-2 RBD Have Different Degrees of Contribution to Binding Affinity

As shown in Figure , ACE2 PD and RBD have stable binding during the whole simulation time, which is consistent with the biological characteristics of the binding domain. Also, the tertiary and quaternary structures of the complex were properly maintained during CGMD simulations from time evolutions of root-mean-square deviations (RMSDs, Figure S1). To capture the detailed interaction dynamics between different residues in these two proteins, we calculated the 2D average distance map. As shown in Figure , the distance of residues between ACE2 PD and SARS-CoV-2 RBD keeps changing to a mild extent (see in the Supporting Information), but the key residue pairs of the binding interface barely change. Residues in the α1-helix (I21-T55) of ACE2 PD contribute the majority of contact residues to recognize RBD, and those from the α2-helix (E56-T92) and the linker (H345-M360) of the β3- and β4-sheets have a minor contribution to the recognition of RBD. This is consistent with the result of the cryo-EM structure[6] and insights of potential interacted sites.[7] Han et al.[15] also proved that the peptide inhibitors, formed by the α1-helix and the α2-helix extracted from ACE2 PD, provide a highly specific and stable binding to SARS-CoV-2 RBD. From the amino acid sequence of SARS-CoV-2 RBD, the contact area with ACE2 is divided into five regions: (1) I402-E406, (2) G416-N422, (3) S443-S459, (4) I472-Q474, and (5) E484-Y505. Clearly, regions 3 and 5 contribute more contact residues than other regions of RBD. From the spatial structure of the interface of ACE2 PD and SARS-CoV-2 RBD, the contact residues are mainly distributed in the loop regions, two of which are located at both ends of the dimer contact region, denoted CR1 and CR3; the middle region in the protein–protein interface (CR2) consists of two short strands of β-sheets bridging across the N-terminal helix of ACE2.[13]

Figure 2

Time evolution of simulation systems with (A) two, (B) one, and (C) zero RBD molecules. The coloring style is the same as Figure .

Figure 3

Two-dimensional (2D) distance heat map between residues of ACE2 PD and SARS-CoV-2 RBD. The warmer the color, the closer the distance between two amino acids, and vice versa.

Critical Residues in the Binding Interface between ACE2 PD and SARS-CoV-2 RBD

ACE2 is a homodimer protein, in which the PD is located at the N-terminus, and is the main binding domain to SARS-CoV-2 RBD. In each of our 10 simulation systems (Figure C), we have the dimer structure of the SARS-CoV-2 RBD-ACE2 complex. In other words, we have 2 × 10 independent contact data for our statistics. Based on this data, we counted the frequency of the residues on the interface between ACE2 PD and SARS-CoV-2 RBD. We picked up the residues of ACE2 presenting more than 15 times (maximum: 20), which include E23, Q24, T27, F28, D30, K31, H34, E37, Y41 (located at the α1-helix of ACE2 PD), Y83 (located at the α2-helix of ACE2 PD), G353, G354, and G356 (located at the linker of the β3- and β4-sheets of ACE2 PD; Figure A). The high frequency of Q24, T27, F28, K31, Y41, Y83, K353, E23, and H34 are consistent with their binding free energy.[2,14] However, M82, which forms a hydrophobic pocket[13] together with F28 (α1-helix), L79, and Y83 (α2-helix), only appeared once (Table S1), which may be related to its poor stability. Because M82 is not a core residue of the hydrophobic pocket (e.g., Y83), and not at the stable secondary structure that is preserved (e.g., F28 and L79) along the MD simulations. The interfacial interactions in the middle region across the N-terminal helix of ACE2 are dominated by hydrophobic interactions involving residues D30, K31, H34, and E37, which all have a high frequency. These hydrophobic contacts have an essential role in enhancing the binding affinity in the interacting state.

Figure 4

(A) Frequency of residues (>15) in ACE2 PD in contact with SARS-CoV-2 RBD. (B) The frequency of residues (>15) in SARS-CoV-2 RBD in contact with ACE2 PD. (C) Structural details of ACE2 (the number in parentheses is the contact frequency). (D) Structural details of SARS-CoV-2 RBD.

The residues, occurring more than 15 times (maximum: 20), on SARS-CoV-2 RBD are shown in Figure B. All of the residues in Figure B possess a high level of binding free energy[2,14,17] except for three residues (Y495, E484, Y473). Y495 is an essential residue of SARS-CoV-2 RBD to contact with the hydrophobic area formed by D30, K31, H34, and E37 on the middle region across the N-terminal helix of ACE2. However, it does not contact a specific residue in ACE2 PD but multiple residues (Figure A). E484 and Y473 mainly connect to residues of the terminal of the α1-helix because the terminal residues are more flexible than the nonterminal residues. In general, Y449, L455, N487, Y489, Q493, Y495, and Y505 interact with the hydrophobic area formed by D30, K31, H34, and E37 in the middle region of ACE2. It was proposed that the polar interactions between CR2 of the RBD and the middle region of the N-terminal helix on ACE2 provide large contributions to the binding of the SARS-CoV-2 RBD.[13,18,19] F486 is anchored to the hydrophobic pocket, formed by F28, L79, and Y83. Q498 interacts with the linker of the β3- and β4-sheets of ACE2, while F456, Y473, and E484 mainly interact with the terminal of the α1-helix.

Figure 5

(A) Frequency of stable residue pairs (>10) in the binding interface between ACE2 PD and SARS-CoV-2 RBD. (B, C) Structural details of two pliers. (D) Structural details of the middle region of the binding interface. (E) Two pliers (the ellipses with blue and red colors) are located at the different extremes of the binding interface.

Critical Residue Pairs in the Binding Interface between ACE2 PD and SARS-CoV-2 RBD

During our μs-scale MD simulations, the residue pair in the binding interface is dynamic. Based on our statistics over the 2 × 10 RBD-ACE2 complex, the residue pairs that appeared more than 10 times were considered stable and are displayed in Figure A. F486 forms one end of the plier and is anchored to the hydrophobic pocket composed of F28, L79, and Y83 in the extreme of ACE2 PD (Figure B). Another end of the plier is composed of Y505 and Q498 that are responsible for binding to the linker of the β3- and β4-sheets of ACE2 PD and residues of the middle hydrophobic area of the α1-helix (Figure C). Y505 was defined as the viral determinants for the specific recognition of SARS-CoV-2 by ACE2.[10] A previous study reported that Q498 formed an interaction with the K353 of ACE2, and this interaction is sustained during MD simulations.[20] The plier feature of RBD is consistent with the binding process and binding characteristics that in the first step of the encounter complex of ACE2 and SARS-CoV-2 RBD, only CR3 (including Y495-Y505) contacts with ACE2, while in the second step, both CR1 (including E471-Y489) and CR3 contact with ACE2.[12] The two ends of the plier determine the key position of binding with ACE2, and then the other residues located on the RBD further bind with ACE2 to transfer the encounter state to the interacting state, enhancing the stability of the binding interface (the position of the plier is shown in Figure E).

Y489/Y495/L455/Q493 connect with the middle hydrophobic region of the α1-helix on ACE2 PD (Figure D), and these connections are the most important part of CR2. CR2 was involved in the interacting state, which is the next step of the encounter state mediated by the plier’s structure of RBD. Among these residues, Y495 interacts closely with multiple residues of the hydrophobic region of the α1-helix on ACE2 (Figure A), and its frequency is top three together with Y505 and F486, compared to other residues in RBD (Figure B). Hence, it should be the key residue to connect to the middle hydrophobic region across the N-terminal helix of ACE2. Besides, Y473, N487, E484, and F456 mainly contact with the terminal of the α1-helix. The connection of this part is mainly due to the higher flexibility of residues on the terminal of the α1-helix. K417 (19 times in Figure B) does not establish a stable connection, which occurs more than 10 times, with any specific residue in ACE2, but it connects to multiple residues in ACE2 PD, such as K26 (9 times), T92 (8 times), and D30 (7 times, Table S2). This result means that K417 is a flexible residue in connecting to residues in ACE2 PD and thus not the suitable target for the design of inhibitors. Compared to connections of residues in two pliers, connections between other residues are more flexible and changeable.

From the perspective of the spatial structure, the binding interface of the ACE2 PD side is a relatively rigid structure containing multiple secondary structures, especially the α1-helix and α2-helix of the N-terminal of ACE2, while the SARS-CoV-2 RBD is composed of several β-sheets followed by a flexible region RBM, a domain similar to the sucker tentacles of a gecko. The flexible region of RBD leads to a higher possibility of binding to ACE2 PD and provides a greater binding affinity. These features make the combination between SARS-CoV-2 RBD and ACE2 PD very stable. Because of the flexibility of RBM, the key amino acids on the contact interface obtained by researchers always have some differences, which is not conducive to the design of antiviral drugs based on the structural characteristics of RBM. However, the plier structure of RBD seems to be the first touchpoint in the process of binding to ACE2 PD. Hence, if we can design drugs and vaccines against the structure of the plier, it is possible to block the virus recognition process. The current mutations of SARS-CoV-2 RBD are concentrated in the non-plier region,[21] which means that the therapeutics targeted at the plier structure could bypass the mutation site and relieve the pressure of the epidemic to a large extent.

Absence of B0AT1 in the Simulation System Would Significantly Increase the Possibility that ACE2 PD Approaches the Cell Membrane

When there is B0AT1 in the simulation system, with zero, one, or two RBDs binding to ACE2, the distance from ACE2 PD to the 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) membrane does not change significantly, mainly at 9.0–10.5 nm, and the ratios of the main region are 80.65, 77.94, and 84.3%, respectively (Figures S3 and 6A–C). However, when B0AT1 was removed from the simulation system, the membrane orientation dynamics of ACE2 was dramatically changed (Figure ). In the absence of B0AT1, the distance between ACE2 PD and the POPC membrane can be much closer, and there is another local minimum free energy region on the left side in Figure D–F, which indicates that in this condition, ACE2 PD exhibits a significant trend of approaching the membrane, and ACE2 neck has no such trend. In addition, it can also be found that as the number of RBDs in the system increases from 0 to 2, a transition zone between the left and right minimum free energy regions gradually appears. According to the distribution of points in the 2D free energy heat map, we selected points with Gibbs free energy less than 4.48 kJ mol–1 (the point distribution less than this free energy value is more concentrated and characteristic) for ratio statistics (Figure S2). The statistical results show that the proportion of transition zone points is positively correlated with the number of RBDs in the simulation system without B0AT1. However, in the case of 0 RBD, the ratio of the transition zone plus the minimum free energy region on the left is mildly greater than that in the case of 1 RBD or 2 RBD, which is consistent with the phenomenon that the presence of RBD in the system with B0AT1 does not increase the probability of ACE2 PD approaching the membrane. To ensure that the contributions of different simulation groups to the results are similar, we counted the percentages of points in different regions within the group (Figures S4–S6), which indicates that there is no extreme contribution in each simulation to Figure D–F. The significance test (Figure S7) shows that there is no significant difference in the right local minimum free energy region of Figure D–F, but there are significant differences in the left local minimum free energy region and the middle transition region of Figure D–F. In general, when there is no B0AT1 in the system, ACE2 has a greater tendency to approach the membrane, and with the increase of the number of binding RBDs in the system, the proportion of the transition region between the left and right minimum free energy regions gradually increases.

Figure 6

Two-dimensional Gibbs free energy map using ACE2-PD-membrane and ACE2-neck-membrane distances as two reaction coordinates. Cases of zero (A, D), one (B, E), or two (C, F) RBD binding were considered. (A–C) Systems with B0AT1 and (D–F) systems without B0AT1.

Conclusions

Through multiple μs-scale coarse-grained MD simulations of ACE2-RBD complexes, we found that the main parts mediating the combination of the flexible SARS-CoV-2 RBD and the rigid ACE2 PD are two ends of the interface. One end is F486 anchored to the hydrophobic pocket formed by F28, L79, and Y83 of ACE2 PD. The other end is made of Q498 and Y505, both of which are mainly responsible for the hydrophobic connection with the linker of the β3- and β4-sheets and the α1-helix in ACE2 PD. Among all of the amino acids that contribute to the binding free energy between ACE2 and SARS-CoV-2 RBD, F486, Q498, and Y505 are relatively ranked in the top part.[8] These two end structures are of great significance for the initial recognition of ACE2 by RBD. If the recognition of the corresponding structure on ACE2 by these two end structures can be curbed, the recognizing capability of the whole RBD to ACE2 may be reduced a lot. Besides, our results indicated that the absence of B0AT1 greatly increases the probability of ACE2 PD approaching the membrane. This means that the presence of B0AT1 is important for ACE2 to maintain its PD away from the cell membrane, which will prevent the binding with SARS-CoV-2 RBD and thus the virus invasion process.

Methods

ACE2-B0AT1-RBD Complex

The atomistic coordinates of the ACE2-B0AT1-RBD complex were taken from PDB ID 6M17.[6] This complex contained a dimer of the ACE2 receptor in complex with the RBD and also the B0AT1 transporter. Our simulation system contained the POPC membrane-embedded dimeric form of the ACE2-B0AT1 with one or two or without RBD molecules in different simulation setups, and we supplemented the control simulations without B0AT1 to illustrate the role of B0AT1.

Martini Force Field

As a popular CG model, the Martini force field (version 2.2)[22−25] was used in the current work. In this model, generally, four heavy atoms are mapped into one interaction site, including four main types: polar (P), nonpolar (N), apolar (C), and charged (Q). Four different subtypes (d = donor, a = acceptor, da = both, 0 = none) were introduced to bead types of N and Q to mimic hydrogen bonding capacities and allow a fine representation of the chemical nature. For bead types of P and C, five different subtypes (from 1, low polarity, to 5, high polarity) were used to describe the degree of polarity.

Molecular Dynamics Simulations

The CGMD simulation of all systems was performed using the GROMACS program v5.1.2[26] and the Martini force field,[22,24,27] while the visualization of system snapshots was done using VMD[28] and Pymol (Version 2.4.0).[29] Unless stated otherwise, systems were first energy-minimized using the steepest descent algorithm for 40000 steps. Then, the system was relaxed after 1 ns simulations using a 20 fs time step and the NPT ensemble. Last, simulations were run for 1.2 μs with a time step of 30 fs. For all simulations, periodic boundary conditions were applied in three dimensions. A v-rescale thermostat[30] with a relaxation time τ = 1.0 ps was used to maintain a constant temperature of 310 K and a constant pressure of 1 bar was kept by Berendsen pressure coupling[31] (coupling constant of 3.0 ps and compressibility was 3 × 10–4 bar–1) in the NPT ensemble. Semi-isotropic pressure coupling was applied to our systems containing a planar lipid bilayer comprising POPC and a dimeric form of ACE2-B0AT1 with one, two, or none of RBD (control simulations without B0AT1 were obtained by removing B0AT1 from these systems). A 1.1 nm cutoff was applied for van der Waals interactions, where the LJ potential was shifted to zero smoothly from 0.9 to 1.1 nm to reduce the cutoff noise. For the Coulombic potential, a 1.1 nm cutoff was used for short-range electrostatic interactions while shifting to zero from 0 to 1.1 nm smoothly. The neighbor list for nonbonded interactions was updated every 10 steps with a cutoff of 1.4 nm. Every simulation was repeated 10 times for the later contact information collections of amino acids.

System Setup

The protein was obtained from the PDB and converted from the all-atom model (Figure A) to the CG model (Figure B) with the automated workflow. This program first converts the atomistic protein to the Martini CG model using script martinize.py (Version 2.6),[32] and the elastic network was applied (spring force constant = 500 kJ mol–1 nm–2, the lower and upper limits of the cutoff distance were 0.5 and 0.9 nm, respectively) to each subunit separately. After which one, two, or none of RBD remained and B0AT1 would be removed or remained according to different setup systems. Subsequently, insane.py[33] was used to generate a bilayer membrane, containing 885 POPC in the upper leaflet and 887 POPC in the lower leaflet, and the solvent in a box of 24 × 24 × 35 nm2 (Figure C). The solvent comprised 14 1960 Martini water molecules and 72 Na+ ions. At the same time as generating the bilayer membrane system, the CG protein was embedded into the bilayer membrane.

2D Contact Map

To determine the specific binding amino acid pairs between the ACE2 PD and RBD, the GROAMCS program gmx mindist was used to calculate the average distance between the amino acids of RBD and the amino acids of ACE2 PD at intervals of 10 frames through the whole trajectory in the CG model. The two-dimensional heat map was generated based on the average distance between amino acids on ACE2 PD and RBD by in-house scripts. Then, to recognize the pivotal amino acids and amino acid pairs, we picked up the frequency of the amino acids and amino acid pairs, whose distance is lower than 0.9 nm, in multiple simulations. Then, to prevent the interactions between adjacent amino acids from affecting the reliability of the results, the amino acid pair with the smallest distance was selected as the vital residue pair in one-to-many amino acid pairs. Moreover, once the residue–residue distance was smaller than 0.5 nm, it was considered a vital residue pair.

2D Gibbs Free Energy Map

A two-dimensional Gibbs free energy map with two reaction coordinates based on the MD trajectories in the CG model was calculated using the equation ΔG = −RT ln(Ω/Ω0), where R is the gas constant (8.31 J·K–1 mol–1) and Ω is the density of states. One reaction coordinate is the center of mass distance along the POPC membrane normal (Z-axis) to ACE2 PD, and another reaction coordinate is the center of mass distance along the POPC membrane normal (Z-axis) to ACE2 neck.

Backmapping

For exhibiting the contact details of ACE2 PD with SARS-CoV-2 RBD with a better view, we performed the backmapping of ACE2 PD and SARS-CoV-2 RBD from the CG model to the all-atom model by backward.py.[34]