Molecular Dynamics Reveals Complex Compensatory Effects of Ionic Strength on the Severe Acute Respiratory Syndrome Coronavirus 2 Spike/Human Angiotensin-Converting Enzyme 2 Interaction

The SARS-CoV-2 pandemic has already killed more than one million people worldwide. To gain entry, the virus uses its Spike protein to bind to host hACE-2 receptors on the host cell surface and mediate fusion between viral and cell membranes. As initial steps leading to virus entry involve significant changes in protein conformation as well as in the electrostatic environment in the vicinity of the Spike/hACE-2 complex, we explored the sensitivity of the interaction to changes in ionic strength through computational simulations and surface plasmon resonance. We identified two regions in the receptor-binding domain (RBD), E1 and E2, which interact differently with hACE-2. At high salt concentration, E2-mediated interactions are weakened but are compensated by strengthening E1-mediated hydrophobic interactions. These results provide a detailed molecular understanding of Spike RBD/hACE-2 complex formation and stability under a wide range of ionic strengths.


Molecular dynamics (MD) simulations
The initial coordinates for RBD/hACE-2 complex were obtained from the Protein Data Bank (PDB ID 6M0J, resolution 2.45 Å) 1 and prepared using the UCSF Chimera 2 by removing co-crystallized hetero groups (except Zn 2+ ) and water molecules. The system conditions were prepared using GROMACS v. 5.1.4 [3][4][5][6] with the OPLS-AA force field 7 . We used the Swiss-Param web-based service to build the Zn 2+ topology 8 . All systems were explicitly solvated with TIP3P water models in a triclinic box (81.06 X 91.06 X 135.02 Å³), neutralized keeping NaCl concentration of 150 mM (114 Na + and 90 Clions), 500 mM (324 Na + and 300 Clions) and 1,000 mM (624 Na + and 600 Clions) and minimized until a maximum force of 10.0 kJ.mol -1 or 5,000 steps. The systems were consecutively equilibrated in isothermal-isochoric (NVT) and isothermal-isobaric (1 bar; NpT) ensembles at 310 K for 1000 ps. MD simulations in the NpT ensemble were carried out for 100 ns in a periodic cubic box considering the minimum distance of 10 Å between any protein atom and box walls. The Zn 2+ positions were restricted throughout the simulations. Distance calculations, root-mean-square deviation (RMSD), and root-mean-square fluctuation (RMSF) calculations were performed using GROMACS. The opening degree of the RBD/hACE-2 interface was quantified as the angle between two vectors chosen using the Cα coordinates along ACE-2 α-helix (residues 38 -49), and the Spike RBD β-strand (residues 492 -494). Analysis of Cα cross-correlated displacements were performed using the R-based package Bio3d 9 . The covariances ranged from -1 (anti-correlated motions) to 1 (correlated motions). The Cα distances were calculated using GROMACS.

Hydrophobicity scale calculations
ClustalW/X 10 was used to compute multiple amino acid sequence alignments. The hydropathy indexes were calculated using Kyte and Doolitle 11 .

Induced dissociation, free diffusion association and Umbrella Sampling simulation
The induced dissociation, free diffusion association, and Umbrella Sampling simulations were done with GROMACS 2018.6 3-6 using the RBD/hACE-2 complex structure (PDB ID 6M0J, resolution 2.45 Å) 1 as initial coordinates and the CHARMM36 force field [12][13][14] . The systems were neutralized keeping the NaCl concentration at 150 mM (223 Na + and 199 Clions), 500 mM (650 Na + and 626 Clions) or 1,000 mM (1276 Na + and 1252 Clions) in a triclinic box (100 X 100 X 210 Å³) filled with 67,980, 62,954, and 61,752 TIP3P waters, respectively. The energy minimization was set to a maximum force of 1,000 kJ mol -1 nm -1 and a maximum number of steps of 50,000. The systems were equilibrated consecutively in isothermal-isochoric (NVT) and isothermal-isobaric (1 bar; NpT) ensembles at 310 K for 2000 ps, followed by the calculation of a 100 ns MD trajectory in the NpT ensemble.

Induced dissociation
Starting from the last frame of 100 ns MD trajectory, a reaction coordinate was calculated as the average of the Cα (RBD) -Cα (hACE-2) intermolecular vectors of all residue pairs at the interface 1 . The intermolecular vectors were computed using an inhouse python script. The RBD/hACE-2 complex was rotated to align the z-axis of the structure with the reaction coordinate vector. A harmonic potential of 600 kJ mol −1 .nm −2 constant force and 0.005 nm.ns −1 pulling speed was applied on the center of mass of the RBD, while the positions of the hACE-2 atoms were restricted using a harmonic force of 1000 kJ.mol −1 .nm −2 . The pulling trajectory time was 1 ns.

Free diffusion simulation
The RBD and hACE-2 were pulled as described above up to a center-of-mass distance of ∼60 Å. The separation distance should be large to capture diffusion effects, but not as large to avoid interactions between them. To equilibrate the system, it was relaxed during a 3 ns annealing, in which the temperature was increased from 310 K to 350 K in 500 ps, followed by an NpT equilibrium dynamics at 350 K for 2 ns, and a S5 temperature decrease from 350 K to 310 K in 500 ps. The two molecules were then released of all restriction forces, and a MD simulation of 20 ns was performed. When this experiment was performed at 500 mM NaCl, there was no recovery of the RBD/hACE-2 complex after 20 ns of simulation. The analysis was based on the 2,000 snapshots saved every 10 ps during the trajectory. The number of Na + and Cl − ions within a radius of 6 Å from each amino acid at the binding interface was computed using an in-house python script.

Umbrella Sampling
Initial configurations for a biased sampling simulation were recorded at every step of ~ 0.1 nm in the dissociation pathway. Umbrella sampling (US) windows were equilibrated in the NpT ensemble using the simulated annealing protocol described above, followed by 10 ns simulation. In total, 33 configurations were calculated at 150 and 500 mM NaCl concentrations, while 30 configurations were computed at 1,000 mM NaCl solution. The three additional configurations at 150 and 500 mM NaCl concentrations generated overlapping biased distribution functions. The weighted histogram analysis method (WHAM) 15 was used to compute the potential of mean force (PMF) along the reaction coordinate (ξ) with the g_wham GROMACS module.
Free binding energies (ΔG) were estimated from the difference between the minimum and maximum of the PMF profile in a given NaCl concentration. Dissociation constants (KD) were determined assuming = ( ), R and T are gas constant and temperature, respectively. PMF errors were calculated using 100 rounds of bootstrapping analysis.

Protein-protein binding assays using surface plasmon resonance (SPR)
The experimental KD (M -1 ) of the spike/hACE-2 complex were obtained by surface plasmon resonance (SPR) assays using a Biacore 8K system (GE Healthcare). The assays as well as protein cloning, expression and purification were performed by the μL.min -1 during 120 s. Resulting data were fitted to a 1:1 binding model using Biacore Evaluation Software (GE Healthcare).
Multiple amino acid sequence alignment of Spike S1 subunit within coronavirus members is represented by WebLogo. The sequence shown below the logo belongs to SARS-CoV-2. The residues in bold belong to RBD and residues highlighted in S9 yellow interact with the human ACE-2 receptor 1 . The cysteines colored in red form a thiosulfate bond in the E1 region probably to maintain the 7-8 hairpin stable in the three-dimensional structure (shown in the Figure  We also used the sequence of bat RaTG13 (GenBank: QHR63300.2) and pangolin coronaviruses (GenBank: QIA48632.1). The multiple sequence alignment was created using ClustalW 16 . The sequence logo was created with the WebLogo server 17 . The residues belonging to the RBD are inside of the dotted blue square.   The two distances that were monitored during the simulations are indicated: G476 RBD -S19 hACE-2 in E1 and G502 RBD -T52 hACE-2 in E2. b) The initial S19 hACE-2 -G476 RBD distance of ~5-6 Å displays a longer lag phase at 1,000 mM when compared to 150 mM NaCl, a difference of ~0.3 ns. At 500 mM, the initial distance is ~13 Å, and the lag phase is significantly shorter,    Table 2 in the main text and Table S9. At least five different RBD concentrations were used to calculate the KD values. Figure S10. Surface plasmon resonance (SPR) data using immobilized recombinant dimer hACE-2 and different concentrations of purified recombinant trimeric form of the Spike (residues 16 to 1213, with mutations at arginines 683 and 685 to alanines: R683A, R685A), ranging from 2 to 500 nM. The experimental data were fitted to a 1:1 S18 binding model using Biacore Evaluation Software (GE Healthcare) to calculate the kinetics parameters (ka and kd) shown in Table 2 in the main text and Table S9. At least six different trimeric form of the Spike concentrations were used to calculate the KD values. Figure S11. Surface plasmon resonance (SPR) data using immobilized recombinant dimer hACE-2 and different concentrations of purified recombinant S1 subunit of Spike (residues residues 16 to 685), ranging from 4 to 250 nM. The experimental data were fitted to a 1:1 binding model using Biacore Evaluation Software (GE Healthcare) to calculate the kinetics parameters (ka and kd) shown in Table 2 in the main text and Table S9. At least five different spike S1 subunit concentrations were used to calculate the KD values. Figure S12. Surface plasmon resonance (SPR) data using immobilized recombinant dimer hACE-2 and different concentrations of purified recombinant monomeric form of the Spike (residues 16 to 1213, with mutations at arginines 683 and 685 to alanines: R683A, R685A), ranging from 4 to 25 nM. The experimental data were fitted to a 1:1 S19 binding model using Biacore Evaluation Software (GE Healthcare) to calculate the kinetics parameters (ka and kd) shown in Table 2 in the main text and Table S9. Seven different Spike concentrations were used to calculate the KD values. Table S1.
The hydropathy scale calculated from Kyte and Doolittle (1982) 11 for each residue of RBDs of SARS-CoV and SARS-CoV-2. Residues located in the E1 region are colored in green and those located in the E2 region in orange. Residues in asterisks make part of the hACE-2 interface.   positive  negative  positive  negative  positive  negative  ACE2 RBD ACE2 RBD ACE2 RBD ACE2 RBD ACE2 RBD ACE2