key: cord-0299249-g4ul6fuw
authors: Lan, Jun; He, Xinheng; Ren, Yifei; Wang, Ziyi; Zhou, Huan; Fan, Shilong; Zhu, Chenyou; Liu, Dongsheng; Shao, Bin; Liu, Tie-Yan; Wang, Qisheng; Zhang, Linqi; Ge, Jiwan; Wang, Tong; Wang, Xinquan
title: Structural and computational insights into the SARS-CoV-2 Omicron RBD-ACE2 interaction
date: 2022-01-04
journal: bioRxiv
DOI: 10.1101/2022.01.03.474855
sha: 5148f45b90aadfc8d04fe4698412c4605883b059
doc_id: 299249
cord_uid: g4ul6fuw

Since SARS-CoV-2 Omicron variant (B.1.1.529) was reported in November 2021, it has quickly spread to many countries and outcompeted the globally dominant Delta variant in several countries. The Omicron variant contains the largest number of mutations to date, with 32 mutations located at spike (S) glycoprotein, which raised great concern for its enhanced viral fitness and immune escape[1–4]. In this study, we reported the crystal structure of the receptor binding domain (RBD) of Omicron variant S glycoprotein bound to human ACE2 at a resolution of 2.6 Å. Structural comparison, molecular dynamics simulation and binding free energy calculation collectively identified four key mutations (S477N, G496S, Q498R and N501Y) for the enhanced binding of ACE2 by the Omicron RBD compared to the WT RBD. Representative states of the WT and Omicron RBD-ACE2 systems were identified by Markov State Model, which provides a dynamic explanation for the enhanced binding of Omicron RBD. The effects of the mutations in the RBD for antibody recognition were analyzed, especially for the S371L/S373P/S375F substitutions significantly changing the local conformation of the residing loop to deactivate several class IV neutralizing antibodies.

fight against the COVID-19 is still in great uncertainty due to the emergence of SARS-CoV-2 46 variants, especially the variants of concern (VOC) with changed pathogenicity, increased 47 and Extended Data Table 3 ). We also utilized Molecular Dynamics (MD) simulation and 114 Molecular Mechanics Generalized Born Surface Area (MM/GBSA) method to estimate the 115 binding free energies (Fig. 2d) . The binding free energies of the Omicron and WT RBDs to ACE2 116 are -42.19 ± 6.61 and -26.66 ± 9.78 kcal/mol, respectively, which also supports the enhanced 117 binding of the Omicron RBD to ACE2. We also analyzed the contribution of each mutation in the 118 binding free energy change (Fig. 2d) . Four mutations including S477N, G496S, Q498R and 119 N501Y have positive impacts on the binding of the Omicron RBD to ACE2 by further reducing 120 the free energy by 0.77, 0.36, 2.85 and 5.84 kcal/mol respectively, which is consistent with the 121 above described newly formed interactions of at these four sites, especially around the 496, 498 122 and 501 cluster. We also observed the negative impacts of Q493K/R and Y505H mutations. The 123 negative impact of the mutation at 493 site is a little unexpected because a new salt bridge to 124 ACE2 E35 is formed by Q493K/R mutation. The loss of Q493 to ACE2 K31 interaction and the 125 close position of positively charged ACE2 K31 and Omicron RBD K493 or R493 may exceed the 126 positive effect of the salt bridge to ACE2 E35 and finally determine the net negative impact of 127 Q493K/R mutation. 128

To explore the variations between the WT and Omicron RBD-ACE2 interactions from a 129 dynamic view, we run microsecond level of MD simulations for the WT and Omicron systems and 130 utilized Markov State Model (MSM) to identify representative states and illustrate the transition 131 process among them. Time-lagged independent component (tIC) analysis by projecting the 132 conformations onto the low-dimensional surface showed that the WT and Omicron RBD-ACE2 interface contacts have different distributions in the first two components, tIC 1 and tIC 2 (Fig.  134   3a) . With a novel lag time chosen algorithm recently designed by us, all conformations were first 135 clustered into microstates and then into macrostates, which are representative conformations for a 136 set of points on the energy landscape (Fig. 3a) . Finally, three macrostates 1, 2, and 3 were sorted 137 out according to their proportions in all conformations (Fig. 3a) . The state 3 of the WT and 138

Omicron locate at a similar position around (0, -1) on the energy landscape. Although the 139 positions of state 2 in WT and Omicron are also similar, the WT system shows a larger value (> 140 2.0) on the tIC 2 dimension. The state 1 of WT locates at (-2, -1) of the energy landscape that is 141 not visited by Omicron, and the coordinate (2,1) of the state 1 in the Omicron is also not touched 142 by the WT either. We further extracted the representative structures for each macrostate and 143 analyzed the transition process between each pair of macrostates by transition path theory ( Fig. 3b  144 and Extended Data Fig. 1-2) . In the WT system, state 3 is the dominant one and the transition 145 from states 1 and 2 to it is very quick (< 1 ns). State 2 is the second large macrostate which covers 146 around 8.7 % conformational space and is relatively easy to be reached from states 1 and 3 (both 147 are 4.90 ns). In the Omicron system, the distributions of states 1, 2 and 3 are relatively balanced 148 with a proportion of 16.6%, 31.1%, and 52.3%, respectively, and similar to the WT system, state 3 149 can be quickly reached from states 1 and 2 with 1.07 ns and 0.89 ns, respectively. Furthermore, it 150 can also easily transfer back to the other two states (4.12 ns to state 1 and 1.70 ns to state 2). 151

Considering the differences in the distributions and kinetics of representative states 152 between the WT and omicron systems, we further evaluated the binding abilities for each 153 macrostate by employing MM/GBSA (Fig. 3c) . For the WT simulation system, state 2 is the least 154 suitable conformation for ACE2 binding in all macrostates (ΔG = -3.88 ± 19.64 kcal/mol), state 155 1 has a better binding ability ( Δ G = -26.25 ± 9.76 kcal/mol), and the dominant state 3 is the 156 best-fit conformation for ACE2 binding ( Δ G = -36.20 ± 8.64 kcal/mol). As for the Omicron 157 system, state 1 shows the highest binding affinity with ACE2 (ΔG = -40.84 ± 6.93 kcal/mol), and 158 states 2 and 3 have similar binding affinities (ΔG = -35.54 ± 8.11 and -36.41 ± 10.26 kcal/mol, 159 respectively). Notably, in the WT system, the dominant state 3 has the lowest binding free energy, 160

whereas the other two minor states 1 and 2 have significantly higher binding free energies. In 161 contrast, all three macrostates in the Omicron system have low binding free energies and the 162 differences among them are not as significant as those in the WT system. We also decomposed 163 and compared the binding free energy to residues and picked up the residues with highly different 164 energy values among macrostates in WT system. As shown in Fig. 3c showed that most of the studied antibodies, especially Class I and Class II antibodies, lost their 176 neutralizing activities due to mutations in the epitopes such as K417N and E484A, which is 177 consistent with previous studies [19] [20] [21] . Most of the studied Class III and Class IV antibodies still 178 maintained neutralizing activities. However, a subset of Class III antibodies were reported to be 179 sensitive to the G339D, N440K, and G446S mutations in the Omicron RBD [1, 2] . We paid special 180 attention to the S371L/S373P/S375F mutations within the epitopes of a subset of Class IV 181 antibodies including S2X35, S304, SA24, H104 and CR3022 [22] [23] [24] [25] (Fig. 4a-c) . Indeed, the 182 significantly reduced neutralization activities have been reported for S2X35 and S304 against the 183 Omicron variant [18] . Structural studies have revealed hydrogen-bonding interactions occurred 184 between these two antibodies and the hairpin loop (residue Y369-C379) ( Fig. 4d and 4e) . 16 185 hydrogen bonds are formed at the interface of RBD and these two antibodies. Specially, half of the 186 hydrogen bonds are formed with main chain atoms of RBD residues and the remaining ones are 187 formed between the side chain atoms (Fig. 4d, 4e) . The S371L/S373P/S375F substitutions not 188 only changed the side chains of residues, but also induced a main-chain conformational change, 189 which would disrupt the specific binding of the antibodies to the hairpin loop (residue 190 Y369-C379). It is expected that antibodies including the hairpin loop in their recognizing epitopes Synchrotron Research Facility. Diffraction data were processed using the HKL3000 software [26] 339 and the data-processing statistics are listed in Extended Data Table 1 . 340

The structure was determined using the molecular replacement method with PHASER in the 342 CCP4 suite [27] . The search models used included the ACE2 extracellular domain and 343 SARS-CoV-2 RBD (PDB: 6M0J). Subsequent model building and refinement were performed 344 using COOT [28] and PHENIX [29] , respectively. Final Ramachandran statistics: 96.95% favored, 345 3.05% allowed and 0.00% outliers for the final structure. The structure refinement statistics are 346 listed in Extended Data Table 1 . All structure figures were generated with PyMol [30] . 347

Surface plasmon resonance 348 Binding kinetics of ACE2 and SARS-CoV-2 RBDs were determined by surface plasmon 349 resonance using a Biacore S200 (GE Healthcare). All experiments were performed in a running

HKL-3000: the 428 integration of data reduction and structure solution--from diffraction images to an initial 429 model in minutes

Phaser crystallographic software

Coot: model-building tools for molecular graphics

PHENIX: building new software for automated crystallographic 437 structure determination

PyMod 2.0: improvements in protein 440 sequence-structure analysis and homology modeling within PyMOL

Amino-Acid-Specific Protein Backbone Parameters Trained 443 against Quantum Mechanics Energy Surfaces in Solution

MSMBuilder: Statistical Models for Biomolecular Dynamics

Slow dynamics in protein fluctuations revealed by 448 time-structure based independent component analysis: the case of domain motions. The 449

Describing protein folding kinetics by molecular dynamics 451 simulations. 2. Example applications to alanine dipeptide and beta-hairpin peptide

State Models: From an Art to a Science

Automatic 456 discovery of metastable states for the construction of Markov models of macromolecular 457 conformational dynamics

py: An Efficient Program for End-State Free Energy 460 Calculations

Assessing the performance of MM/PBSA and MM/GBSA methods

Entropy effects on the performance of end-point binding free energy calculation 464 approaches

buffer composed of 10 mM HEPES, pH 7.2, 150 mM NaCl, and 0.005% Tween-20 (v/v). ACE2 351 was immobilized on a CM5 sensor chip (GE Healthcare) to a level of~500 response units. A 352 2-fold dilution series ranging from 50 to 3.125 nM of the SARS-CoV-2 WT RBD and omicron 353 RBD were injected at a flow rate of 30 µl/min (association 60s, dissociation 180s), and the 354 immobilized ACE2 was regenerated using 5mM NaOH for 10s. The resulting data were fit to a 355 1:1 binding model using Biacore Evaluation Software (GE Healthcare). 356 357 The Omicron RBD-ACE2 structure and WT RBD-ACE2 complex (PDB ID: 6M0J) were used to 358 build the MD simulation systems. Besides, K493 in Omicron RBD-ACE2 structure was mutated 359 to R493 to follow the recent sequence of Omicron. Therefore, three simulation systems including 360 WT, Omicron (Q493K) and Omicron(Q493R) were constructed. To keep the consistency among 361 simulation systems, S19-D615 in ACE2 and T333-G526 in RBD were maintained in WT and the 362 two Omicron systems. 363

The FF19SB force field was applied to model the systems [31] . The initial structures were solvated 364 in a truncated octahedral transferable intermolecular potential three point (termed as "TIP3P") 365 water box with a buffer of 10 Å around it. Then, counterions Na+ or Cl-were added to the 366 systems for neutralization and 0.15 mol•L -1 NaCl was added to solvents. 367After construction, the systems were firstly minimized for 15,000 cycles with a restraint of 500 368 kcal•mol −1 •Å −2 on the RBD and ACE2. Then, all atoms encountered 30,000 cycles of minimization. 369Next, the systems were heated from 0 to 300 K in 300 ps and equilibrated for 700 ps with 10 370 kcal•mol −1 •Å −2 positional restraint on non-solvent atoms. Finally, the WT, Omicron (Q493K) and 371Omicron (Q493R) simulation systems encountered 8 parallel rounds of 200 ns production MD 372 simulations, respectively. During simulations, the temperature (300 K) and pressure (1 atm) was 373 controlled by Langevin thermostat and Berendsen barostat, respectively. Long-range electrostatic 374 interactions were treated by the Particle mesh Ewald algorithm with a grid size of 1 Å, and a 375 cutoff of 10 Å was employed for short-range electrostatic and van der Waals interactions. The 376 SHAKE algorithm was applied to restrain the bond with hydrogens. MD simulations were 377 performed on Amber20. 378

Starting with the code base of the current stable version 3.8.0 of MSMBuilder [32] , we developed a 380 more robust algorithm to describe the transition process in Markov state model. Our algorithm 381 modifies the fixed lag time into a random one by a kernel function, which is further used to count 382 transition matrix and build MSM model. As a consequence, the MSM model based on our 383 algorithm exhibits a more robust and powerful representation ability for describing the protein 384 conformational space. We will discuss this method in-depth in our future publications. 385From the trajectories of WT and Omicron (Q493R) system, the time-lagged independent 386 component (tIC) analysis was firstly applied to decrease the dimension of the conformational 387 space [33] . We selected the residue pairs, in which one residue was from RBD and the other was 388 from ACE2. For each residue pair, we measured any pair of distances between the heavy atoms on 389 RBD and those on ACE2 and kept the distances less than 5 Å as inputs for ContactFeaturizer. WT, Omicron (Q493K) and Omicron (Q493R) simulation systems [37] . In total, the binding free 405 energy (ΔGbind) of RBD towards ACE2 is expressed as equation (2). 406Meanwhile, the second law of thermodynamics reflects that ΔGbind equals to enthalpy changes 407 (ΔH) minus the product of entropy changes and temperature (TΔS), as equation (3) expresses. 408The system conformation entropy (termed as "−T∆S") is usually estimated by normal mode 409 analysis with a quasi-harmonic model, however, accurate estimation of the conformation entropy 410 for the protein-protein interactions remains challenging. Notably, the item could be omitted here 411 considering that the differences of enthalpy (termed as "ΔH") are large enough and the similarity 412 of system conformation entropy among simulation systems [38] . Therefore, we omitted the 413 calculation of the −T∆S term and only concentrated on the relative ordering of the free energy 414 changes. 415In the simulation process, ΔH is transformed into the sum of the molecular mechanical energies 416 (ΔEMM) and the solvation free energy (ΔGsolv), according to equation (4). 417In addition, ΔEMM consists of the intramolecular energy (ΔEint including bond, angle and dihedral 418 energies), van der Waals energy (ΔEvdw) and electrostatic energy (ΔEele), while ΔGsolv consists of 419 the polar (ΔGP) and the non-polar items (ΔGnp). Equations (5) and (6) represent them. 420The Generalized Born model was used to calculate ΔGP. ΔGnp was calculated based on the 421 solvent-accessible surface area (SASA) in equation (7). 422The solvation parameter γ and b were 0.00542 kcal (mol −1 ·Å −2 ) and 0.92 kcal/mol, respectively. 423The decomposition of the free energy into residues was subsequently carried out by the 424 MMPBSA.py plugin [37] . During the decomposition process, 1-4 interactions were added to electric 425 interactions and van der Walls interactions. 426

We thank the SSRF BL02U1 and BL10U2 beam line for data collection and processing. We thank 468 the X-ray crystallography platform of the Tsinghua University Technology Center for Protein 469Research for providing the facility support. This work was supported by funds from the National 

We declare no competing interest. 482

The coordinates and structure factors for the Omicron RBD-ACE2 complex were deposited in 484 Protein Data Bank with accession code 7WHH. 485