key: cord-0829242-i30kia1b
authors: Taka, Elhan; Yilmaz, Sema Z.; Golcuk, Mert; Kilinc, Ceren; Aktas, Umut; Yildiz, Ahmet; Gur, Mert
title: Critical Interactions Between the SARS-CoV-2 Spike Glycoprotein and the Human ACE2 Receptor
date: 2021-05-12
journal: J Phys Chem B
DOI: 10.1021/acs.jpcb.1c02048
sha: 7db335fe786cd8c6dca9c268712786b7d9646b2c
doc_id: 829242
cord_uid: i30kia1b

[Image: see text] Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) infects human cells by binding its spike (S) glycoproteins to angiotensin-converting enzyme 2 (ACE2) receptors and causes the coronavirus disease 2019 (COVID-19). Therapeutic approaches to prevent SARS-CoV-2 infection are mostly focused on blocking S-ACE2 binding, but critical residues that stabilize this interaction are not well understood. By performing all-atom molecular dynamics (MD) simulations, we identified an extended network of salt bridges, hydrophobic and electrostatic interactions, and hydrogen bonds between the receptor-binding domain (RBD) of the S protein and ACE2. Mutagenesis of these residues on the RBD was not sufficient to destabilize binding but reduced the average work to unbind the S protein from ACE2. In particular, the hydrophobic end of RBD serves as the main anchor site and is the last to unbind from ACE2 under force. We propose that blocking the hydrophobic surface of RBD via neutralizing antibodies could prove to be an effective strategy to inhibit S-ACE2 interactions.

Coronavirus disease 2019 (COVID-19) pandemic is caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), which is a positive-sense RNA betacoronavirus. Phylogenetic analyses demonstrated that the SARS-CoV-2 genome shares ∼79% sequence identity with severe acute respiratory syndrome coronavirus (SARS-CoV), and ∼52% with the Middle East respiratory syndrome coronavirus (MERS-CoV). 1 Despite these similarities, SARS-CoV-2 is much more infectious and fatal than SARS-CoV and MERS-CoV together. 2 SARS-CoV-2 consists of a 30 kb single-stranded RNA genome that is encapsulated by a lipid bilayer and three distinct structural proteins that are embedded within the lipid membrane: envelope (E), membrane (M), and spike (S). Host cell entry is primarily mediated by homotrimeric S glycoproteins located on the viral membrane ( Figure 1a ). 3 Each S protomer consists of S1 and S2 subunits that mediate binding to the host cell receptor and fusion of the viral envelope, respectively. 3, 4 The receptor-binding domain (RBD) of S1 undergoes a rigid body motion to bind to angiotensinconverting enzyme 2 (ACE2). In the closed state, all RBDs of the S trimer are in the down position, and the binding surface is inaccessible to ACE2. The switching of one of the RBDs into a semi-open intermediate state is sufficient to expose the ACE2 binding surface and stabilize the RBD in its up position ( Figure  1b) . 5 The S protein binds to the human ACE2 receptor, a homodimeric integral membrane protein expressed in the epithelial cells of the lungs, heart, kidneys, and intestines. 7 Each ACE2 protomer consists of an N-terminal peptidase domain (PD), which interacts with the RBD of the S protein through an extended surface (Figure 1a ,c). 7−9 Upon ACE2 binding, proteolytic cleavage of the S protein by the serine protease transmembrane serine protease 2 (TMPRSS2) separates the S1 and S2 subunits. 10 The S2 protein exposes fusion peptides that insert into the host membrane and promote fusion with the viral membrane. 4 To prevent SARS-CoV-2 infection, there is a global effort to design neutralizing antibodies, 11 nanobodies (single-domain antibodies), 12 peptide or miniprotein inhibitors, 13−15 and small molecules 16, 17 that target the ACE2 binding surface of the S protein. Yet, only a limited number of studies have been performed to investigate critical interactions that facilitate S-ACE2 binding using molecular dynamics (MD) simulations. Initial studies have constructed a homology model of SARS-CoV-2 RBD in complex with ACE2, based on the SARS-CoV crystal structure 9, 18 and performed conventional MD (cMD) simulations to estimate binding free energies 19−21 and interaction scores. 22 More recent studies used the crystal structure of SARS-CoV-2 RBD in complex with ACE2 to perform coarse-grained 23 and all-atom MD simulations, in the presence of an explicit solvent 24−28 and implicit solvent 29 to investigate binding free energy, 23,25−27 binding energy, 29 and unbinding work. 28 The effect of the mutations that disrupt close contact residues between SARS-CoV-2 RBD and ACE2 PD on binding free energy was investigated by post-processing of the MD trajectories 20,30 or using bioinformatic methods. 24 The work required to unbind the SARS-CoV-2 S protein from ACE2 has been estimated via steered MD (SMD) simulations, 28 but these simulations were performed at high pulling velocities in the absence of glycans, and without satisfying the stiff-spring approximation. 31 Structural, 7, 8, 32 biochemical, 11, 33 and computational 20, 27, 30 studies identified the critical residues that stabilize S-ACE2 binding, but the contribution of these interaction pairs to the S-ACE2 binding energy is not well understood. It also remains controversial whether the SARS-CoV-2 S protein binds to ACE2 more strongly than the SARS-CoV S protein. 19, 23, 26, 28, 29, 34, 35 In this study, we performed a comprehensive set of all-atom MD simulations totaling 23.95 μs in length using the recentlysolved structure of the RBD of the SARS-CoV-2 S protein in complex with the PD of ACE2. 8 Simulations were performed in the absence and presence of external force to investigate the binding characteristics and estimate the binding strength.

These simulations showed additional interactions between RBD and PD to those observed in the crystal structure. 8 An extensive set of alanine substitutions and charge reversal mutations of the RBD amino acids involved in ACE2 binding were performed to quantify how mutagenesis of these residues weaken binding in the presence and absence of force in simulations. We showed that the hydrophobic end of RBD primarily stabilizes S-ACE2 binding, and targeting this site could potentially serve as an effective strategy to prevent SARS-CoV-2 infection. ■ METHODS MD Simulations System Preparation. For cMD simulations, the crystal structure of SARS-CoV-2 S protein RBD bound with ACE2 PD at 2.45 Å resolution (PDB ID: 6M0J) 8 was used as a template. The chloride ion, zinc ion, glycans, and water molecules in the crystal structure were kept in their original positions. The protonation states of the titratable residues of both ACE2 and RBD proteins were predicted using the PROPKA web server, 36,37 titratable residues were left in the dominant protonation state at pH 7.0. Single and double point mutants were generated using the Mutator Plugin in Visual Molecular Dynamics (VMD). 38 Each system was solvated in a water box (using the TIP3P water model) having 35 Å cushion in the positive x-direction and 15 The structure of the full-length S protein in complex with ACE2. The S protein is a homotrimer (green, purple, and gray) and embedded into the viral membrane. ACE2 is a homodimer (blue and orange) and embedded into the host cell membrane. The full-length structure of the S protein in complex with ACE2 was modeled using the full-length S protein model 6 and the crystal structure of the S protein RBD in complex with ACE2 (PDB ID: 6M17). Both proteins were manually inserted into the membrane by their transmembrane domains. (b) The structure of an S protomer with its RBD in the down and up positions. S1/S2 and S2′ are the cleavage sites of the S protomer upon ACE2 binding. (c) MD simulations were performed for RBD of the S protein in complex with the PD of ACE2. Catalytic residues of ACE2, glycans, and Zn 2+ and Cl − ions are shown in brown, red, yellow, and purple, respectively. Å cushions in other directions. This puts a 50 Å water cushion between the RBD−PD complex and its periodic image in the x-direction, creating enough space for unbinding simulations. Ions were added to neutralize the system and the salt concentration was set to 150 mM to construct a physiologically relevant environment. The size of each solvated system was ∼164 000 atoms. All system preparation steps were performed in VMD. 38 cMD Simulations. All MD simulations were performed in NAMD 2.13 39 using the CHARMM36 40 force field with a time step of 2 fs. MD simulations were performed under N, P, T conditions. The temperature was maintained at 310 K using Langevin dynamics with a damping coefficient of 1 ps −1 . The pressure was maintained at 1 atm using the Langevin Nose− Hoover method with an oscillation period of 100 fs and a damping time scale of 50 fs. Periodic boundary conditions were applied. A cutoff distance of 12 Å was used for van der Waals interactions. Long-range electrostatic interactions were calculated using the particle-mesh Ewald method. For each system; first, 10 000 steps of minimization followed by 2 ns of equilibration were performed by keeping the protein fixed. The complete system was minimized for additional 10 000 steps, followed by 4 ns of equilibration by applying constraints on C α atoms. Subsequently, these constraints were released and the system was equilibrated for an additional 4 ns before initiating the production runs. The length of the equilibrium steps is expected to account for the structural differences due to the radically different thermodynamic conditions of crystallization solutions and MD simulations. 41 MD simulations were performed in supercomputers Comet and Stampede2 using ∼12 million core-hours in total.

Root Mean Square Fluctuation (RMSF) Calculations. RMSF values were calculated as ⟨ΔR i 2 ⟩ 1/2 = ⟨(R i − ⟨R i ⟩) 2 ⟩ 1/2 , where, ⟨R i ⟩ is the mean atomic coordinate of the ith C α atom and R i is its instantaneous coordinate.

SMD Simulations. SMD 42 simulations were used to explore the unbinding process of RBD from ACE2 on time scales accessible to standard simulation lengths. SMD simulations have been applied to explore a wide range of processes, including domain motion, 5,43 molecule unbinding, 44 and protein unfolding. 45 In SMD simulations, a dummy atom is attached to the center of mass of "steered" atoms via a virtual spring and pulled at a constant velocity along the "pulling direction", resulting in force F to be applied to the steered atoms along the pulling vector 39

where U is the guiding potential, k is the spring constant, v is the pulling velocity, t is the time, and R and R 0 are the coordinates of the center of mass of steered atoms at time t and 0, respectively, and n is the direction of pulling. 39 Total work (W) performed for each simulation was evaluated by integrating F over displacement ξ along the pulling direction as ∫ 0 ξ F(ξ)dξ. In SMD simulations of SARS-CoV-2, C α atoms of ACE2 residues S19-S43, T78-P84, Q325-N330, G352-I358, and P389-R393 were kept fixed, whereas C α atoms of RBD residues K417-I418, G446-F456, Y473-A475, and N487-Y505 were steered. Steered atoms were selected as the region comprising the interacting residues. For SARS-CoV SMD simulations the same ACE2 residues were kept fixed. However, two slightly different steered atom selections were applied: (i) using the same residue positions as for SARS-CoV-2, which are V404-I405, T433-L443, F460-S461, and N473-Y491, and (ii) selecting the region comprising the interacting residues, which are T433-L443, F460-D463, and N473-Y491. The total number of fixed and steered atoms was identical in all simulations. The pulling direction was selected as the vector pointing from the center of mass of fixed atoms to the center of mass of steered atoms. The pulling direction also serves as the reaction coordinate ξ for free energy calculations. SMD simulations were performed for 22.5 ns using a pulling velocity of 2 Å/ns (660 SMD simulations totaling 14850 ns in length), whereas 6 SMD simulations were performed for 300 ns each using a 0.1 Å/ns pulling velocity. To select the spring constant that satisfies the stiff-spring approximation, 31 we performed SMD simulations (v = 2 Å/ns) with various spring constants in the range of 25−125 kcal/(mol Å 2 ) ( Figure S1 ). At a spring constant of 100 kcal/(mol Å 2 ), the center of mass of the steered atoms followed the dummy atom closely while the spring was still soft enough to allow small deviations, hence satisfying the stiff-spring approximation.

For each system, 20 conformations were sampled with a 10 ns frequency from their cMD simulations (10 conformers from each set of the cMD simulations listed in Table S1 , MD1− 33a,b). These conformations served as 20 separate starting conformations, R 0 , for each set of SMD simulations (Table S1 , MD1−33c,d). In the SMD simulations deformation of the RBD structure was not observed; change in root-mean-square deviation (RMSD) was 1.0 ± 0.2 Å ( Figure S2 ).

Potential of Mean Force for Unbinding of RBD. Work values to unbind RBD from ACE2 at low pulling velocities along the reaction coordinate were analyzed using Jarzynski equality, which provides a relation between equilibrium free energy differences and the work performed through nonequilibrium processes 31, 46, 47 

where ΔA is the change in the free energy, k B is the Boltzmann constant, and T is the temperature. Because work values sampled in our SMD simulations differ more than 1 k B T, the average work calculated in eq 3 will be dominated by small work values that are only rarely sampled. For a finite (N)

Thus, eq 3 provides an upper bound on ΔA, which was used as an estimate of the potential of mean force (PMF). 31 

Interaction Sites between the S Protein and ACE2. To model the dynamic interactions of the S-ACE2 binding interface, we used the co-structure of RBD of the SARS-CoV-2 S protein in complex with the PD of human ACE2 8 (Figure 1c) . The structure was solvated in a water box that contains a physiologically relevant salt (150 mM NaCl) concentration. Two sets of cMD simulations, each 100 ns in length (Table S1) , were performed to determine the formation of salt bridges 48 and hydrogen bonds, as well as electrostatic and hydrophobic interactions between RBD and PD. A cutoff distance of 6 Å between the basic nitrogens and acidic oxygens was used to score a salt bridge formation. 48 For hydrogen bond formation, a maximum distance of 3.5 Å between hydrogen

The Journal of Physical Chemistry B pubs.acs.org/JPCB Article bond donor and acceptor and a 30°angle between the hydrogen atom, the donor heavy atom, and the acceptor heavy atom was used. 49 Interaction pairs that satisfy the distance, but not the angle criteria were analyzed as electrostatic interactions. For hydrophobic interactions, a cutoff distance of 8 Å between the side chain carbon atoms was used. 50−52 Using these criteria, we identified eleven hydrophobic interactions (Figure 2a) , eight hydrogen bonds (Figure 2b) , two salt bridges, and six electrostatic interactions (Figure 2c ) between RBD and PD. Observation frequencies were classified as high and moderate for interactions that occur in 49% and above and between 15 and 48% of the total trajectory, respectively. F486 and Y489 of RBD formed hydrophobic interactions with F28, L79, M82, and Y83 of PD, while L455, F456, Y473, and A475 of RBD formed hydrophobic interactions with T27 of PD at high frequencies ( Figure 2d ).

Salt bridges between K417-D30 (RBD−PD) and E484-K31, and hydrogen bonds among N487-Y83, T500-D355, and Q493-E35 were observed at high frequencies, whereas hydrogen bonds among Y449-D38, Q498-K353, T500-Y41, Y505-E37, and Q493-E35 were observed at moderate frequencies ( Figure 2d ). Residue pairs Y453-H34, N487-Q24, T500-Y41, N501-K353, Q493-K31, and Y449-Q42 exhibited electrostatic interactions throughout the simulations (Figure 2d ). The interaction network we identified in our cMD simulations was mostly consistent with reported interactions in the RBD−PD crystal structure. 8 However, our simulations identified four hydrogen bonds (Q498-K353, T500-D355, Y505-E37, and Q498-Q42), one hydrophobic interaction (L455-T27), and two electrostatic interactions (Y453-H34 and N501-K353) that are not present in the crystal structure. In turn, we did not detect frequent hydrogen bonding between G446-Q42, G502-K353, and Y505-R393 and electrostatic interaction of G496-K353 observed in the crystal structure. 8 This discrepancy may be due to radically different thermodynamic conditions between crystallization solutions and cMD simulations. 41 Our results are more consistent with recent cryo-EM studies, 7,32 which comprised all of the missing interactions, except Y505-E37. We divided the RBD−PD interaction surface into three regions (CR1−3, Figure 2a −c). 27 CR2 comprised significantly fewer interactions than the ends of the RBD binding surface (CR1 and CR3). Remarkably, 10 out of 13 interactions we detected in CR1 were hydrophobic, which were proposed to play a central role in the anchoring of RBD to PD. 27 Unlike CR1, CR2 formed only a single hydrophobic interaction with PD, whereas CR3 did not form any hydrophobic interactions.

Unbinding of the S Protein from ACE2 under Force. To estimate the binding strength of the S protein to ACE2, we performed SMD simulations to pull RBD away from PD at a constant velocity of 2 Å/ns along the vector pointing away from the binding interface (Figure 3a) . Steering forces were applied to the C α atoms of the RBD residues on the binding interface, whereas C α atoms of PD residues at the binding interface were kept fixed. Simulations at a pulling velocity of 2 Å/ns were also repeated in the absence of ACE2 to account for the work done against the viscous drag of water ( Figure S3 ). The calculated average work against viscous drag was subtracted from all SMD work values. In 20 SMD simulations (22.5 ns each, a total of 450 ns in length, Table S1), the average work applied to unbind SARS-CoV-2 RBD from PD was 64.4 ± 5.4 kcal/mol (mean ± s.d.). We also used the Jarzynski equality 46,47 to estimate the free energy profiles as a function of a reaction coordinate, referred to as PMF. 31 Binding free energy was estimated as −55.5 kcal/mol based on 20 SMD simulations performed at 2 Å/ns. Therefore, our SMD simulations demonstrate that the S protein binds stably to ACE2 (Figure 3b ).

Because part of the work applied is lost to the irreversible processes as we pull RBD away from PD at a finite velocity, our simulations provide the upper bound estimate of the free energy of S-ACE2 binding. To obtain a closer estimate of the unbinding work, we lowered the SMD pulling speed to 0.1 Å/ The Journal of Physical Chemistry B pubs.acs.org/JPCB Article ns, which is comparable to those applied experimentally in high-speed atomic force microscopy. 53 We performed 3 SMD simulations at this velocity (300 ns each, a total of 900 ns in length, Table S1 , 34c,d). The average work of these trajectories was 36.8 ± 9.4 kcal/mol (mean ± s.d.).

Mutagenesis of the S-ACE2 Binding Interface. To investigate the contribution of the interacting residues to the overall binding strength, we introduced point mutations on the RBD. Salt bridges were eliminated by charge reversals (K417E and E484K). We also replaced each interacting amino acid with alanine (except A475, Table S1) to disrupt the pairwise interactions, 54 with minimal perturbations to the protein backbone. 55 Two sets of cMD simulations (a total of 3.4 μs in length) were performed for each point mutant. We first quantified the RMSF of the C α atom of the RBD residues located on the PD binding surface (Figure 3c ). The rigid body motions were eliminated by aligning the RBD interacting surface of PD for each conformer (see Methods). Thirteen out of 17 mutations increased the residue fluctuations compared to WT ( Figure S4a ), suggesting that disrupting the interactions between RBD and PD results in floppier binding. The largest fluctuations were observed for 2 mutations in CR1 (F486A and N487A), 2 mutations in CR3 (Y449A and Y505A), and 1 mutation in CR2 (L455A) (Figure 3c ). Mutation of these residues also increased the fluctuations in their neighboring region. While mutations in CR1 increased fluctuations in CR3 significantly, mutations in CR3 had little to no effect on the fluctuations in CR1 (Figures 3d and S4b) .

We next performed SMD simulations at 2 Å/ns pulling speed to model unbinding of each point mutant from PD (20 simulations for each mutant, a total of 7.65 μs in length, Table  S1 ) and provide relative changes in the binding free energy of WT and mutant RBD under the same velocity and thermodynamic conditions. F486A, Y489A, E484A, E484K, Y505A, and N487A, mutations decreased the work require-ment to unbind RBD−PD by 6−13% (Figures 3e, f and S5 ). Based on the estimated PMF using the Jarzynski equality ( Figure S6 ), these mutations resulted in a decrease in the binding energy by 10−27% compared to WT. We note that most of these mutations also led to the largest increase in residue fluctuations on the binding surface (Figure 3c ). Five of these mutations (F486A, Y489A, E484A, E484K, and N487A) are located in CR1, whereas Y505A is located in CR3. These results highlight the primary role of hydrophobic interactions in CR1 to stabilize S-ACE2 binding.

To further characterize critical interactions of the S-ACE2 binding interface, we introduced double mutants to neighboring residues of RBD that form critical interactions with PD. We performed a total of 2.8 μs of cMD and 6.3 μs of SMD simulations for 14 double mutants (Table S1 ). In particular, double mutants in CR1 resulted in 4 out of 6 highest increases in RMSF (Figures 4a and S4a) . The F486A/N487A mutation at CR1 resulted in the largest increase in fluctuations in both CR1 and CR3 (Figures 4a and S4b) . In SMD simulations, 12 out of 14 double mutations also further decreased the average work to unbind RBD from PD (Figures 4c,d and S7 ). (Figures 4d and  S8) . A charge reversal of K417E in combination with either Q493A or Y453A also resulted in a large decrease in work values (Figure 4d) . Collectively, these results show that two salt bridges (E484-K31 and K417-D30) and the network of 27 we investigated the order of events that result in detachment of RBD from PD in SMD simulations. The unbinding process appears to perform a zipper-like detachment starting from CR3 and ending at CR1 in 80% of the simulations when RBD was pulled at 2 Å/ns pulling velocity (Figure 5a and Movie S1). In only 20% of the simulations, CR3 was released last from PD (Figure 5a and Movie S2). Similarly, CR1 released last in 2 out of the 3 SMD simulations performed using 0.1 Å/ns pulling velocity. Because unbinding simulations can reveal features characteristic of the reverse process of binding, 56−60 these results suggest that CR1 binding is the first and critical event for the S protein binding to ACE2. Mutagenesis of the critical residues in CR1, in general, resulted in a substantial decrease in the percentages of unbinding events that terminate with the release of CR1 from PD. In alanine replacement of the hydrophobic residues (F456A, Y473A, F486A, and Y489A), CR1 was released last in 75, 80, 45, and 75% of the SMD simulations, respectively (Figure 5b) . The probability of CR1 to release last under force was further reduced in double mutants of F486A/N487A (65%), E484A/ Y489A (55%), E484A/F486A (50%), and L455A/F456A (60%) (Figure 5b ). Unlike these mutants, E484A and E484K mutants in CR1 increased the probability of CR1 to release last. These results indicate that single and double mutants of the critical residues in CR1 substantially reduce the binding free energy of this region to ACE2.

Comparison of the ACE2 Binding Strengths of SARS-CoV-2 and SARS-CoV S Proteins. It remains unclear whether higher infectivity of SARS-CoV-2 than SARS-CoV can be attributed to stronger interactions between S and ACE2 in SARS-CoV-2. 2, 19 There are both computational 23, 26, 29 and experimental 2 studies which reported similar binding strengths for SARS-CoV and SARS-CoV-2 S proteins to ACE2, whereas others have reported that SARS-CoV-2 S protein interacts more strongly with ACE2 than the SARS-CoV S protein. 28, 34, 61 In addition, using protein pull-down assays it was reported that SARS-CoV-2 RBD shows stronger ACE2 binding affinity than the SARS-CoV RBD, whereas the entire SARS-CoV-2 S protein shows ACE2 binding affinity comparable to or lower than the entire SARS-CoV-2 S protein. 35 To test whether SARS-CoV and SARS-CoV-2 S proteins have distinguishable binding strengths, we performed two sets of MD simulations for the RBD of SARS-CoV S protein bound to the PD of ACE2 (PDB ID: 2AJF 9 ) and compared these results to those of SARS-CoV-2. Similar to SARS-CoV-2, RBD of SARS-CoV makes an extensive network of interactions with PD. We identified eleven hydrophobic interactions (Figure 6a ), six hydrogen bonds (Figure 6b) , and seven electrostatic interactions (Figure 6c ). Only 6 of these interacting amino acids are conserved in SARS-CoV-2 and the following mutations have taken place: L443/F456 (SARS-CoV/SARS-CoV-2), F460/Y473, P462/A475, P470/E484, L472/F486, V404/K417, N479/Q493, Y484/Q498, and T487/N501. Similar to SARS-CoV-2, L472 and Y475 of SARS-CoV RBD formed a total of seven hydrophobic interactions at a high frequency with the hydrophobic pocket of ACE2 (Figure 6d ). Unlike SARS-CoV-2, SARS-CoV RBD did not form any salt bridges with ACE2.

We next modeled the unbinding of RBD of SARS-CoV from PD by performing 20 SMD simulations with a pulling velocity of 2 Å/ns (a total of 450 ns in length, Table S1 ). The average work applied to unbind SARS-CoV RBD from PD was 68.4 ± 8.7 kcal/mol (mean ± s.d., Figure 6e ). Furthermore, we modeled SARS-CoV RBD unbinding with a pulling velocity of 0.1 Å/ns (a total of 900 ns in length, Table S1 , 35c,d), which resulted in an average work of 37.1 ± 2.4 kcal/mol (mean ± s.d.) (Figure 6e ). These average unbinding work values are indistinguishable from that of SARS-CoV-2 under the same pulling speed conditions (two-tailed Student's t-test, p = 0.99 and 0.09 for 2 and 0.1 Å/ns pulling velocities, respectively). Using the work values for SMD pulling velocity of 2 Å/ns and the Jarzynski equality, binding free energy was calculated as −54.5 kcal/mol for SARS-CoV. Unlike SARS-CoV-2, CR1 released last from PD in only 50% of the unbinding events of RBD of SARS-CoV, whereas the unbinding of CR3 was the last event in the remaining 50% (Figure 6f) . These results indicate that the S protein binds stably to ACE2 in both SARS-CoV and SARS-CoV-2 with similar binding strengths. The absence 

We performed an extensive set of in silico analysis to identify critical residues that facilitate binding of the RBD of the SARS-CoV-2 S protein to the human ACE2 receptor. Mutagenesis of these residues and pulling the RBD away from PD at a low velocity enabled us to estimate the free energy of binding and the order of events that result in the unbinding of RBD from PD. There is currently no consensus on the exact binding free energy of the S protein to ACE2. Binding free energies ranging from −23.19 up to −140.0 kcal/mol were reported by postprocessing all-atom MD trajectories (generalized Born and surface area continuum solvation approach (MM-GBSA) 19, 25, 26, 34 and molecular mechanics/Poisson−Boltzmann surface area (MM-PBSA) 20, 30 ) and 12.6 kcal/mol by umbrella sampling simulations in which the protein backbones were kept fixed in addition to the applied umbrella potential. 62 Furthermore, the average unbinding work was evaluated 28 as 135.0 ± 4.9 kcal/mol using 5 SMD simulations, which were performed by pulling the center of mass of RBD at a pulling Our simulations showed that the PD interacting surface of RBD can be divided into three contact regions (CR1−3). Hydrophobic residues of CR1 strongly interact with the hydrophobic pocket of PD in both SARS-CoV and SARS-CoV-2. CR1 of SARS-CoV-2 also forms a salt bridge with ACE2 that is not present in SARS-CoV. Based on our SMD simulations, we did not observe a major difference in the binding strength of the S protein to ACE2 between SARS-CoV and SARS-CoV-2. These results are consistent with the recent MD simulations, 26, 29 coarse-grained simulations, 23 and biolayer interferometry, 2 while inconsistent with other MM-GBSA studies, 19, 34 and also an SMD study 28 which reported different binding strengths for SARS-CoV-2 and SARS-CoV S proteins. These differences may be attributed to different selection of simulation parameters and force fields, sampling sizes, simulation lengths, number of replicas, the presence of glycans and explicit solvent, ion concentrations, applied external constraints, and free energy calculation methods. Furthermore, an experimental study reported that SARS-CoV-2 RBD has higher ACE2 binding affinity than SARS-CoV RBD, while its entire spike has ACE2 binding affinity comparable to or lower than that of SARS-CoV. 35 Collectively, there is consensus that S proteins of both SARS-CoV and SARS-CoV-2 bind tightly to ACE2, but it remains controversial whether SARS-CoV-2 S binds more tightly to ACE2 and whether this increases the infectivity of SARS-CoV-2.

Our analysis suggests that CR1 is the main anchor site of the SARS-CoV-2 S protein for binding with ACE2, and blocking the CR1 residues F456, Y473, E484, F486, N487, and Y489 could significantly reduce the binding affinity. Consistent with this prediction, llama-based nanobodies H11-H4 and H11-D4 neutralize SARS-CoV-2 12 by interacting with 50% of the critical residues we identified in CR1 and CR2. Furthermore, alpaca-based nanobody Ty1 neutralizes SARS-CoV-2 66 and among its primary interactions, E484 in CR1, and Q493 and Y449 in CR3 were also determined as critical residues in our study. Similarly, the human neutralizing antibodies CV30, 67 B38, 67 CB6, 67 and VH3−53 68,69 interact with 50−100% of the critical residues we identified. Starr et al. 33 performed a deep mutational scanning of SARS-CoV-2 RBD amino acids using flow cytometry and demonstrated that RBD residues Y449, L455, F486, and Y505 are required for ACE2 binding. Two of these residues (Y505 and F486) were determined as critical in our study. Their mutagenesis results also showed that mutations in Q493, Q498, and N501 residues enhanced the affinity to ACE2, consistent with our alanine mutagenesis results.

Experimental studies revealed that antibodies against SARS-CoV induce limited neutralizing activity against SARS-CoV-2. 11, 27 This may be attributed to the low sequence conservation of the CR1 between SARS-CoV and SARS-CoV-2. In particular, the S protein of SARS-CoV-2 contains critical phenylalanine (F486) and glutamate (E484) residues not present in SARS-CoV, that form hydrophobic interactions and a salt bridge with ACE2, respectively. It remains to be determined whether this difference plays a role in higher infectivity of SARS-CoV-2 than SARS-CoV.

Our simulations show that single and double mutants of CR1 are not sufficient to disrupt the binding of RBD to ACE2, but reduce the binding free energy of this region. Because RBD makes multiple contacts with ACE2 through an extended surface, small molecules or peptides that target a specific region in the RBD−ACE2 interaction surface may not be sufficient to prevent binding of the S protein to ACE2. Instead, blocking a larger surface of the CR1 with a neutralizing antibody or nanobody is more likely to introduce steric constraints to prevent S-ACE2 interactions. (Table  S1 ) (PDF) CR1 in 80% of the simulations when RBD was pulled at 2 Å/ns pulling velocity (Movie S1) (AVI) 20% of the simulations, CR3 was released last from PD (Movie S2) (AVI)

contact region; fs, femtosecond; MD, molecular dynamics

nanoscale molecular dynamics; ns, nanosecond; PD, peptidase domain

PMF, potential of mean force; ps, picosecond; RBD, receptor-binding domain; RMSD, rootmean-square deviation

RMSF, root mean square fluctuation

SMD, steered molecular dynamics; TMPRSS2, transmembrane serine protease 2

VMD, visual molecular dynamics; WT, wild-type ■ REFERENCES

Identification of a novel coronavirus causing severe pneumonia in human: a descriptive study

Structure, function, and antigenicity of the SARS-CoV-2 spike glycoprotein

Mechanisms of coronavirus cell entry mediated by the viral spike protein

Stabilized coronavirus spikes are resistant to conformational changes induced by receptor recognition or proteolysis

Conformational transition of SARS-CoV-2 spike glycoprotein between its closed and open states

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

Structural basis for the recognition of the SARS-CoV-2 by full-length human ACE2

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

Structure of SARS coronavirus spike receptor-binding domain complexed with receptor

SARS-CoV-2 cell entry depends on ACE2 and TMPRSS2 and is blocked by a clinically proven protease inhibitor

Key residues of the receptor binding motif in the spike protein of SARS-CoV-2 that interact with ACE2 and neutralizing antibodies

Computational Design of ACE2-Based Peptide Inhibitors of SARS-CoV-2

De novo design of picomolar SARS-CoV-2 miniprotein inhibitors

Targeting the SARS-CoV-2-spike protein: from antibodies to miniproteins and peptides

Repurposing approved drugs as inhibitors of SARS-CoV-2 S-protein from molecular modeling and virtual screening

High-throughput virtual screening of drug databanks for potential inhibitors of SARS-CoV-2 spike glycoprotein

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

Molecular Mechanism of Evolution and Human Infection with SARS-CoV-2. Viruses

SARS-CoV-2, an evolutionary perspective of interaction with human ACE2 reveals undiscovered amino acids necessary for complex stability

Computational Prediction of Mutational Effects on SARS-CoV-2 Binding by Relative Free Energy Calculations

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

Critical Differences Between the Binding Features of the Spike Proteins of SARS-CoV-2 and SARS-CoV

Effect of mutation on structure, function and dynamics of receptor binding domain of human SARS-CoV-2 with host cell receptor ACE2: a molecular dynamics simulations study

Dynamics of the ACE2−SARS-CoV-2/ SARS-CoV spike protein interface reveal unique mechanisms

The MERS-CoV receptor DPP4 as a candidate binding target of the SARS-CoV-2 spike

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

Does SARS-CoV-2 Bind to Human ACE2 More Strongly Than Does SARS-CoV?

Comparing the binding interactions in the receptor binding domains of SARS-CoV-2 and SARS-CoV

Critical Sequence Hotspots for Binding of Novel Coronavirus to Angiotensin Converter Enzyme as Evaluated by Molecular Simulations

Free energy calculation from steered molecular dynamics simulations using Jarzynski's equality

Receptor binding and priming of the spike protein of SARS-CoV-2 for membrane fusion

Deep mutational scanning of SARS-CoV-2 receptor binding domain reveals constraints on folding and ACE2 binding

Is the rigidity of SARS-CoV-2 spike receptor-binding motif the hallmark for its enhanced infectivity? Insights from all-atom simulations

Cell entry mechanisms of SARS-CoV-2

Improved treatment of ligands and coupling effects in empirical calculation and rationalization of p K a values

PROPKA3: consistent treatment of internal and surface residues in empirical p K a predictions

VMD: visual molecular dynamics

Scalable molecular dynamics with NAMD

Jr Optimization of the additive CHARMM all-atom protein force field targeting improved sampling of the backbone ϕ, ψ and side-chain χ1 and χ2 dihedral angles

Why protein conformers in molecular dynamics simulations differ from their crystal structures: a thermodynamic insight

Steered molecular dynamics simulation of the Rieske subunit motion in the cytochrome bc1 complex

Computational design of new peptide inhibitors for amyloid beta (Aβ) aggregation in Alzheimer's disease: application of a novel methodology

Unfolding of titin immunoglobulin domains by steered molecular dynamics simulation

Equilibrium free-energy differences from nonequilibrium measurements: A master-equation approach

Nonequilibrium equality for free energy differences

Zipping and unzipping of adenylate kinase: atomistic insights into the ensemble of open↔ closed transitions

HBonanza: a computer algorithm for molecular-dynamics-trajectory hydrogen-bond analysis

Direct and quantitative AFM measurements of the concentration and temperature dependence of the hydrophobic force law at nanoscopic contacts

A study of the preferred environment of amino acid residues in globular proteins

Molecular Dynamics Simulation of Antimicrobial Peptide Arenicin-2: b-Hairpin Stabilization by Noncovalent Interactions

High-speed force spectroscopy unfolds titin at the velocity of molecular dynamics simulations

Rapid mapping of protein functional epitopes by combinatorial alanine scanning

Comparing experimental and computational alanine scanning techniques for probing a prototypical protein−protein interaction

Molecular dynamics simulations suggest that electrostatic funnel directs binding of Tamiflu to influenza N1 neuraminidases

Molecular dynamics study of unbinding of the avidin-biotin complex

Steered molecular dynamics simulations reveal the likelier dissociation pathway of imatinib from its targeting kinases c-Kit and Abl

Unbinding of nicotine from the acetylcholine binding protein: steered molecular dynamics simulations

Computational insights into the mechanism of ligand unbinding and selectivity of estrogen receptors

Cryo-EM structure of the 2019-nCoV spike in the prefusion conformation

Quantifying the adhesive strength between the SARS-CoV-2 S-proteins and human receptor and its effect in therapeutics

Alsteens, D. Molecular interaction and inhibition of SARS-CoV-2 binding to the ACE2 receptor

Structural and functional basis of SARS-CoV-2 entry by using human ACE2

Potent binding of 2019 novel coronavirus spike protein by a SARS coronavirus-specific human monoclonal antibody

An alpaca nanobody neutralizes SARS-CoV-2 by blocking receptor interaction

Structural basis for potent neutralization of SARS-CoV-2 and role of antibody affinity maturation

An alternative binding mode of IGHV3-53 antibodies to the SARS-CoV-2 receptor binding domain