Insights into SARS-CoV-2's Mutations for Evading Human Antibodies: Sacrifice and Survival.

The highly infectious SARS-CoV-2 variant B.1.351 that first emerged in South Africa with triple mutations (N501Y, K417N, and E484K) is globally worrisome. It is known that N501Y and E484K can enhance binding between the coronavirus receptor domain (RBD) and human ACE2. However, the K417N mutation appears to be unfavorable as it removes one interfacial salt bridge. Here, we show that despite the decrease in binding affinity (1.48 kcal/mol) between RBD and ACE2, the K417N mutation abolishes a buried interfacial salt bridge between the RBD and neutralizing antibody CB6. This substantially reduces their binding energy by 9.59 kcal/mol, thus facilitating the process by which the variant efficiently eludes CB6 (including many other antibodies). Our theoretical predictions agree with existing experimental findings. Harnessing the revealed molecular mechanisms makes it possible to redesign therapeutic antibodies, thus making them more efficacious.

Introduction

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2)[1] is a positive-sense single-stranded RNA virus that causes the ongoing coronavirus disease 2019 (COVID-19) pandemic. Because SARS-CoV-2 is a highly transmissible and pathogenic novel coronavirus, the global morbidity and mortality rates of COVID-19 have been quite substantial since it was first detected at the end of 2019.[2] This has brought together researchers from the international scientific community to combat COVID-19 with extraordinary effort. Although we have come a long way in a short period of time, it is still challenging to have the pandemic under full control especially with the recent discovery of several new variants of SARS-CoV-2, such as B.1.1.7 (or 501Y.V1) from the United Kingdom, B.1.351 (or 501Y.V2) from South Africa, and P.1 from Brazil. With the increased transmissibility and hence rapid spreading of these newly identified lineages, serious concerns over whether they could undermine the currently available vaccines or natural immunity of people previously infected with COVID-19 have been raised.[3] Unfortunately, a recent experiment provided evidence that variant B.1.351 eludes natural and vaccine-induced sera.[4]

Experimental studies have demonstrated that angiotensin converting enzyme 2 (ACE2), a protein expressed on the surface of human cells in various organs, plays a crucial role in the viral infection of SARS-CoV-2.[5−7] As a prelude to SARS-CoV-2’s cell entry, the receptor binding domain (RBD) of the spike glycoprotein (S-protein)[8] on the virion surface binds ACE2 on a host cell. Therefore, the S-protein of SARS-CoV-2 becomes a primary target of the host immune system and is used as the leading antigen for vaccine development. Despite being detected from three different continents, the fact that the three aforementioned SARS-CoV-2 variants share some of the same mutations at the S-protein’s RBD suggests that these mutations might have conferred an evolutionary advantage to the virus. While various mutations as well as deletions outside the RBD of S-protein are also important for the enhanced fitness (of variants) for entering host cells, here we focus on mutations in the RBD that harbors the binding site of ACE2 and the epitopes for neutralizing monoclonal antibodies (mAbs).

For the U.K. variant (with the N501Y mutation on the RBD), experimental studies have demonstrated that the reproductive number that measures its infectiousness is approximately 0.4–0.7 higher than those of other strains of the virus[9] and determined recently it is unlikely to escape BNT162b2 vaccine-mediated protection.[10] Our previous work[11] showed that this mutation in the U.K. variant can increase the RBD’s binding affinity for ACE2 but has no obvious effect on mAbs. In addition to the N501Y mutation, the South Africa and Brazil variants also contain the K417N and E484K mutations. Recent studies showed that the E484K mutation not only can enhance RBD–ACE2 binding[12] but also can help the virus escape the therapeutically relevant mAbs.[13] However, reflected in the paucity of existing work, the significance of the K417N mutation so far is still elusive, which prevents us from fully understanding the infection mechanism of these new variants. As shown in Figure a, in the wild-type RBD, K417 forms a salt bridge with D30 in ACE2 (human). Consistently, animals such as pangolin, sheep, Siberian tiger, horseshoe bat, hamster, cat, horse, cow, and ferret having either D30 or E30 (after sequence alignment) in their ACE2 are susceptible to SARS-CoV-2, while other animals such as chicken (A30), rat (N30), and duck (A30) are not.[14] This suggests that the salt bridge formed by K417 (in the RBD) and D30/E30 (in ACE2) is important for RBD–ACE2 recognition. Thus, the K417N mutation resulting in the abolishment of this favorable interfacial interaction (i.e., reducing the RBD–ACE2 binding affinity as verified in experiment[12]) is highly unintuitive. Here, we are motivated to investigate the molecular mechanism of the K417N mutation, to unveil its key benefits for the virus to evolve through this path.

Figure 1

MD simulation systems. (a) Illustration of the salt bridge formed by D30 in ACE2 and K417 in the RBD at the interface. (b) Simulation setup for the RBD–ACE2 complex. (c) Simulation setup for the RBD–CB6 complex. Proteins are shown in cartoon representation, with RBD colored gray, ACE2 colored blue, and the heavy (VH and CH) and light (VL and CL) chains of the CB6’s Fab colored purple and green, respectively. Na+ and Cl– are shown as yellow and cyan spheres, respectively. Water is shown transparently.

Complementary to experimental efforts, all-atom molecular dynamics (MD) simulations with sophisticated and well-calibrated force fields have been widely used to image nanoscale events and investigate the molecular mechanism of proteins.[15−18] In this work, we conducted a computational analysis of the K417N mutation in the South Africa variant using MD simulations with explicit solvent, aiming to gain a better understanding of its underlying molecular mechanism. In addition to the RBD–ACE2 interaction, we also investigated the RBD’s interaction with mAbs. In particular, we explored the mAb CB6 that recognizes an epitope site in the RBD overlapping the binding site of ACE2 and investigated its binding competition with ACE2 over RBD. Our results may help provide invaluable insights into why K417N has been selected during viral evolution and inspire a better design of more efficacious mAbs for treating COVID-19 patients infected with the new SARS-CoV-2 variants.

Results

Figure b illustrates the simulation system for modeling the interaction between ACE2 and RBD (see the Experimental Section for detailed simulation protocols). Briefly, taken from the crystallographic structure (PDB entry 6M0J), the structure of the RBD–ACE2 complex was solvated in a 0.15 M NaCl electrolyte. Similarly, as shown in Figure c, we built the simulation system for the complex of RBD and one Fab of human antibody CB6 with their atomic coordinates taken from the crystallographic structure (PDB entry 7C01). Hereafter, we simply refer to the Fab as CB6.

During the ∼350 ns MD simulation, the RBD–ACE2 complex starting from the structure in the crystal environment was equilibrated in the physiology-like environment (a 0.15 M electrolyte). Figure S1 shows the root-mean-square deviation (RMSD) of protein-backbone atoms in the RBD bound with either ACE2 or CB6. After ∼40 ns, the RMSD values saturated at around 1.7 Å for the RBD bound on CB6 and 1.5 Å for the RBD bound on ACE2. These small RMSD values suggest that the secondary structure of the modeled RBD (with four internal disulfide bonds) was very stable and properly equilibrated. We also calculated the interfacial contact areas for the complexes using the solvent accessible surface area method.[19] On average, the contact area is ∼8.6 nm2 between the RBD and ACE2 while that for the RBD–CB6 contact is larger and is ∼10.3 nm2, as shown in Figure a. During the entire simulation time, values of contact areas fluctuated around the mean, indicating that these two complex structures were stable and equilibrated in the electrolyte.

Figure 2

MD simulations of both RBD–ACE2 and RBD–CB6 complexes. (a) Time-dependent contact area S in simulated RBD–ACE2 (blue) and RBD–CB6 (orange) complexes. (b) Time-dependent numbers of heavy (or non-hydrogen) atoms in ACE2 (blue) or CB6 (orange) that are within 5 Å of K417 in the RBD. (c) Illustration of the salt bridge formed by K417 in the RBD and D30 in ACE2 that is present on the complex’s surface (i.e., being exposed to water). (d) Illustration of the salt bridge formed by K417 in the RBD and D104 in CB6 (VH) that is buried inside the complex. The RBD, ACE2, VH, and VL of CB6 are shown in molecular surface representation. VL of CB6 is transparent to allow viewing of the buried salt bridge.

By examining the simulation trajectories, we found that K417 indeed played an important role in stabilizing the RBD–ACE2 and RBD–CB6 complexes. In Figure b, we highlight the time-dependent number N of heavy (or non-hydrogen) atoms in ACE2 or CB6 that were within 5 Å of K417 in the RBD. Notably, K417 in RBD interacts with many more atoms in CB6 than in ACE2. The average saturated number N for the RBD–CB6 complex is ∼23 atoms (orange line in Figure b), while that for the RBD–ACE2 complex is only ∼6 atoms (blue line in Figure b). Additionally, the number N for the RBD–ACE2 complex fluctuates much more than that for the RBD–CB6 complex, which indicates that the interaction in the latter complex is more stable. Due to the K417–D104 salt bridge, K417 in the RBD stably contacted nearby residues (including D104) in the heavy chain. However, the distance between K417 in the RBD and residues (such as T94 and P95) in the light chain increased slightly to accommodate water molecules that entered into the hydrophilic interfacial area, resulting in a decrease in the contact number N at the beginning of the MD simulation (see the orange line in Figure b).

Among those residues (in ACE2 or CB6) in contact with K417 in the RBD, there are two key salt bridges. One is formed by K417 in the RBD and D30 in ACE2[18] (as shown in Figure a). We illustrate in Figure c that this salt bridge is on the surface of the RBD–ACE2 complex. Thus, this salt bridge is exposed to water and constantly disrupted by polar water molecules, accounting for the observed fluctuations in N (blue line in Figure b). The other one is formed by K417 in the RBD and D104 in CB6[20] as shown in Figure d. Remarkably, this salt bridge is buried among the heterotrimer composed of the RBD and two variable domains of heavy chain VH and light chain VL of CB6. This salt bridge buried inside the protein complex was very stable during the simulation, consistent with the nearly constant number N in Figure b (orange line).

To further demonstrate the stability of these two salt bridges, we define distance D between the NZ atom in the lysine (K) and the CG atom in the aspartate (D), as shown in the inset of Figure a. The histograms in Figure a show the probability distributions of distance D. For the salt bridge in the RBD–CB6 complex, the most probable distance in the single-peak distribution is ∼3.2 Å, corroborating the idea that this buried salt bridge is very stable. For the salt bridge in the RBD–ACE2 complex, the most probable distance is in the first peak and is ∼3.5 Å. This larger distance (as compared with that for the buried salt bridge) reflects the weakened electrostatic interaction inside the water-exposed salt bridge. Furthermore, an additional peak at larger distances (∼5.5 Å) is an indication of a state of the broken salt bridge. Indeed, we observed in the simulation trajectory that this salt bridge is susceptible to polar water molecules and can break and re-form from time to time.

Figure 3

Characterizing salt bridges at the RBD–CB6 and RBD–ACE2 interfaces. (a) Distribution of distance D between a pair of residues in a salt bridge, at the RBD–CB6 interface (orange) and at the RBD–ACE2 interface (blue). The inset shows the definition of distance D. (b–e) Thermal dynamic cycle to obtain the binding free energy change induced by the K417N mutation. The mutations in the bound and free states are shown in panels b and c and panels d and e, respectively.

Therefore, in the wild-type RBD, K417 interacts more strongly with D104 in CB6 than with D30 in ACE2. The calculations of the binding free energy change induced by the K417N mutation were accomplished using the free energy perturbation (FEP) method.[21] As required in FEP calculations, we performed 177 ns MD simulations of the RBD in a 0.15 M NaCl electrolyte (a free state), as shown in Figure S2. With protein structures for both bound (or complex) and free states in respective MD simulations, we employed the FEP alchemy method to calculate the binding free energy difference for the K417N mutation on the RBD, using the thermodynamic cycle shown in panels b–e of Figure . By definition, binding free energy changes for the RBD with either ACE2 or CB6 (due to the K417N mutation) can be obtained as ΔΔG = ΔG2 – ΔG1 (see Figure ). In practice, it is not easy to directly calculate ΔG1 and ΔG2, which can be circumvented by computing ΔGA and ΔGB instead using the thermodynamic cycle (see Figure ). Therefore, ΔΔG = ΔGA – ΔGB. Through the ensemble average,[21]ΔGA and ΔGB can be calculated theoretically (see the Experimental Section) as , where kB is the Boltzmann constant, T is the temperature, and Hi and Hf are the Hamiltonians for the initial (i) stage with K417 (see Figure b,d) and the final (f) stage with N417 (see Figure c,e), respectively.

The results from the FEP calculations are summarized in Table . In the free state (Figure d,e), the K417N mutation yielded a free energy change ΔGB of −32.08 kcal/mol. In the bound state for the RBD–CB6 complex (Figure b,c), free energy change ΔGA = −22.49 kcal/mol. Altogether, the obtained value of ΔΔG is 9.59 kcal/mol, suggesting that the K417N mutation significantly reduced the binding affinity between the RBD and CB6. Similarly, for the RBD–ACE2 complex, we have a ΔGA of −30.60 kcal/mol, and consequently, ΔΔG is 1.48 kcal/mol. Thus, the K417N mutation also reduced the binding affinity between the RBD and ACE2, but the reduction is ∼6.5 times less than that between the RBD and CB6. Previously, it was demonstrated in an experiment that the K417N mutation weakened the binding affinity between the RBD and ACE2,[12] which is consistent with our simulation results. The noticeably small errors listed in Table manifest the accuracy and convergence afforded by the FEP methodology.

Table 1

Values of ΔΔG for the K417N Mutation for the Binding of the RBD to ACE2 and CB6 (mAb)

	ΔGA,B (kcal/mol)	ΔΔG (kcal/mol)
RBD only	–32.08 ± 0.24	–
RBD–ACE2	–30.60 ± 0.33	1.48 ± 0.41
RBD–CB6	–22.49 ± 0.20	9.59 ± 0.31

The free energy loss due to the K417N mutation is mainly caused by the change in electrostatic energy. Dielectric constant εp inside a protein is approximately 4,[22] while dielectric constant εw of water is ∼90 for the TIP3P model.[23] If we approximate the protein–water interface with a planar surface separating two dielectric materials (εp and εw), the effective dielectric constant on the interface can be calculated as (εp + εw)/2, from a simple continuum electrostatics calculation. Thus, the ratio γ of electrostatic energy changes for salt bridges on the protein surface and in water can be roughly estimated as (εp + εw)/2εw ∼ 11.8 (larger than our simulation result γ = 6.5). However, the more realistic explicit model with atomic details[24] predicts that the average dielectric constant on the protein–water interface is approximately 20–30, yielding γ values of ∼5–7.5, which are in excellent agreement with our numerical result.

More importantly, the phenomenon of weakening the binding affinity of the RBD–CB6 complex observed in the K417N mutation is not unique. After comparison of some crystal structures for the complex of the RBD (of the SARS-CoV-2’s spike protein) and human antibodies, surprisingly we noticed that K417 in the RBD can form a buried salt bridge with either a glutamate or an aspartate in four other human antibodies as highlighted in IGHV-53 (Figure a), C1AB3 (Figure b), BD-236 (Figure c), and CC12.1 (Figure d). Thus, it is likely that all of these antibodies may not be able to neutralize the 501Y.V2 variant. It is worth noting that during the preparation of this paper, two experimental preprints posted on bioRxiv demonstrating that the South Africa variant can evade human antibodies CB6 (Figure c, also known as LyCoV016) and CC12.1 (Figure d) that are known to neutralize wild-type SARS-CoV-2.[13,25] These experimental results further confirm our predictions. Because K417 in RBD can form a buried salt bridge with a glutamate/aspartate in a class of human antibodies, we speculate that K417 is the Achilles’ heel of the virus and can be easily targeted by many human antibodies. Therefore, the K417N mutation could be evidence of viral adaptation to the human immune system or natural selection under the “pressure” of human neutralizing mAbs.

Figure 4

Illustrations of various salt bridges buried at the RBD–antibody interfaces. (a) Neutralizing antibody IGHV-53 (PDB entry 7JMO) for SARS-CoV-2’s RBD. (b) Neutralizing antibody C1A-B3 (PDB entry 7B3O) for SARS-CoV-2’s RBD. (c) Neutralizing antibody BD-236 (PDB entry 7CHB) for SARS-CoV-2’s RBD. (d) Neutralizing antibody CC12.1 (PDB entry 6XC2) for SARS-CoV-2’s RBD. The RBD, VL, and VH of various antibodies are colored gray, green, and purple, respectively. Residues forming a salt bridge at the RBD–antibody interface are shown as sticks. VL and VH of various antibodies are shown transparently to allow viewing of buried salt bridges.

Discussion and Conclusions

Although it is counterintuitive to learn that the K417N mutation actually decreases the RBD’s binding affinity with ACE2 by 1.48 kcal/mol, it is in fact beneficial for RBD to escape human antibody CB6 by dramatically reducing the binding affinity by 9.59 kcal/mol. Because of the loss of binding affinity, the variant with this single mutation is less competitive than the wild-type SARS-CoV-2 virus, in terms of the ability to enter host cells. So far, there is no evidence that this particular mutation occurs alone in any known variants of SARS-CoV-2. In the South Africa variant (501Y.V2), actually there are two more mutations, N501Y and E484K. To explore their contributions, we carried out MD simulation for the RBD–ACE2 complex with all three mutations (K417N, N501Y, and E484K), as shown in Figure .

Figure 5

MD simulation of the RBD–ACE2 complex with the N501Y, K417N, and E484K mutations in the RBD. (a) Simulation system. ACE2 (blue) and the RBD (gray) are shown in cartoon representation. The glycosylated residues are not shown. Na+ and Cl– ions are shown as yellow and cyan spheres, respectively. Water is shown transparently. (b) Equilibrated atomic coordinations around the three mutated residues. Protein residues are shown in stick representation. (c) Time-dependent distance D between E75 in ACE2 and K484 in the RBD. D is defined as the distance between the CD atom in E75 and the NZ atom in K484.

During the entire 340 ns MD simulation with the triple mutations (Figure a), Y501 interacted stably with residues Y41 and K353 in ACE2 (see Figure b and Movie S1). Namely, Y501 in the RBD formed a T-shaped contact with Y41 in ACE2, and simultaneously, Y501 formed a hydrophobic interaction with the alkane chain of the amphipathic K353. We did not observe any coordination formed by N417 in the RBD and any residue in ACE2. Before the E484K mutation, E484 forms a hydrogen bond with F490 (Figure S5), which attaches the loop (residues 477–486) to the remainder of the RBD. After the mutation, this hydrogen bond was abolished and the loop was able to detach from the remainder of the RBD. Interestingly, K484 was gradually pulled toward E75 in ACE2 due to the electrostatic attraction (see Figure b and Movie S1). In Figure c, we show that the initial distance D between K484 in the RBD and E75 in ACE2 was on average ∼11 Å. After ∼190 ns, the distance was decreased to ∼3.5 Å, suggesting the formation of a salt bridge. Due to the exposure to water, this salt bridge can break and re-form frequently. This new state persisted for the remaining 150 ns (Figure c). These dynamic interfacial interactions are illustrated in Movie S1, as well. Therefore, the gain from the improved RBD–ACE2 binding that resulted from both N501Y and E484K mutations might be enough to compensate for the loss caused by the K417N mutation and is likely to yield an interaction with ACE2 even stronger than that of the wild-type virus.

In summary, we have shown that K417 plays an important role in the binding between the RBD and ACE2 or CB6 by forming interfacial salt bridges. The salt bridge between K417 in the RBD and D104 in CB6 is buried inside the complex and is therefore much more stable than the water-exposed one formed between K417 in the RBD and D30 in ACE2. Thus, the K417N mutation can weaken the binding of the RBD with CB6 much more than with ACE2. More dramatically, the K417N mutation allows the variant to escape from many human antibodies other than CB6 by removing a salt bridge buried in the RBD–antibody interface. Interestingly, the virus with the K417N mutation seems to sacrifice its binding affinity for ACE2 to survive the attack of the antibodies. In fact, after the sequence alignment, the residue in SARS-CoV corresponding to K417 in SARS-CoV-2 is valine (V404).[26]

On the contrary, a recent study has demonstrated that the human immune systems may also quickly respond to these new variants and produce corresponding neutralizing antibodies that contain more somatic hypermutation, increased potency, and resistance to RBD mutations.[27] The continued evolution of the humoral response to recent SARS-CoV-2 variants warrants further studies. Additionally, with the revealed mechanism underlying the K417N mutation, it is possible to design more efficacious antibody cocktails to treat COVID-19 patients infected with variant 501Y.V2 as well as the recently discovered variant P.1 in Brazil.

Experimental Section

MD Simulations

We carried out all-atom MD simulations for the bound (the RBD–ACE2 complex) and free (stand-alone RBD) states using the NAMD2.13 package[28] running on the IBM Power Cluster. To model the RBD–ACE2 and RBD–CB6 complexes, we first obtained the crystal structures with PDB entries 6M0J(29) and 7C01,[30] respectively, and then solvated them in a rectangular water box that measured ∼80 Å × ∼80 Å × ∼135 Å. Bound Zn2+ ions were included in the RBD–ACE2 complex. Na+ and Cl– were added to neutralize the entire simulation system and set the ion concentration to 0.15 M (Figure b,c). Each built system was first minimized for 10 ps and further equilibrated for 1000 ps in the NPT ensemble (P ∼ 1 bar, and T ∼ 300 K), with atoms in the backbones harmonically constrained (spring constant k = 1 kcal mol–1 Å–2). The production run was performed in the NVT ensemble, and only atoms in the backbones of ACE2 that are far from the RBD (residues 110–290, 430–510, and 580–615 in ACE2) were constrained, preventing the whole complex from rotating out of the water box. The similar approach was applied in the production run for the RBD–CB6 complex, with atoms in the backbones of the CH (residues 122–218) and CL (residues 109–215) domains constrained. We also performed MD simulation for the RBD alone in the 0.15 M NaCl electrolyte (a free state) using the same protocol (Figure S1).

To simulate the binding between the RBD of the South Africa variant and ACE2, we obtained the crystal structure (PDB entry 6M0J) for the RBD–ACE2 complex and introduced three mutations (N501Y, K417N, and E484K) when preparing the .psf and .pdb files for the complex. We retained all glycosylations in both the RBD and ACE2 as well as one Zn2+ bound on ACE2 resolved in the crystal structure. The complex was solvated in a water box that measures ∼131.8 Å × ∼131.8 Å × ∼131.8 Å; 231 Na+ and 208 Cl– ions were added to neutralize the whole system and set the ion concentration to 0.15 M. The final simulation system (Figure a) comprises 232337 atoms.

After equilibrating the structures in bound and free states, we carried out free energy perturbation (FEP) calculations.[21] In the perturbation method, many intermediate stages (denoted by λ) whose Hamiltonian H(λ) = λHf + (1 – λ)Hi are inserted between the initial and final states to yield a high accuracy. With the softcore potential enabled, λ in each FEP calculation for ΔGA or ΔGB varies from 0 to 1.0 (forward) in 20 perturbation windows (lasting 300 ps in each window), yielding gradual (and progressive) annihilation and exnihilation processes for K417 and N417, respectively. We tested the convergence of the forward FEP runs with the backward ones and with the forward FEP runs with a longer simulation time of 1.5 ns per window; the obtained results are similar (Figure S3). The same approach has been applied in our previous work in which the obtained FEP results are consistent with experimental ones.[31] We followed our previous FEP protocol[20] to obtain the means and errors for ΔGA and ΔGB from the lowest five energy changes among ten data in total. In Figure S4, we also show that the means and errors for ΔGA and ΔGB obtained from all ten data are similar.

In FEP runs for the K417N mutation, the net charge of the MD system changed from 0 to −1 e (where e is the elementary charge). It is important to have similar sizes of the simulation systems for the free and bound states,[32,33] so that the energy shifts from the Ewald summation approximately cancel out when calculating ΔΔG.

We applied the CHARMM36m force field[34] for proteins, the TIP3P model[35,36] for water, and the standard force field[37] for Na+, Zn2+, and Cl–. The periodic boundary conditions (PBC) were used in all three dimensions. Long-range Coulomb interactions were calculated using particle-mesh Ewald (PME) full electrostatics with a grid size of ∼1 Å in each dimension. The pairwise van der Waals (vdW) energies were calculated using a smooth (10–12 Å) cutoff. The temperature T was maintained at 300 K by applying the Langevin thermostat,[38] while the pressure was kept constant at 1 bar using the Nosé–Hoover method.[39] While the SETTLE algorithm[40] enabled us to keep all bonds rigid, the simulation time step was set to 2 fs for bonded and nonbonded (including vdW, angle, improper, and dihedral) interactions, and electric interactions were calculated every 4 fs, with the multiple-time step algorithm.[41]