key: cord-0683979-p550fip1 authors: Socher, Eileen; Heger, Lukas; Paulsen, Friedrich; Zunke, Friederike; Arnold, Philipp title: Molecular dynamics simulations of the delta and omicron SARS-CoV-2 spike – ACE2 complexes reveal distinct changes between both variants date: 2022-02-26 journal: Comput Struct Biotechnol J DOI: 10.1016/j.csbj.2022.02.015 sha: 858328a5a14f09d911107984b363f1f4bab4fbaf doc_id: 683979 cord_uid: p550fip1 SARS-CoV-2, the virus causing the COVID-19 pandemic, changes frequently through the appearance of mutations constantly leading to new variants. However, only few variants evolve as dominating and will be considered as “Variants of Concern” (VOCs) by the world health organization (WHO). At the end of 2020 the alpha (B.1.1.7) variant appeared in the United Kingdom and dominated the pandemic situation until mid of 2021 when it was substituted by the delta variant (B.1.617.2) that first appeared in India as predominant. At the end of 2021, SARS-CoV-2 omicron (B.1.1.529) evolved as the dominating variant. Here, we use in silico modeling and molecular dynamics (MD) simulations of the receptor-binding domain of the viral spike protein and the host cell surface receptor ACE2 to analyze and compare the interaction pattern between the wild type, delta and omicron variants. We identified residue 493 in delta (glutamine) and omicron (arginine) with altered binding properties towards ACE2. The COVID-19 pandemic caused by the novel severe acute respiratory syndrome coronavirus 2 39 (SARS-CoV-2) is driven by newly emerging variants that arise from the original virus. Although most 40 of the genetic changes have little to no impact on the virus' properties, some virus variants arise with 41 higher transmissibility and/or increased virulence. These variants are classified as "Variants of 42 Concern" (VOC) by the WHO. In January 2022, the WHO listed five VOCs labeled with Greek letters 1 : 43 (1) the alpha variant ( 47 All these VOCs listed by the WHO acquired some distinctive mutations ( Fig. 1a lists the mutations of 48 the predominant delta and omicron variants) in the trimeric spike glycoprotein of SARS-CoV-2, which 49 specifically binds with its receptor-binding domain (RBD) to the host cell receptor angiotensin-50 converting enzyme 2 (ACE2) [3] [4] [5] . This binding step is crucial for viral entry into a host cell to initiate 51 infection ( Fig. 1b) 6, 7 . Therefore, understanding the interaction of the virus with ACE2 at the cell 52 surface is fundamentally important in the fight against SARS-CoV-2. The protein structure of the RBD 53 of the wild type spike protein in complex with ACE2 is known (e.g. PDB ID code: 7KMB 8 ) and allows 54 characterization of receptor binding at the atomic level, greatly improving our understanding of host-55 pathogen interactions at this site. By introducing the observed mutations into the wild type RBD-56 ACE2 complex structure with in silico modeling and subsequently conducting molecular dynamics 57 (MD) simulations of these RBD-ACE2 complexes, the influence of specific mutations in the RBD on 58 binding towards the human receptor can be strongly enhanced 9-11 . Experimentally determined 59 protein structures provide the basis for MD simulations, which can add information on the protein 60 dynamics, the flexibility and residue interactions that would be difficult to access with experimental 61 methods. Additionally, MD simulations can today often be performed and analyzed much faster than 62 experimentally determined structures can be obtained and therefore it can serve as a predictive tool. The mentioned alpha, beta, gamma and omicron variants harbor as common feature all a 64 substitution of asparagine-to-tyrosine at position 501 (N501Y) in the RBD, which is not present in the 65 delta variant 12 . In addition to the N501Y mutation, omicron carries 14 other mutations within the 66 RBD (Fig. 1a) . A considerable number of them can be found at the direct binding interface with ACE2 67 (Fig. 1b) . In contrast to omicron, the delta variant has a characteristic leucine-to-arginine substitution carrying this L452R mutation, stronger affinity of the spike protein for the ACE2 receptor was 71 described 13, 14 . The delta and omicron variants share a common T478K mutation and for delta 72 reduced affinity of neutralizing antibodies at this position was shown [15] [16] [17] . In omicron many of the 73 mutations found at the RBD are located within epitopes recognized by neutralizing antibodies. An 74 exemplary array of such neutralizing antibodies is shown in Supplementary Fig. 1 and experimental 75 data suggests reduced binding of antibodies from patient sera 18, 19 . Of note, especially antibodies 76 binding at the RBD at the site of contact formation with ACE2 show high neutralization properties 20 . The present study compares the original wild type (wt), the delta variant and the similar variants 78 epsilon and kappa with the omicron variant. We identify residue 493 within the RBD as a major 79 difference between all three variants. In wt and delta residue 493 is a glutamine (Gln493) Additionally, omicron also carries an amino acid exchange from lysine 417 to asparagine (K417N). 89 The same mutation was already observed in the beta variant, while in the gamma variant lysine 417 90 is mutated to a threonine (K417T). In wild type, Lysine 417 forms a very stable salt bridge with 91 aspartate 30 from the ACE2 receptor and we could show that exchange to asparagine or threonine 92 largely disrupts binding at this position 9 . Here, we show that major changes in interaction with ACE2 occur at positions 417, 493, 501 and at 94 position 505 where a tyrosine is mutated to a histidine (Y505H). We also show, that many of the 95 remaining amino acid exchanges occurring within the RBD do not largely influence binding to ACE2, 96 but rather change epitope areas for neutralizing antibodies (e.g. N440K, G446S, T478K, E484A and 97 G496S) 20 , which might explain the reduced neutralization capacity identified for serum samples from 98 vaccinated patients 18 . 99 Results The contact analyses and electrostatic potential of the viral receptor-binding domain The SARS-CoV-2 delta and omicron variants have both acquired a threonine-to-leucine (T478K) 102 substitution in the RBD. The delta variant carries an additional leucine-to-arginine (L452R) 103 substitution in the RBD and the omicron variant has, in addition to this T478K mutation, 14 other 104 mutations within the RBD (Fig. 1a) . In order to compare individual virus variants, we performed four 105 independent MD simulation runs over 500 ns for every variant's RBD in complex with ACE2. 106 Subsequently, we analyzed the overall intermolecular contacts formed between the human ACE2 107 and the viral RBD of the wild type and the delta and omicron variants. We found a reduced number 108 of contacts between the RBD and ACE2 for the omicron variant, while the number of contacts for the 109 delta variant was not markedly changed (Fig. 1c) . Then, we assigned the individual number of 110 contacts to all RBD amino acids within 8 Å distance to ACE2 and found marked differences for 111 residues 493, 501 and 505 (Fig. 1d, Supplementary Fig. 1a ). Wild type and delta variant express an 112 asparagine at position 501, whereas the omicron variant harbors a tyrosine (N501Y) at this position. 113 The functional consequences of this N501Y mutation were already described 22 and a cryo-electron 114 microscopy structure carrying this N501Y exchange in complex with ACE2 was also experimentally 115 solved (PDB ID code: 7MJN 23 ). In the omicron variant, the tyrosine at position 505 (wt and delta 116 variant) is mutated to histidine and we identified a reduced number of contacts to ACE2 for histidine 117 at this position ( Fig. 1d, Supplementary Fig. 2a and Supplementary Fig. 2b ). In detail, residue 505 118 interacts with the carbon atoms of the side chain of lysine 353 on ACE2 ( Supplementary Fig. 2c ). Here 119 a pi-electron ring system (as for tyrosine expressing variants) might have more contacts than the 120 inserted histidine in omicron. For the epsilon (B.1.427/B.1.429) and kappa (B.1.617.1) variant, we 121 also performed contact analyses and found no marked differences on the individual amino acid level 122 ( Supplementary Fig. 2d) . As the omicron variant shows a considerable number of positively charged 123 amino acids at sites of mutation, we analyzed the electrostatic potential at the RBD-ACE2 interface 124 and found increased electro positivity in the delta variant (L452R and T478K) and even more 125 pronounced in the omicron variant (T478K, Q493R, Q498R and Y505H; Fig. 1e ). However, the 126 omicron variant also loses a positively charged amino acid at the interaction interface at position 417 127 (K417N; Fig. 1b, Fig. 1e ). The interaction between residue 493 and ACE2 is altered in delta and omicron variant As numerous positively charged amino acids appear at the interaction interface, we also analyzed the 141 linear electrostatic interaction energy to identify newly formed or resolved intermolecular salt 142 bridges. Largest differences between wt, delta and omicron were found at positions 417 and 493 143 (Fig. 2a, Supplementary Fig. 3 ). Analyzing the electrostatic potential at the interaction interface on 144 ACE2 revealed aspartate 30, glutamate 35 and aspartate 38 as potential interaction partners for 145 these two residues (Fig. 2b) . For other variants that exchanged lysine 417 to asparagine or threonine, 146 we already identified a strongly reduced electrostatic interaction energy with aspartate 30 9 . Thus it 147 was not very surprising that we found the same reduction for the omicron variant (Fig. 2c ) and a salt 148 bridge in wt (Fig. 2d ) that also remains stable over time (Fig. 2e , Table S1 ). Next, we analyzed contact 149 formation for residue 493 (Fig. 2f ). We found that wt and delta (both with Gln493) behave similar, 150 but that omicron shows some major differences. The delta variant shows a small increase in contact 151 formation with lysine 31 and glutamate 35 (Fig. 2f ). In contrast, omicron changes its interaction 152 profile completely. For one MD simulation run we identified interaction with aspartate 30 (Asp30) 153 and in three runs increased interaction with glutamate 35 (Glu35) and with aspartate 38 (Asp38). 154 Additionally, the interaction with lysine 31 (Lys31) is reduced in all four simulation runs of the 155 omicron variant (Fig. 2f ) probably due to electrostatic repulsion between the positively charged 156 arginine 493 and the also positively charged lysine 31. 169 170 To analyze changes in the delta variant, we pooled it with four individual MD simulation runs of the 172 epsilon and four runs of the kappa variant (all these variants carry the L452R mutation) and 173 compared it to twelve 500 ns long simulation runs of the wt RBD-ACE2 complex (eight of them were 174 already described earlier 9,10 ). Within the neighboring beta sheet of glutamine 493, the L452R 175 exchange can be found ( Fig. 3a, b) , thereby changing the structural environment of glutamine 493. 176 We found that the number of intermolecular contacts between glutamine 493 and the human 177 receptor ACE2 increases in the different L452R variants compared to wild type RBD-ACE2 complexes 178 (Fig. 3c) . Decomposed on the individual amino acid level, although not significant, we identified the 179 ACE2 residues glutamate 35 and to a minor extend lysine 31 with an increased number of contacts to 180 glutamine 493 (Fig. 3d ,e). Contact formation between glutamine 493 and histidine 34 of ACE2 was 181 unchanged, however (Fig. 3d ,e). As glutamine residues engage in hydrogen bonds, we also analyzed 182 the number of intermolecular hydrogen bonds of glutamine 493 with ACE2 and found a significant 183 increase for variants that express the L452R mutation (Fig. 3f) . Detailed analysis showed an increased 184 occurrence of hydrogen bonds with lysine 31 and with glutamate 35 ( The omicron variant forms two new stable salt bridges with ACE2 residues 214 As described above (Fig. 2f) , we identified an entirely altered interaction pattern for arginine at 215 position 493. The number of contacts with glutamate 35 increased and contacts with aspartate 38 216 were newly formed and not present in variants expressing a glutamine at position 493 (Fig. 2f) . For 217 three of our MD simulations, these two newly formed salt bridges can also be appreciated in time 218 course distance plots of arginine 493 with aspartate 30, glutamate 35 and aspartate 38 ( Fig. 4a ; 219 Supplementary Fig. 5a ). However, one of our MD simulation runs showed a different picture and we 220 found a stable intermolecular salt bridge between arginine 493 and aspartate 30 on ACE2 (Fig. 4b) . 221 Analyzing the average distances to favored residues on ACE2 within single runs provided insight into 222 the stability of the individual salt bridges (Fig. 4c) . Especially interactions between arginine 493 and 223 aspartate 30 and glutamate 35 show low standard deviations with average distances below 5 Å, thus 224 strongly supporting ionic interaction. Analysis of the linear electrostatic interaction energy at position 225 493 showed a strong increase for omicron when compared to wt and delta (Fig. 4d) between arginine 498 (only present in omicron) and aspartate 38 (Fig. 4f, Supplementary Fig. 5b ). 232 Together, omicron carrying an arginine at position 493 (Q493R) shows a higher flexibility in contact 233 formation than wt or delta in the same region. As there was a sequence conflict at position 493 at 234 the beginning, we also simulated the RBD-ACE2 complex with a lysine at this position 493. All runs 235 showed stable salt bridges between lysine 493 and glutamate 35 ( Supplementary Fig. 6a ,b) and 236 aspartate 38 ( Supplementary Fig. 6c ,d) of ACE2. A lysine side chain at this position might be too short 237 to interact with aspartate 30 (Supplementary Fig. 6e ). The increase in linear electrostatic interaction 238 energy, however, remains similar ( Supplementary Fig. 6f ). One selection advantage of the omicron 239 over the omicron+Q493K variant might be the more flexible interaction regime at this position. 252 Discussion The appearance of a L452R exchange within the RBD of the SARS-CoV-2 spike protein rendered the 254 affected variants to former "Variants under Investigation" (epsilon, kappa) or current "Variant of 255 Concern" (delta variant) 1 . The delta variant dominated the pandemic worldwide from mid-2021. All 256 three variants carry additional amino acid exchanges within the other parts of the spike protein, 257 including the D614G variant that we identified as a potential molecular switch to liberate the fusion 258 peptide 10 . Here, we show that the L452R mutation induces an increased number of hydrogen bonds 259 formed by glutamine 493 between the RDB (delta variant) and ACE2. Glutamine 493 was also 260 identified as an important interaction residue in the original wt SARS-CoV-2 virus 24-26 . The increased 261 number of hydrogen bonds in delta might induce a stronger interaction between RBD and ACE2 and 262 thus support increased viral infectivity. Comparing the omicron variant to wt and delta reveals major 263 differences in RBD-ACE2 interaction. As in the alpha variant, omicron carries a N501Y mutation and, 264 as the beta variant, carries the K417N mutation. Both residues are important for RBD-ACE2 265 interaction in wt 24,26 and the N501Y mutation apparently increases the binding energy between RBD 266 and ACE2 thus residue 498 is consistently involved in RBD-ACE2 interaction 9,24-26 . As a side note, we could also 276 show that a lysine at this position has a less flexible interaction pattern. Overall, three mutations 277 within the RBD of omicron can be attributed to direct changes in interaction with ACE2 and these are 278 K417N, Q493R and N501Y. With the exchange of tyrosine 505 to histidine, a reduced number of 279 contacts is associated. However, the effect of this exchange cannot be fully judged. When comparing 280 the delta and omicron variant to wt an additional trend becomes evident. Amino acid exchanges 281 within the RBD render it more electropositive. Delta is more positive than wt and omicron more 282 positive than delta. This might enable a better passive adhesion to the negatively charged glycocalyx 283 and here especially to heparin sulfate, which is a critical factor for SARS-CoV-2 binding 27-29 . Thus, we 284 can identify amino acid exchanges which (i) influence direct interaction with ACE2, (ii) change 285 epitopes of neutralizing antibodies and/or (iii) change the electrostatic surface potential. In omicron, 286 five of the fifteen mutations within the RBD insert positively charged amino acids (N440K, T478K, 287 Q493R, Q498R and Y505H) and only one is lost (K417N). All of these six positions are within epitope 288 regions for neutralizing antibodies 20 . In conclusion, omicron started to dominate the pandemic at the 289 end of 2021 with patient numbers rising steeply in affected areas. MD simulation and structural 290 analysis of the RBD-ACE2 complex shows clear differences of the omicron to the delta variant. 291 Omicron has progressed in terms of ACE2 binding, in terms of immune escape and presents a more 292 positively charged interface area. 294 Methods Generation of the starting structures To investigate the interface between the RBD of the spike protein and ACE2, the respective wild type 297 start structure was taken from the PDB database (PDB ID code: 7KMB 8 G339D, S371L, S373P, S375F, 301 K417N, N440K, G446S, S477N, T478K, E484A, Q493R, G496S, Q498R, N501Y, Y505H and for an 302 omicron variant with a lysine at position 493: G339D, S371L, S373P, S375F, K417N, N440K, G446S, 303 S477N, T478K, E484A, Q493K, G496S, Q498R, N501Y, Y505H) were introduced with Swiss-PdbViewer 304 4.1.0 30 (http://www.expasy.org/spdbv/). For the electrostatic analyses of the protein surface, "PQR" output files were generated by using 306 charge parameters assigned with the APBS-PDB2PQR software suite 31 307 (https://server.poissonboltzmann.org/). 309 Molecular dynamics simulations were performed exactly as described before 9,10 . By using version 20 310 of the Amber Molecular Dynamics software package (ambermd.org) 32 , the ff14SB force field 33 and 311 the Amber Tool LEaP, all systems were electrically neutralized with Na + ions and solvated with 312 TIP3P 34 water molecules. The receptor-binding domain complexed with ACE2 was solvated in a water 313 box with the shape of a truncated octahedron and a distance of at least 25 Å from the borders to the 314 solute. Minimization was carried out in three consecutive parts to optimize the geometry of the initial 316 structures. In the first minimization part, all water molecules were minimized, while all other atoms 317 were restrained at the initial positions by using a constant force of 10 kcal·mol -1 ·Å -2 . During the 318 second part, additional relaxation of the sodium ions and the hydrogen atoms of the protein was 319 allowed, while the remaining protein was restrained with 10 kcal·mol -1 ·Å -2 . In the third part, the 320 entire protein, ions, and water molecules were minimized without any restraints. All three 321 minimization parts started with 2500 steps using the steepest descent algorithm, followed by 2500 322 steps of a conjugate gradient minimization. After minimization, the systems were equilibrated in two 323 successive steps. In the first step, the temperature was increased from 10 to 310 K within 0.1 ns and 324 the protein was restrained with a constant force of 5 kcal·mol -1 ·Å -2 . In the second step (0.4 ns 325 length), only the C α atoms of the protein were restrained with a constant force of 5 kcal·mol -1 ·Å -2 . In 326 both equilibration steps, the time step was 2 fs. Minimization and equilibration were carried out on 327 CPUs, while the subsequent production runs were performed using pmemd.CUDA on Nvidia A100 328 GPUs 35-37 . Subsequent 500 ns long production runs were performed without any restraints and at 329 310 K (regulated by a Berendsen thermostat 38 ). Furthermore, the constant pressure periodic 330 boundary conditions with an average pressure of 1 bar and isotropic position scaling were used. For 331 bonds involving hydrogen, the SHAKE algorithm 39 was applied in the equilibration and production 332 phases. To accelerate the production phase of the MD simulations, hydrogen mass repartitioning 333 (HMR) 40 was used in combination with a time step of 4 fs. For statistical analyses, four independent 334 500 ns long MD simulation runs were performed for the epsilon, delta, kappa and omicron spike 335 protein RBD variants in complex with ACE2. For the wild type RBD in complex with ACE216 336 independent 500 ns long MD simulation runs were used (eight completely new MD simulations and 337 in addition also eight MD simulation runs conducted for two earlier studies 9,10 were re-evaluated). -5.599 -4.376 -5.468 -4.906 -6.438 -5.540 -4.076 -3.840 -1.855 -2.118 -11.387 -3 Sequence analysis of the emerging SARS-CoV-2 variant Omicron in 359 South Africa SARS-CoV-2 Cell Entry Depends on ACE2 and TMPRSS2 and Is Blocked by 361 a Clinically Proven Protease Inhibitor Molecular interaction and inhibition of SARS-CoV-2 binding to the ACE2 363 receptor Angiotensin-converting enzyme 2 is a functional receptor for the SARS 365 coronavirus Mechanisms of SARS-CoV-2 entry into cells A pneumonia outbreak associated with a new coronavirus of probable bat 369 origin Cryo-EM Structures of SARS-CoV-2 Spike without and with ACE2 Reveal a pH-371 Dependent Switch to Mediate Endosomal Positioning of Receptor-Binding Domains Computational decomposition reveals reshaping of the SARS-CoV-2-ACE2 374 interface among viral variants expressing the N501Y mutation Mutations in the B.1.1.7 SARS-CoV-2 Spike Protein Reduce Receptor-Binding 377 Affinity and Induce a Flexible Link to the Fusion Peptide Key Interacting Residues between RBD of 379 SARS-CoV-2 and ACE2 Receptor: Combination of Molecular Dynamics Simulation and Density 380 Functional Calculation Impact of the Double Mutants on Spike Protein of SARS-CoV-2 B Lineage on the Human ACE2 Receptor Binding: A Structural Insight Transmission, infectivity, and antibody neutralization of an emerging SARS-385 CoV-2 variant in California carrying a L452R spike protein mutation. medRxiv (2021) Antibody-Mediated Neutralization of Authentic SARS-CoV-2 B.1.617 387 Variants Harboring L452R and T478K/E484Q Neutralising antibody escape of SARS-CoV-2 spike protein: Risk 389 assessment for antibody-based Covid-19 therapeutics and vaccines Long-Time-Step Molecular 444 Dynamics through Hydrogen Mass Repartitioning Software for Processing Molecular Dynamics Trajectory Data py: An Efficient Program for End-State Free Energy 448 Calculations UCSF Chimera-A visualization system for exploratory research and 450 analysis SARS-CoV-2 neutralizing antibody structures inform therapeutic strategies. 452 Studies in humanized mice and convalescent humans yield a SARS-CoV-2 454 antibody cocktail Gln493-OE1 Lys31-HZ3 Gln493-OE1 Lys31-HZ1 Gln493-OE1 Lys31-HZ2 Gln498-OE1 Lys353-HZ1 Lys353-NZ Gln498-OE1 Lys353-HZ2 Lys353-NZ Gln498-OE1 Lys353-HZ3 Lys353-NZ