key: cord-0795900-dhdrveky
authors: Socher, Eileen; Conrad, Marcus; Heger, Lukas; Paulsen, Friedrich; Sticht, Heinrich; Zunke, Friederike; Arnold, Philipp
title: Decomposition of the SARS-CoV-2-ACE2 interface reveals a common trend among emerging viral variants
date: 2021-05-28
journal: bioRxiv
DOI: 10.1101/2021.05.28.446149
sha: 63eaffdc36c0675d782db9460bb24a9b3aa6af47
doc_id: 795900
cord_uid: dhdrveky

New viral variants of the SARS-CoV-2 virus show enhanced infectivity compared to wild type, resulting in an altered pandemic situation in affected areas. These variants are the B.1.1.7 (United Kingdom), B.1.1.7 with the additional E484K mutation, the B.1.351 variant (South Africa) and the P.1 variant (Brazil). Understanding the binding modalities between these viral variants and the host cell receptor ACE2 allows depicting changes, but also common motifs of virus-host cell interaction. The trimeric spike protein expressed at the viral surface contains the receptor-binding domain (RBD) that forms the molecular interface with ACE2. All the above-mentioned variants carry between one and three amino acid exchanges within the interface-forming region of the RBD, thereby altering the binding interface with ACE2. Using molecular dynamics simulations and decomposition of the interaction energies between the RBD and ACE2, we identified phenylalanine 486, glutamine 498, threonine 500 and tyrosine 505 as important interface-forming residues across viral variants. We also suggest a reduced binding energy between RBD and ACE2 in viral variants with higher infectivity, attributed to residue-specific differences in electrostatic interaction energy. Importantly, individual amino acid exchanges not only influence the affected position, but also alter the conformation of surrounding residues and affect their interaction potential as well. We demonstrate how computational methods can help to identify changed as well as common motifs across viral variants. These identified motifs might play a crucial role, in the strategical development of therapeutic interventions against the fast mutating SARS-CoV-2 virus. Significance Statement The COVID-19 pandemic caused by the SARS-CoV-2 virus has significantly changed our lives. To date, there is a lack of neutralizing drugs that specifically target SARS-CoV-2. Hope lies in newly developed vaccines that effectively prevent severe cases of acute respiratory syndrome. However, emerging viral variants escape vaccine-induced immune-protection. Therefore, identification of appropriate molecular targets across viral variants is important for the development of second- and third-generation vaccines and inhibitory antibodies. In this study, we identify residues across viral variants that are important for viral binding to the host cell. As such residues cannot be replaced without diminishing infectivity of the virus, these residues represent primary targets for intervention, for example by neutralizing antibodies.

and inhibitory antibodies. In this study, we identify residues across viral variants that are important for viral binding to the host cell. As such residues cannot be replaced without diminishing infectivity of the virus, these residues represent primary targets for intervention, for example by neutralizing antibodies.

The COVID19 pandemic caused by the SARS-CoV-2 virus is having a major impact on human lives worldwide (1, 2). For cellular infection, the virus engages the cell surface protein angiotensin converting enzyme 2 (ACE2) via its trimeric spike protein (3, 4) . Within the spike protein, the receptor-binding domain (RBD) interfaces with ACE2 (5, 6) (Fig. 1A) , mediating binding of the virus to the host cell surface. Priming at the S2' cleavage site within the spike protein by the serine protease TMPRSS2 releases the fusion peptide from the protein backbone (3, 7, 8) . After insertion of the fusion peptide into the host cell membrane and the activation of a conformational switch leading to dissociation of the spike protein's S1 and S2 domains (9) (10) (11) , the virus and host cell membrane come close and fuse (7, 12) . To allow the two membranes to approach, all non-cleaved subunits of the spike protein must dissociate from the ACE2 receptors, to avoid steric hindrance (5, 13, 14) . Thus, optimized binding efficiency between the RBD and ACE2 is important for viral entry (14) . In several viral variants that showed superior infectivity over the wild type (wt) SARS-CoV-2 virus, amino acid exchanges at the interface of the RBD and ACE2 have been reported (15, 16) . In this study, we examined the RBD variants of B.1.1.7 (originated in the UK), B.1.1.7+E484K (also UK), B.1.351 (South Africa) and P.1 (Brazil) in complex with ACE2 using molecular dynamics (MD) simulations. The B.1.1.7 variant accommodates an N501Y exchange within the RBD-ACE2 interface and the B.1.1.7+E484K exhibits an additional glutamate-tolysine exchange at position 484 (15) . The B.1.351 and P.1 variants carry both of these mutations and have additional exchanges at amino acid position 417: in B.1.351, the lysine at position 417 is exchanged for an asparagine (K417N) and in P.1 for a threonine (K417T) (15) . Notably, the exchange of lysine for glutamate at position 484 is associated with a reduced efficacy for certain vaccines (17, 18) . In this study, we demonstrate that progressive viral variants occur with reduced binding energies to ACE2. Additionally, we show how decomposition of the RBD-ACE2 interface can identify residues on the RBD that are critical for virus host interactions across different SARS-CoV-2 variants. Analyses of electrostatic and van der Waals linear interaction energies reveal the amino acids phenylalanine 486, glutamine 493, and tyrosine 505 as critical residues for binding of the RBD to ACE2 across all RBD variants examined. Our data helps to identify epitopes that are common in all viral variants and might serve as targets for the development of neutralizing antibodies.

To calculate the interaction energy between the cell surface receptor ACE2 and different RBD variants, the structure of the wild type RBD and ACE2 complex (PDB ID: 7kmb (19) ) was used as starting point. Variant-specific amino acid substitutions were introduced in the RBD sequence before applying MD simulations: B.1.1.7 (N501Y), B.1.1.7+E484K (N501Y, E484K), B.1.351 (N501Y, E484K, K417N) and P.1 (N501Y, E484K, K417T). After simulation for 500 ns, we calculated the overall linear interaction energy between the different RBDs and the host cell receptor ACE2 (Fig. 1B) . We found the strongest interaction of ACE2 with the wild type, which was significantly stronger compared to variants B.1.1.7, B.1.351 and P.1. We also calculated the binding energy using molecular mechanics/Generalized Born surface area (MM/GBSA) and found a similar trend (Fig. S1A ). When we plotted the sum of the electrostatic and the van der Waals parts of the linear interaction energy for all RBD residues within 8 Å to ACE2 (Fig. S1B) , we identified eight residues that we analyzed in more detail (Fig. 1C) : amino acid positions 417 (lysine, asparagine, or threonine), 484 (glutamate or lysine), 486 (phenylalanine), 493 (glutamine), 500 (threonine), 501 (asparagine or tyrosine), and 505 (tyrosine). The distribution of the overall interaction energy from the RBD in complex with ACE2 revealed the central part of the RBD as the main contact site (Fig. 1D for wild type; Fig. S1C Illustrating the interaction interface of the wild type RBD with all residues as spheres with radii and colors according to their overall linear interaction energy (yellow (5 kcal/mol)blue (-50 kcal/mol)), demonstrated that lysine 417 and glutamine 493 are major interaction partners with ACE2 (Fig. 1E ). Both amino acids are flanked by phenylalanine 486 on one side and glutamine 498, threonine 500, asparagine 501, and tyrosine 505 on the other side (Fig. 1E) . Depiction of the mutant RBD binding interfaces revealed significant structural changes within the ACE2 interaction site. For example, B.1.1.7 showed a major change in the RBD as glutamine 498 loses its interaction energy compared with the wild type variant (Fig. 1E, F) . In the B.1.1.7+E484K variant, the newly inserted lysine at position 484 introduces interaction energy in an area close to phenylalanine 486 (Fig. 1G) . Exchange of lysine 417 to asparagine (B.1.351, Fig. 1H ) or threonine (P.1, Fig. 1I ) reduced the importance of amino acid 417 in ACE2 binding. This resulted in glutamine 493 as the most important interacting residue flanked by two balanced interaction sites formed by phenylalanine 486 and lysine 484 on one side and threonine 500, tyrosine 501, and tyrosine 505 on the other side (Fig. 1H, I) . We also evaluated the root mean square fluctuation (RMSF) values for all RBD variants in complex with ACE2 to assess the overall structural flexibility of the different RBD variants, but we found no significant differences in structural dynamics (Fig. S1D) .

Analyses of the electrostatic interaction energies between the RBD and ACE2 showed reduced energies for all mutated variants, with significant differences for the B.1.1.7, B.1.351, and P.1 variants compared to wild type ( Fig. 2A) . Decomposition of the electrostatic linear interaction energy to the level of individual amino acids (Fig. S2A , B) revealed residues 417, 484, 493, and 498 as strongest interaction sites for different spike variants (Fig. 2B ). In contrast, the van der Waals linear interaction energy was similar between the different RBD variants analyzed (Fig. 2C) , and the decomposition revealed only minor differences for individual residues (Fig. 2D, Fig. S2C , D). To analyze the consequences of individual residues on ACE2 binding at the molecular level, we measured interatomic distances, calculated the number of contacts and linear interaction energies for specific amino acids. Interestingly, we discovered that the lysine at position 417, found in wild type and RBD variants B.1.1.7 and B.1.1.7+E484K, formed a salt bridge with aspartate 30 of the ACE2 receptor (Fig. 2E) . Accordingly, the distance plot of the two amino acids lysine 417 (RBD wt) and aspartate 30 (ACE2) showed that both residues were in close proximity (<4 Å) underlining their tight interaction as characteristic for salt bridges (Fig. 2F ). This indicates the importance of the amino acid at position 417 in RBD-ACE2 interaction. In the following, we therefore analyzed the effect of mutations at this position as reported for variants B.1.351 and P.1. Calculation of the average numbers of contacts per frame showed that a lysine at position 417 (wt, B.1.1.7, B.1.1.7+E484K) induces ~35 contacts, an asparagine (P.1.351) ~5 contacts and a threonine (P.1) ~2 contacts (Fig. 2G ). This is reflected by a strongly reduced electrostatic interaction energy with ACE2 for this residue when comparing RBD variants carrying a lysine (wt, B.1.1.7, B.1.1.7+E484K) or lacking a lysine residue (Asn/Thr; P.1.351, P.1) at position 417 (Fig. 2h ). An additional effect that we noticed was the increase in electrostatic interaction energy for glutamine 493 in cases where lysine was absent at position 417, as found for variants B.1.351 and P.1 (Fig. 2I ). Since all variants lacking a lysine at position 417 also exhibit an exchange of glutamate to lysine at position 484 (E484K), we assume a compensatory effect counteracting the reduced binding energy.

As indicated by mutation and interaction energy data, residue 484 plays a critical role in the RBD-ACE2 binding properties. The exchange from glutamate to lysine at position 484 (as found in B.1.1.7+E484K, B.1.351 and P.1) changes a negatively charged residue to a positive one, allowing salt bridge formation with negatively charged residues in adjacent ACE2 regions. In this region, we identified glutamate 75 expressed on ACE2 as a potential interaction partner (Fig. 3A) . Plotting the total fraction of contacts for residue 484 (as glutamate or lysine) with residues expressed on ACE2 revealed increased interaction with glutamate 75 for variants expressing a lysine at position 484 (B.1.1.7+E484K, B.1.351 and P.1) (Fig. 3B ). All variants showed interaction with leucine 79. Comparison of variants expressing a glutamate (wt and B.1.1.7) or a lysine (B.1.1.7+E484K, B.1.351 and P.1) at position 484 also showed an increase in the total number of contacts per frame (Fig. 3C) . Comparison of the electrostatic linear interaction energies for variants expressing glutamate or lysine at position 484 demonstrated a significantly increased interaction energy for lysine-expressing variants. Thus, our analysis supports the idea that the exchange of a lysine at position 484 for a glutamate increases the interaction energy and thereby stabilizes the RBD-ACE2 interaction. It also introduces electrostatic interaction in close proximity to the important van der Waals interaction formed by phenylalanine 486 (Fig. 1G -I).

Variants expressing a tyrosine at position 501 instead of an asparagine (all except for wild type) showed a significant decrease in the electrostatic interaction energy between the residue at position 501 and ACE2 (Fig. 4A ). Although the van der Waals interaction energy increased for these variants compared to wild type (Fig. 4B) , the resulting overall interaction energy for this residue decreased when asparagine was exchanged for tyrosine (Fig. 4C) . Analysis of the interaction energy of residues in the vicinity of residue 501 revealed that glutamine 498 lost almost all of its interaction energy (from about -22 kcal/mol to -4 kcal/mol) when a tyrosine was present at position 501 (Fig. 1C ). This loss of overall interaction energy can be attributed to a loss of electrostatic interaction energy. Comparative analyses of variants with identical residues at this position (asparagine for wild type and tyrosine for all other variants), revealed a reduced electrostatic interaction energy from about -18 kcal/mol to -2 kcal/mol (Fig. 4D) . Analysis of the RBD-ACE2 interface around residue 501, which includes glutamine 498, showed a conformational change for this glutamine induced by the insertion of the tyrosine at position 501 (Fig.  4E, F) . Analysis of residues interacting with glutamine 498 on ACE2 revealed major interactions with aspartate 38, tyrosine 41, glutamine 42, leucine 45, and lysine 353 (Fig. 4E, G) . In variants expressing a tyrosine at position 501, glutamine lost its contacts completely to aspartate 38 and lysine 353, and partially to tyrosine 41. This can be explained by the bulky tyrosine inserted for an asparagine displacing glutamine 498 from the binding interface. With these results, we demonstrate that single amino acid changes also affect surrounding residues and alter the conformation in their immediate vicinity.

Understanding the interaction between SARS-CoV-2 and the host cell is important to develop therapeutic antibodies that work across different viral variants. Recently emerged viral variants, such as B.1.1.7, B.1.1.7+E484K, B.1.351, and P.1 display mutations within the RBD at the interface with ACE2 (15, 20) . Here we demonstrate how MD simulations allows us to track these changes in interaction energies between RBD and ACE2. A general trend toward reduced interaction energies between the RBD and ACE2 indicated that not only binding, but also subsequent dissociation of the virus from the cell surface receptor serves as a driving force in viral evolution (21, 22) . Decomposition of the interaction energies into electrostatic and van der Waals interactions demonstrates that the decrease in electrostatic interaction is induced primarily by the exchange of asparagine 501 for tyrosine and the exchange of lysine 417 to asparagine or threonine. For both mutations at position 417, reduced binding affinity to ACE2 but increased expression of the RBD was measured in a yeast system (23) . The E484K mutation induces a local increase in electrostatic interaction energy adjacent to phenylalanine 486. Comparing the spatial distribution of the overall interaction energies for the different variants, it appears that the insertion of N501Y, K417N/T, and E484K as present in B.1.351 and P.1 equalizes the interaction energy uniformly across over the ACE2 binding interface. For the N501Y variants, we can also describe local conformational changes that partially exclude glutamine 498 from its interaction partners on ACE2. These local conformational changes may also play a role in reduced binding of therapeutic antibodies or inadequate protection after vaccination (18, 22, 24) .

Our decomposition analysis also allows us to point to residues that seem to have a critical function for RBD-ACE2 interaction across different spike protein variants. These residues are phenylalanine 486, glutamine 493, threonine 500, asparagine/tyrosine501, and tyrosine 505. Further combinatorial mutation experiments on these residues will determine whether they are dispensable for ACE2 binding. From our data, we hypothesize that a loss of glutamine 493 may prevent sufficient binding and thus viral entry. As no spike variant with higher infectious potential has yet been reported to replace glutamine 493, this residue might be a primary target for a neutralizing antibody acting across all SARS-CoV-2 variants. Our data and data from other groups indicate that MD simulation is a powerful tool to evaluate the consequences of individual amino acid exchanges on the interaction energies between the RBD and ACE2 and beyond (25, 26) . In particular, for variants with multiple changes or rapidly mutating strains, MD simulations serve as a rapid and powerful tool to understand the impact of each mutation and unravel its mechanistic properties. Validation of these results in biochemical or cell biological experiments and circling back into the simulation enables rapid identification of suitable target sites that can be used for translational approaches.

To investigate the interface between the RBD of the spike protein and ACE2, the respective wild type start structure was taken from the PDB database (PDB ID code: 7KMB (19) ). To also generate the starting structures for the MD simulations of the different spike variants, the amino acid substitutions 

Molecular dynamics simulations were performed using version 20 of the Amber Molecular Dynamics software package (ambermd.org) (27) and the ff14SB force field (28) . Using the Amber Tool LEaP, all systems were electrically neutralized with Na+ ions and solvated with TIP3P (29) water molecules. The receptor-binding domain complexed with ACE2 was solvated in a water box with the shape of a truncated octahedron and a distance of at least 25 Å from the borders to the solute. The simulations followed a previously applied protocol (30) . First, minimization was carried out in three consecutive steps to optimize the geometry of the initial structures. In the first step of the minimization, the water molecules were minimized, while all other atoms were restrained at the initial positions with a constant force of 10 kcal·mol -1 ·Å -2 . In the second step, additional relaxation of the sodium ions and the hydrogen atoms of the protein was allowed, while the remaining protein was restrained with 10 kcal·mol -1 ·Å -2 . In the last step, no restraints were used, so the entire protein, ions, and water molecules were minimized. All three minimization parts started with 2500 steps using the steepest descent algorithm, followed by 2500 steps of a conjugate gradient minimization. After minimization, the systems were equilibrated in two successive steps. In the first step, the temperature was increased from 10 to 310 K within 0.1 ns and the protein was restrained with a constant force of 5 kcal·mol -1 ·Å -2 . In the second step (0.4 ns length), only the Cα atoms of the protein were restrained with a constant force of 5 kcal·mol -1 ·Å -2 . In both equilibration steps, the time step was 2 fs. Minimization and equilibration were carried out on CPUs, while the subsequent production runs were performed using pmemd.CUDA on Nvidia A100 GPUs (31) (32) (33) . Subsequent 500 ns long production runs were performed without any restraints and at 310 K (regulated by a Berendsen thermostat (34)). Furthermore, the constant pressure periodic boundary conditions with an average pressure of 1 bar and isotropic position scaling were used. For bonds involving hydrogen, the SHAKE algorithm (35) was applied in the equilibration and production phases. To accelerate the production phase of the molecular dynamics (MD) simulations, hydrogen mass repartitioning (HMR) (36) was used in combination with a time step of 4 fs. For all different spike protein variants in complex with ACE2, the MD simulations were performed four times. Trajectory analysis (analysis of root-mean-square fluctuations (RMSF), analysis of contacts (always with distance criterion of ≤5 Å between any pair of atoms; total fraction of contacts for residue pairs), measurement of interatomic distances, calculation of linear interaction energy (electrostatic and van der Waals interactions) was performed using the Amber tool cpptraj (37). For the calculation of interaction energies between ligand and receptor based on the MM/GBSA method, the mm_pbsa.pl script with default parameters was used (38-40).

Statistical analyses were performed with GraphPad Prism (version 8.0.0 for Windows, GraphPad Software, San Diego, California USA, www.graphpad.com) and statistical tests were applied as indicated below the figure. Plots were created in GraphPad and Gnuplot (version 5.2). All structure images were made with UCSF Chimera 1.15 (41). (19)). (B) Overall linear interaction energy as calculated for the different spike protein variants (wt = wild type). Statistical analysis was performed using one-way ANOVA (n=4, differences assumed significant for ** p<0.01, *** p<0.001). (C) Heat map of receptor-binding domain (RBD) expressed residues important for binding of ACE2. (D) Wild type RBD in complex with ACE2 (aquamarine). All residues within a maximum distance of 8 Å to ACE2 are displayed according to their interaction energy with different colors and sphere radii. The interacting residues within a maximum distance of 8 Å are shown according to their linear interaction spheres of different diameter according to their electrostatic linear interaction energy with ACE2. The residues interacting with ACE2 are also shown without ACE2 as a separate view (right). (C) Heat map for van der Waals linear interaction energies for all residues of the receptor-binding domain (RBD) that were in a maximum distance of 8 Å with ACE2. Colors are according to the van der Waals linear interaction energy from strong (blue) to weak (yellow). (D) The different RBD variants are shown in complex with ACE2 (left) and residues are shown in different colors (blue-yellow) and spheres of different diameters according to their van der Waals linear interaction energy with ACE2. The residues interacting with ACE2 are also shown without ACE2 as a separate view (right). 

Epidemiological and clinical characteristics of 99 cases of 2019 novel coronavirus pneumonia in Wuhan, China: a descriptive study

World Health Organization, WHO Coronavirus (COVID-19) Dashboard

SARS-CoV-2 Cell Entry Depends on ACE2 and TMPRSS2 and Is Blocked by a Clinically Proven Protease Inhibitor

A pneumonia outbreak associated with a new coronavirus of probable bat origin

Molecular interaction and inhibition of SARS-CoV-2 binding to the ACE2 receptor

Biomechanical characterization of SARS-CoV-2 spike RBD and human ACE2 protein-protein interaction

A Multibasic Cleavage Site in the Spike Protein of SARS-CoV-2 Is Essential for Infection of Human Lung Cells

Proteolytic activation of the SARScoronavirus spike protein: cutting enzymes at the cutting edge of antiviral research

Activation of the SARS coronavirus spike protein via sequential proteolytic cleavage at two distinct sites

Engineering human ACE2 to optimize binding to the spike protein of SARS coronavirus 2

Host cell entry of Middle East respiratory syndrome coronavirus after two-step, furin-mediated activation of the spike protein

Protease-mediated enhancement of severe acute respiratory syndrome coronavirus infection

Structural investigation of ACE2 dependent disassembly of the trimeric SARS-CoV-2 Spike glycoprotein

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

Centers for Disease Control and Prevention, Science Brief: Emerging SARS-CoV-2 Variants

The E484K mutation in the SARS-CoV-2 spike protein reduces but does not abolish neutralizing activity of human convalescent and post-vaccination sera. medRxiv : the preprint server for health sciences

Escape from neutralizing antibodies by SARS-CoV-2 spike protein variants

Cryo-EM Structures of SARS-CoV-2 Spike without and with ACE2 Reveal a pH-Dependent Switch to Mediate Endosomal Positioning of Receptor-Binding Domains

SARS-CoV-2 Variants of Concern in the United States-Challenges and Opportunities

Evolutionary and structural analysis elucidates mutations on SARS-CoV2 spike protein with altered human ACE2 binding affinity

Impact of South African 501.V2 Variant on SARS-CoV-2 Spike Infectivity and Neutralization: A Structure-based Computational Assessment

Deep Mutational Scanning of SARS-CoV-2 Receptor Binding Domain Reveals Constraints on Folding and ACE2 Binding

SARS-CoV-2 501Y.V2 escapes neutralization by South African COVID-19 donor plasma. bioRxiv : the preprint server for biology

Enhanced Binding of SARS-CoV-2 Spike Protein to Receptor by Distal Polybasic Cleavage Sites

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

Simmerling, ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99SB

Comparison of simple potential functions for simulating liquid water

The conformational stability of nonfibrillar amyloid-β peptide oligomers critically depends on the C-terminal peptide length

Routine Microsecond Molecular Dynamics Simulations with AMBER on GPUs. 2. Explicit Solvent Particle Mesh Ewald

SPFP: Speed without compromise-A mixed precision model for GPU accelerated molecular dynamics simulations

Molecular dynamics with coupling to an external bath

Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes

Long-Time-Step Molecular Dynamics through Hydrogen Mass Repartitioning

The bulky tyrosine 501 seemed to partially exclude glutamine 498 from the RBD-ACE2 interface. (G) Total fraction of contacts for glutamine 498 with residues expressed on ACE2. Contacts were mainly lost towards aspartate 38 (Asp38), tyrosine 41 (Tyr41) and lysine 353 (Lys353)

energy in different colors and sphere diameters. (D) Root mean square fluctuation (RMSF) values for the different RBD variants show only subtle changes across different spike protein variants

Supplemental Figure 2. Decomposition into electrostatic and van der Waals linear interaction energies. (A) Heat map for electrostatic linear interaction energies for all residues of the receptorbinding domain (RBD) that were in a maximum distance of 8 Å with ACE2. Colors are according to the electrostatic linear interaction energy from strong (blue) to weak (yellow). (B) The different RBD variants are shown in complex with ACE2 (left) and residues are shown in different colors

The authors gratefully acknowledge the compute resources and support provided by the Erlangen Regional Computing Center (RRZE) and by NHR@FAU.