The SARS-CoV-2 Exerts a Distinctive Strategy for Interacting with the ACE2 Human Receptor

The COVID-19 disease has plagued over 200 countries with over three million cases and has resulted in over 200,000 deaths within 3 months. To gain insight into the high infection rate of the SARS-CoV-2 virus, we compare the interaction between the human ACE2 receptor and the SARS-CoV-2 spike protein with that of other pathogenic coronaviruses using molecular dynamics simulations. SARS-CoV, SARS-CoV-2, and HCoV-NL63 recognize ACE2 as the natural receptor but present a distinct binding interface to ACE2 and a different network of residue–residue contacts. SARS-CoV and SARS-CoV-2 have comparable binding affinities achieved by balancing energetics and dynamics. The SARS-CoV-2–ACE2 complex contains a higher number of contacts, a larger interface area, and decreased interface residue fluctuations relative to the SARS-CoV–ACE2 complex. These findings expose an exceptional evolutionary exploration exerted by coronaviruses toward host recognition. We postulate that the versatility of cell receptor binding strategies has immediate implications for therapeutic strategies.


Introduction
The coronavirus disease 2019 (COVID- 19), initially detected in the Wuhan seafood market in the Hubei province of China [1], is caused by the SARS-CoV-2 (referred to as the COVID-19 virus for clarity). The COVID-19 virus spread within three months from its appearance to more than 200 countries, resulting in over 200,000 deaths worldwide. The COVID-19 virus is capable of human-to-human transmission and was introduced to humans in a zoonotic event [2].
Currently, seven confirmed coronavirus (CoV) species are known as human pathogens [3]. Four CoVs are endemic species in humans and cause mild respiratory symptoms, mostly in pediatric patients [4]. These are the HCoV-HKU1 and the HCoV-OC43 from the Betacoronavirus (BCoV) genus, and the HCoV-229E and the HCoV-NL63 from the Alphacoronavirus (ACoV) genus. The other human CoVs have caused severe outbreaks. The SARS-CoV (referred to as the SARS-2002 virus for clarity), is a BCoV that emerged in humans in 2002, giving rise to the severe acute respiratory syndrome (SARS) outbreak. The Middle East Respiratory Syndrome (MERS) BCoV caused an outbreak in 2012-2013. Most recently, the SARS-CoV-2, with high homology to the 2002 SARS-CoV, caused the current COVID-19 pandemic [5].
To gain access to host cells, coronaviruses rely on spike proteins, which are membrane-anchored trimers containing a receptor-binding S1 segment and a membrane-fusion S2 segment [6].

Materials and Methods
The structural model of the COVID-19 spike protein receptor-binding domain (RBD) in complex with ACE2 was generated by comparative modeling using MODELLER 9.18 [16] with the COVID-19 sequence (RefSeq: YP_009724390.1). We relied on the crystal structure of the spike protein receptor-binding domain from a SARS coronavirus-designed human strain complexed with the human receptor ACE2 (PDB 3SCI, resolution 2.9 Å) as a template for comparative modeling. The SARS-2002 spike protein RBD and the HCoV-NL63 in complex with ACE2 were taken from PDB 2AJF (resolution 2.9 Å) and 3KBH (resolution 3.3 Å), respectively. The missing residues were added in MODELLER. The MERS RBD structure was taken from the complex with the neutralizing antibody CDC2-C2 (PDB 6C6Z, resolution 2.1 Å) and structurally aligned onto the SARS-2002 RBD in complex with the ACE2 receptor. The designed SARS-CoV variant is from PDB 3SCI.
The MD simulations were performed with GROMACS 2020 software [17] using the CHARMM36m force field [18]. Each of the complexes was solvated in transferable intermolecular potential with 3 points (TIP3P) water molecules, and ions were added to equalize the total system charge. The steepest descent algorithm was used for initial energy minimization until the system converged at Fmax < 1000 kJ/(mol · nm). Then, water and ions were allowed to equilibrate around the protein in a two-step equilibration process. The first part of the equilibration was at a constant number of particles, volume, and temperature (NVT). The second part of the equilibration was at a constant number of particles, pressure, and temperature (NPT). For both of the MD equilibration parts, positional restraints of k = 1000 kJ/(mol · nm 2 ) were applied to the heavy atoms of the protein, and the system was allowed to equilibrate at a reference temperature of 300 K, or reference pressure of 1 bar for 100 ps at a time step of 2 fs. Following the equilibration, the production simulation duration was 100 nanoseconds with 2 fs time intervals. Altogether 10,000 frames were saved for the analysis at intervals of 10 ps. The MD trajectories are available at ftp://ftp.cs.huji.ac.il/users/michall/covid19_trajectories/.
The interaction scores between the virus spike RBD and the ACE2 were calculated for each frame of the trajectory using the statistically optimized atomic potentials (SOAP) [22]. In the interface contact analysis, a residue-residue contact was defined based on the inter-atomic distance, with a cutoff of 4CH. The interaction buried surface area was calculated by subtracting the surface area of the protein complex from the sum of the surface areas of the individual proteins. The surface area was calculated as the solvent accessible surface [23].

Comparable Interaction Energies for ACE2-Binding BCoVs
The COVID-19 RBD (residues 319-529) shares a 72.8% sequence identity and high structural similarity with the SARS-2002 RBD (Table 1). In contrast, the RBD of HCoV-NL63 is only 17.1% identical to that of COVID-19, and there are no significant structural similarities between them ( Figure S2). Remarkably, the RBD of MERS-CoV, which is structurally similar to that of COVID-19 (20.1% sequence identity, 65% structure similarity), recognizes a different host receptor (DPP4) for its cell entry and does not bind to ACE2 [13]. We ran 100 ns molecular dynamic (MD) simulations of ACE2 in complex with the RBDs of the COVID-19, SARS-2002, and HCoV-NL63 viruses to quantify the energetics and the dynamics of the different RBD-ACE2 interactions. The simulation trajectory snapshots at 10 ps intervals (10,000 frames) were analyzed using a statistical potential to assess the probability of the RBD-ACE2 interaction (SOAP score [22]), with lower values corresponding to higher probabilities, and thus higher affinities. The interaction scores (SOAP) for COVID-19 RBD-ACE2 were comparable to those of SARS-2002, medians of −1865.9 and −1929.5, respectively ( Figure 1A,B). HCoV-NL63 has RBD-ACE2 interaction scores that are higher than both of the SARS-CoVs (median of −941.6). MERS, which is structurally similar to COVID-19 (Table 1) does not bind to ACE2. The MERS virus which binds to dipeptidyl peptidase-4 (DPP4, also known as CD26 [13]), has RBD-ACE2 interaction scores that indicate an extremely weak affinity (median of −692.6), as expected from a non-cognate receptor interaction. COVID-19 has the largest buried surface area at the interface (1204 Å 2 ), followed by the interface area for SARS-2002 (998 Å 2 ) and HCoV-NL63 (973 Å 2 ). The number of ACE2 contacting residues maintains the same order, with 30, 24, and 23 for COVID-19, SARS-2002, and HCoV-NL63, respectively ( Figure 1C). The three RBDs exploit specific binding sites on ACE2 based on the analysis of the MD trajectories ( Figure 1C,D, Movie S1). There is a significant overlap of ACE2 interacting residues between COVID-19 and SARS-2002 (at least 73%), while HCoV-NL63 shares only 17% and 36% of contacts with SARS-2002 and COVID-19, respectively. These findings suggest that the coronaviruses exert different interaction strategies with their cognate receptors to achieve the affinity that is required for effective cell entry.

Accelerated Evolution among the Key Anchoring Residues of the RBD-ACE2 Interface
While the sequence identity between the RBDs of COVID-19 and SARS-2002 is 73% (Table 1, Figure S2), we observe a significantly higher residue substitution rate at the interaction interface with the ACE2 receptor. Out of 29 RBD interface residues, only 10 residues (34%) in COVID-19 are conserved with respect to SARS-2002 ( Figure 2A, Table S1, Figure S2). Similarly, only 12 residues (40%) in SARS-2002 are conserved with respect to COVID-19. and HCoV-NL63. An ACE2 residue is considered as part of the interface if one of its atoms is within 4Å from any RBD atom in at least 10% of the 10,000 MD simulation frames; (D) overlay of 50 snapshots for each of the three RBDs. The ACE2 is in surface representation (gray). The frames were aligned using the N-terminal fragment of the ACE2 that contains the two helices participating in the RBDs binding.

Accelerated evolution among the key anchoring residues of the RBD-ACE2 interface
While the sequence identity between the RBDs of COVID-19 and SARS-2002 is 73% (Table 1, Figure S2), we observe a significantly higher residue substitution rate at the interaction interface with the ACE2 receptor. Out of 29 RBD interface residues, only 10 residues (34%) in COVID-19 are conserved with respect to SARS-2002 ( Figure 2A, Table S1, Figure S2). Similarly, only 12 residues (40%) in SARS-2002 are conserved with respect to COVID-19.
To investigate these interface residues, we constructed and overlaid the contact maps for the RBD-ACE2 interfaces for COVID-19 and SARS-2002 ( Figure 2B). We define a residue-residue contact frequency (CF) as the fraction of MD trajectory frames in which the contact appears. Remarkably, only eight out of the total 72 residue-residue interface contacts have comparable (< 50% difference) contact frequencies between the COVID-19-ACE2 and SARS-2002-ACE2 interfaces ( Figure 2B, colored gray). Furthermore, we find two interaction patches unique to COVID-19 ( Figure 2B, patches 1 and 3) and another patch unique to SARS-2002 ( Figure 2B, patch 2). COVID-19 has a significant and unique contact site between the residues 500-505 of the RBD and residues 353-357 of ACE2 ( Figure  2B and C). COVID-19 also creates a new interaction patch with the middle of the N-terminal ACE2 helix ( Figure 2B and C), while SARS-2002 has a unique interaction patch with the end of the same helix ( Figure 2B and C). The rest of the changes in the interface contact frequencies are due to the and HCoV-NL63. An ACE2 residue is considered as part of the interface if one of its atoms is within 4 Å from any RBD atom in at least 10% of the 10,000 MD simulation frames; (D) overlay of 50 snapshots for each of the three RBDs. The ACE2 is in surface representation (gray). The frames were aligned using the N-terminal fragment of the ACE2 that contains the two helices participating in the RBDs binding.
To investigate these interface residues, we constructed and overlaid the contact maps for the RBD-ACE2 interfaces for COVID-19 and SARS-2002 ( Figure 2B). We define a residue-residue contact frequency (CF) as the fraction of MD trajectory frames in which the contact appears. Remarkably, only eight out of the total 72 residue-residue interface contacts have comparable (<50% difference) contact frequencies between the COVID-19-ACE2 and SARS-2002-ACE2 interfaces ( Figure 2B, colored gray). Furthermore, we find two interaction patches unique to COVID-19 ( Figure 2B, patches 1 and 3) and another patch unique to SARS-2002 ( Figure 2B, patch 2). COVID-19 has a significant and unique contact site between the residues 500-505 of the RBD and residues 353-357 of ACE2 ( Figure 2B,C). COVID-19 also creates a new interaction patch with the middle of the N-terminal ACE2 helix ( Figure 2B,C), while SARS-2002 has a unique interaction patch with the end of the same helix ( Figure 2B,C). The rest of the changes in the interface contact frequencies are due to the different interface loop conformations (COVID-19 residue numbers 474-498, SARS-2002 residue numbers 461-484) (Figure 2A,B, Table S1). COVID-19 has a significantly higher number of well-defined contact pairs compared to SARS-2002: 52 vs. 28 contacts (with 44 and 20 unique pairs, excluding the ones with similar CFs) were found for RBD-ACE2 of the COVID-19 and SARS-2002, respectively ( Figure 2B). These results expose the accelerated evolution among the key anchoring residues of the RBD-ACE2 interface and raise the following question: how does SARS-2002 RBD reach an ACE2 binding affinity that is comparable to that of COVID-19 but with fewer contact pairs and a smaller interface area?
Viruses 2020, 12, 497 5 of 10 numbers 461-484) (Figure 2A and B, Table S1). COVID-19 has a significantly higher number of welldefined contact pairs compared to SARS-2002: 52 vs. 28 contacts (with 44 and 20 unique pairs, excluding the ones with similar CFs) were found for RBD-ACE2 of the COVID-19 and SARS-2002, respectively ( Figure 2B). These results expose the accelerated evolution among the key anchoring residues of the RBD-ACE2 interface and raise the following question: how does SARS-2002 RBD reach an ACE2 binding affinity that is comparable to that of COVID-19 but with fewer contact pairs and a smaller interface area?

COVID-19 and SARS-2002 Differ in Their Fluctuation Signatures
The distribution of SOAP scores throughout the simulation trajectory has a larger fluctuation range for SARS-2002, relative to COVID-19 ( Figure 1A,B, Figure S3A) suggesting that the SARS-2002-ACE2 interaction is fluctuating between several structural states. Moreover, the analysis of contact frequencies along the entire trajectory reveals that none of the SARS-2002 contacts are maintained over 90% of the frames, while COVID-19 still maintains about half of its contacts for 90% of the trajectory ( Figure S3B).
To investigate the dynamics of COVID-19 binding compared to SARS-2002, we calculated the root-mean-square fluctuation (RMSF) of each residue with respect to the lowest energy snapshot from their respective 100 ns MD simulation trajectory. The interface region in the RBD contains two loops (loop1: residues 474-489, loop2: residues 498-505; using COVID-19 numbering, Figure 3A,D) that bind to the ACE2 N-terminal helix on both of its ends. These two loops are highly flexible in the SARS-2002 RBD ( Figure 3A,D). While loop1 is also fluctuating in the COVID-19 RBD, albeit much less, loop2 remains relatively rigid in the COVID-19 RBD. In addition, we find that in the COVID-19-RBD, a region centered around K417 leads to further stability relative to the corresponding region in SARS-2002. We attribute this difference to the unique interaction of COVID-19 at position K417 with the middle of the N-terminal ACE2 helix, thus serving as an anchor site to the receptor ( Figures 2C and 3A). The contribution of K417 to ACE2 binding is observed in a recent cryo-EM structure of the COVID-19 spike protein bound to ACE2 [24]. Overall, COVID-19 is more rigid compared to SARS-2002 ( Figure 3A,D).
We investigated the dynamics of a designed SARS (SARS-des) variant [11], which differs from SARS-2002 at only two positions: Y455F and L486F (the corresponding positions in SARS-2002 are 442 and 472, respectively; Table S1). The L486F mutation is of special interest for the COVID-19 RBD because it has this same substitution. Our MD simulation analysis reveals that the SARS-des has a substantially lower interaction score with ACE2 (median of −2199.2, Figure S3), as expected for an optimized human ACE2-binding RBD design. We observed that these two mutations not only enhance the binding affinity to ACE2, but also lead to a substantial stabilization of the interaction interface. The fluctuation signatures along the RBD of SARS-des are surprisingly similar to those recorded for COVID-19 ( Figure 3B,C). Thus, the switch from a flexible binding mode (for SARS-2002) to a stable one (COVID-19 and SARS-des, Figure 3B) highlights the remarkable capacity of the RBD to adopt alternative receptor binding strategies driven by a minimal number of amino acid substitutions. This analysis reveals the critical role of L486F (SARS-des residue F472) for stabilizing the COVID-19-ACE2 interface and producing a reduction in the number of states of the COVID-19 spike protein bound to an ACE2 receptor.

Discussion
The experimental affinity measurements (e.g., surface plasmon resonance (SPR)) confirm the high affinity of SARS-2002 RBD-ACE2 binding, with an equilibrium dissociation constant (K D ) of~1.5-10.0 nM [25][26][27], comparable to the binding affinity of ACE2 and the COVID-19 RBD (K D =~1.2-14.7 nM [27,28]) and consistent with our MD-based calculation ( Figure 1A, Figure S3, and Table S2). While we relied on the modeled COVID-19 RBD-ACE2 interaction, the median interaction score of the complex (−1928.8), based on the MD simulation starting from the recently published X-ray structure (PDB 6LZG [21]), was even closer to the SARS-2002 RBD-ACE2 score (−1929.5). Binding affinity is achieved through a combination of interface contact optimization and protein stability ( Figure 3E). While the RBD-ACE2 complex can be resolved at high-resolution by cryo-EM [24] and X-ray crystallography, MD simulations provide orthogonal information about the interaction dynamics on a nanosecond timescale. In the case of CoVs, MD simulations suggest an exceptional versatility of viral receptor binding strategies ( Figure 3E). COVID-19 adopted a different strategy for achieving comparable affinity to SARS-2002: the interface of COVID-19 is significantly larger than that of SARS-2002 (1204 Å vs. 998 Å) with a remarkable number of interacting residues (ACE2: 30 vs. 24, Figure 1C). In contrast, SARS-2002 is more flexible in its interaction with ACE2, Viruses 2020, 12, 497 7 of 10 interacting through fewer contacts that serve as "hot spots". Therefore, we predict that SARS-2002 RBD neutralizing antibodies will not be effective for COVID-19. The failure of several of these antibodies to neutralize the binding of COVID-19 RBD to its receptor is consistent with our findings [28,29]. The fluctuation from high-to low-affinity conformations in SARS-2002 led to an increased efficacy for inhibiting peptides [30] and high-affinity antibodies [31] compared to COVID-19. This implies that a therapeutic challenge is attributed to the enhanced rigidity of the COVID-19 RBD relative to that of the SARS-2002.

Discussion
The experimental affinity measurements (e.g., surface plasmon resonance (SPR)) confirm the high affinity of SARS-2002 RBD-ACE2 binding, with an equilibrium dissociation constant (KD) of ~1.5-10.0 nM [25][26][27], comparable to the binding affinity of ACE2 and the COVID-19 RBD (KD = ~1.2-14.7nM [27,28]) and consistent with our MD-based calculation ( Figure 1A, Figure S3, and Table S2). While we relied on the modeled COVID-19 RBD-ACE2 interaction, the median interaction score of the complex (-1928.8), based on the MD simulation starting from the recently published X-ray structure (PDB 6LZG [21]), was even closer to the SARS-2002 RBD-ACE2 score (-1929.5). Binding residues (Tyr, Trp, Phe) in the interface with ACE2 ( Figure 2A). Moreover, in the SARS-designed variant (SARS-des [11]), the addition of an aromatic residue (L486F substitution) significantly improved the interaction scores and interface stability ( Figure 3B,D). Our findings shed light on the accelerated evolution of spike protein binding to the ACE2 receptor, similar to the rapid evolution along the antibody-antigen affinity maturation process.