Point mutations in SARS-CoV-2 variants induce long-range dynamical perturbations in neutralizing antibodies

Abstract

Monoclonal antibodies are emerging as a viable treatment for the coronavirus disease 19 (COVID-19). However, newly evolved variants of the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) can reduce the efficacy of currently available antibodies and can diminish vaccine-induced immunity. Here, we demonstrate that the microscopic dynamics of neutralizing monoclonal antibodies can be profoundly modified by the mutations present in the spike proteins of the SARS-COV-2 variants currently circulating in the world population. The dynamical perturbations within the antibody structure, which alter the thermodynamics of antigen recognition, are diverse and can depend both on the nature of the antibody and on the spatial location of the spike mutation. The correlation between the motion of the antibody and that of the spike receptor binding domain (RBD) can also be changed, modulating binding affinity. Using protein-graph-connectivity networks, we delineated the mutant-induced modifications in the information-flow along allosteric pathway throughout the antibody. Changes in the collective dynamics were spatially distributed both locally and across long-range distances within the antibody. On the receptor side, we identified an anchor-like structural element that prevents the detachment of the antibodies; individual mutations there can significantly affect the antibody binding propensity. Our study provides insight into how virus neutralization by monoclonal antibodies can be impacted by local mutations in the epitope via a change in dynamics. This realization adds a new layer of sophistication to the efforts for rational design of monoclonal antibodies against new variants of SARS-CoV2, taking the allostery in the antibody into consideration.

---

Introduction

The coronavirus disease 19 (COVID-19) caused by the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2)¹ has claimed, as of March 2022, more than 6 million lives worldwide, creating a global pandemic and one of the largest public health crises in human history.² To address this urgent problem, scientists across various disciplines are striving to develop drugs, vaccines, and antibodies against the virus.^3–7 The spike proteins, present on the surface of the virus, recognize and bind to the human angiotensin converting enzyme 2 (hACE2) receptor in lung cells and initiate infection.^8,9 The receptor binding domain (RBD) of the spike is a key target for drug development and antibody recognition.^10–16 A large number of monoclonal antibodies (mAb), as well as the natural antibodies (Ab) generated by the immune system, block infection by binding to the RBD and prevent it from attaching to the hACE2 receptor. Moreover, many of the currently available vaccines, such as the ones developed by Pfizer,¹⁷ Moderna,¹⁸ and AstraZeneca,¹⁹ use the spike protein as their epitope and induce immunity by generating antibodies and memory cells that can potentially recognize regions of the spike protein, including the RBD.

The process of antigen–antibody recognition is of fundamental importance for developing immune response against invading pathogens. Since their discovery in 1975,²⁰ monoclonal antibodies (mAb) found a wide range of therapeutic applications, particularly in the treatment of cancer, chronic inflammation and viral infection.^21–24 Unlike natural antibodies which have varying sequences, monoclonal antibodies, engineered in the laboratory, have a single sequence, and are thereby tailored to treat specific diseases. They are generated by identical B-lymphocyte immune cells cloned from a unique parent white blood cell. The use of engineered antibodies as drugs has become increasingly effective due to their high binding affinity and antigen specificity, both modulated by the complementarity determining regions (CDR) within the variable heavy (V_H) and variable light (V_L) chains.²⁵

Monoclonal antibodies have significant advantages over conventional immunization techniques such as convalescent plasma therapy (CPT), because they reduce the risk of blood-borne diseases, shorten the time needed to generate high-affinity and low-risk antibodies, and allow for variable epitope specificity.²⁶ A number of monoclonal antibodies have been reported to neutralize the novel coronavirus and are at various stages of clinical trial, with bamlanivimab and REGN-COV2 having received emergency use authorization (EUA) in the United States of America by the Food and Drug Administration (FDA) as of June 2021.²⁷ Molecular structures of multiple monoclonal antibodies in complex with the spike protein of SARS-CoV-2 have now been resolved, facilitating a detailed structural-although not as yet a dynamical-understanding of the antigen recognition process.^28–35

The major concern for these treatments is the possibility of the appearance of new variants, given mutations in the RBD region, which can evade immune response and reduce or nullify the action of neutralizing antibodies. Multiple such variants of SARS-CoV-2 have emerged throughout the course of the pandemic, threatening the effectiveness of current therapeutics against the coronavirus. One of the earliest well-documented variant, with mutations in the RBD, was the B.1.1.7 strain, a.k.a. the alpha variant. It emerged in the United Kingdom in the second half of 2020 and spread across the entire world vigorously because of its 70–80% higher transmissibility compared to the wild type.³⁶ This variant contains a N501Y mutation in the RBD, which has been attributed to its increased immune evasion capability, although no significant increase in reinfection was observed.³⁷ The next major highly infectious variant was the beta variant or the B.1.351 strain, originating in South Africa in October 2020. Apart from N501Y, it has two more mutations in the RBD: K417N and E484K; this mutant reduces the effectiveness of conventional antibody therapy up to 10-fold.³⁸ Later, a new strain (B.1.617) with two RBD mutations, L452R and E484Q, emerged in India and became officially known as the kappa variant. A lineage of that strain, B.1.617.2, is named as the delta variant by the World Health Organization (WHO). With 50–67% higher infectivity than the alpha strain,³⁹ the delta variant is designated as a variant of concern (along with alpha and beta) because of its potential to make currently available vaccines less effective in preventing serious illness. Evidence of serious decline in vaccine effectiveness against infection (VE-I) has already appeared, along with the proliferation of the delta variant.^40,41 Single molecule force spectroscopy experiments, in combination with atomistic simulations, confirmed the direct relation between the change of the receptor binding affinity and mutations in the spike protein.⁴² The same study also demonstrated the reduction of spike protein neutralization by immunoglobulin G based monoclonal antibodies such as B-R41.

To address this urgent crisis, a number of computational studies have been performed to obtain a molecular level understanding of the effect of mutations in the RBD on its binding with the hACE2 receptor and antibodies.^43–49 A recent work, specifically explored the change in binding free energy and molecular interactions between spike RBD and a range of neutralizing nanobodies and monoclonal antibodies.⁵⁰ Taken collectively, these studies revealed that the residues which got mutated in the variants of concern can play an important role in antibody binding. In addition to directly computing binding affinities, some of these studies proposed the possibility of allosteric communication between the RBD and the receptor or the antibody,^44,46,49 and showed that dynamical perturbations within the hACE2 receptor can arise due to its binding to the RBD. At the same time, simulations have shown that (i) the dynamics of the glycan coat covering the spike affects antibody binding, and (ii) that antibody binding itself affects the motion of the RBD.⁵¹

From the results of these in silico studies, it can be argued that the dynamical changes and allosteric effects can play a key role in determining the relative stability of antigen–antibody complexes or the efficacy of monoclonal antibodies. Although allosteric modulation within the spike trimer has received significant attention throughout the pandemic, the viewpoint that mutations in the spike protein can induce allosteric perturbations inside neutralizing antibodies has remained largely unexplored. This is important because such an allostery can facilitate immune evasion. In the present work, we explore the signatures of RBD mutations in the intrinsic dynamics of the monoclonal antibodies when in conjugation with the spike protein. Such effects can only be manifested in atomistic resolution, making them elusive to most common experimental techniques. To understand these allosteric communication mechanisms in greater detail, we performed extensive all-atom molecular dynamics (MD) simulations to gauge the effect of RBD mutations on the stability, dynamics and the unbinding process of monoclonal antibodies targeted to the SARS-CoV-2 spike protein. Apart from the wild type RBD, we studied four mutant strains with the following mutations: B.1.1.7 (alpha) (N501Y), B.1.351 (beta) (K417N, E484K, N501Y), B.1.617 (kappa) (L452R, E484Q) and B.1.617.2 (delta) (L452R, T478K). The two antibodies, B38 (ref. 31) and BD23,³² were chosen for the current study because their epitopes contain at least two of the four residues involved in the mutations in the studied variants and also because the new mutations in the beta variant have been shown to increase the binding affinity for B38 but decrease it for B23.⁴³ This contrasting property makes it particularly interesting to understand how the molecular interaction, and the allosteric cross-talk with the RBD, differ between these two mAbs.

Both of these antibodies were first extracted from the plasma of COVID-19 convalescent patients.^31,32 Each is a dimeric protein comprising a heavy and a light chain (labeled as chain A and B, respectively, in this manuscript) with clearly distinguishable antigen binding domains (Fig. 1). The BD23-Fab is smaller (126 and 107 residues in the heavy and light chain respectively) compared to the B38-Fab (222 and 218 residues in the heavy and light chain respectively). The epitope of the B38 contains the RBD residues 417 and 501 while the residues 452, 484 and 501 are included in the epitope of BD23. The structures of the RBD–antibody complexes as well as the location of the mutations are depicted in Fig. 1.

Fig. 1 The structure of the RBD in complex with BD23 and B38-Fab, the two neutralizing antibodies studied in this work. Each antibody consists of a heavy (A) and a light chain (B), with chain A in green, chain B in yellow, and the RBD in cyan. The antigen binding domain (ABD) of the Ab is defined as the antibody domain that comes in contact with the RBD; the rest of the antibody is referred to as the distant domain (DD). The glycan moiety on Asn343 is colored in licorice. The residues that are mutated in the studied variants are represented as colored spheres: K417 (red), L452 (orange), T478 (black), E484 (blue) and N501 (pink).

A number of computational studies in the past two years have demonstrated the significance of the glycan shield of the spike protein in receptor binding and immune evasion.^52–55 Casalino et al. showed that the roles of glycosylations can go beyond shielding and can alter the conformational dynamics of the spike trimer.⁵² Such conclusions were reinforced by the work of Sztain et al. where a glycan moiety was identified to work as a gate in allowing the “down” to “up” conformational change in the spike RBD, a mechanistic step, that is essential for SARS-CoV-2 virus to initiate infection.⁵³ As the binding of the hACE2 receptor and neutralizing antibodies generally takes place in the “up” conformation, the effect of glycan moieties in the recognition of the RBD can be less direct. Indeed, Nguyen et al. demonstrated the residue contact signature between the RBD and ACE2 remains unaffected by the presence of the spike glycans. The role of ACE2 glycans, on the contrary, is much more significant in stabilizing the spike receptor complex.⁵⁴ In the present work, we focused on the role of individual mutations on the dynamics of neutralizing antibodies; so we carefully picked two antibodies whose epitopes inside RBD do not include glycosylated residues. Nevertheless, we explicitly modeled the glycan on Asn343, although it has no interaction with either of the antibodies during the course of the simulation.

We studied the effect of RBD mutations on the dynamics of the antigen–antibody complex from three different angles. First, we demonstrated the relative flexibility of individual residues across the antibody, and identified the structural changes taking place in the complex when the RBD of the WT coronavirus is replaced with each of the mutants. Next, we quantified the change in the mutual information-based cross-correlation between all pairs of residues and constructed a protein graph connectivity network to decipher how are the allosteric information pathways modulated by the mutations in the RBD. Finally, we studied the change in the dissociation mechanism of the antigen–antibody complex for all combinations of the mutant RBDs using non-equilibrium pulling simulations, and rationalized our finding based on arguments involving enthalpy and hydrogen bonding. From the diverse dynamical perturbations observed in these systems, we delineated a few common themes of how the mutations in the SARS-CoV-2 antigen can impart local and global effects through the antibody structure and potentially impair the neutralizing efficacy of monoclonal antibodies. This can occur by changing interaction energies, as well as, via large-scale motions, the conformational entropy of the binding.

Results

Modulation of conformational flexibility of RBD–antibody complex

The RBD, with or without the antibody, is a relatively rigid system (RMSD ∼ 1–3 Å) with little conformational variability, likely due to the abundance of beta-sheets in the structure. The antibodies, conversely, are highly flexible with large conformational fluctuations, as indicated by their high root-mean-squared deviation (RMSD) values throughout the production runs (see ESI Fig. S3 and S4†). Unlike the BD23 system where the entire antibody is highly flexible, in the case of B38, the fluctuations are significantly larger in the distant domain (DD), compared to the antigen binding domains for both chains. The RBD–Ab complexes of both the WT and mutants are also capable of exploring more than one free energy minimum during a period of 0.4 μs, as demonstrated by the projection of their dynamics on the first two principal components (ESI Fig. S6 and S7†). However, any discernible difference between the dynamics of the WT and mutant complexes is not apparent from the RMSD and PCA analysis.

In contrast, the root-mean-squared-fluctuations (RMSF) of certain residues of both the RBD and Ab do show significant change when comparing the WT to the mutants (ESI Fig. S5†). To quantify this effect, we computed the ratio of RMSF of each residue for every mutant with that of the WT system. Even in the case of the RBD-only system, the region comprising residues 460–500 shows increased flexibility in the L452R–E484Q variant and decreased flexibility in the N501Y mutant. Interestingly, the E484Q mutation in the kappa variant, and the T478K mutation in delta variant are present within this region, which also interacts directly with both the B38 and BD23 antibodies. The change in flexibility upon mutation is the highest in the loop formed by residues 470–490, going up or down by a factor of 2–3 in the B.1.617 or the B.1.617.2 strains, respectively. The individually mutated residues show a sharp change in RMSF values, particularly residues 452, 484 and 501 in the respective variants. This effect is also visible in the antibody-bound systems, where the mutated residues show sharp peaks in the RMSF ratio plot (Fig. 2 and 3). Residue 484 shows a contrasting behavior in the presence of B38 and BD23 antibody. When in complex with B38-Fab, K/Q484 and adjacent residues show a large increase in flexibility (∼2 times that of E484 in WT), but they become almost twice more rigid when bound to BD23. The RBD of kappa variant shows the largest increase in rigidity compared to WT, on the average, when bound to either antibody, but becomes the most flexible of all without antibody. The N501Y, the K417N-E484K-N501Y and the L452R-T478K mutants show increase in overall flexibility when in complex with the antibody. But their RMSF remain virtually unchanged with respect to the WT RBD in the absence of any antibody.

Fig. 2 Change of root mean squared fluctuation (RMSF) of the RBD-B38 Fab complex due to the mutations in the variants. The RBD is represented by a transparent surface. The mutated residues are shown as cyan spheres. The antibody residues are colored according to the ratio of their RMSF when bound to a mutant RBD with that of the WT case. Blue color indicates increase in RMSF and red color indicates decrease in RMSF. The residue RMSF ratio is projected on the structure of the B38–Ab complex for alpha (N501Y), beta (K417N–E484K–N501Y), kappa (L452R–E484Q) and the delta (L452R–T478K) variants. The RMSF ratio for the RBD and the antibody are shown separately in the right column. See text for the definition of the RMSF ratio. In the RBD RMSF ratio figures, the location of the mutations are shown as vertical dashed lines.

Fig. 3 Same as Fig. 2 but for the RBD-BD23 Fab complex.

More interesting is the change in RMSF in the antibody, because this change indicates that mutations in the antigenic epitope can influence the microscopic dynamics of the antibody it binds. The change in the residue-wise RMSF of both the antibodies are, as expected, less dramatic when compared to that in the RBD, likely because these changes arise from a secondary perturbation caused by mutations outside the protein domain (i.e., via non-covalent interactions). Fig. 2 and 3 show the ratio of the RMSF of the antibodies in complex with the mutant RBDs with respect to their WT counterparts. For B38-Fab, the deviation is small (RMSF ratio within 0.5–1.5), except for specific residues in the distant domain, when bound to the L452R–E484Q mutant. However, the changes in flexibility are much more prominent for the BD23 antibody, for which RMSF increases on average in response to a N501Y mutation in the alpha variant, but the effect is largely compensated by the K417N and E484K mutations in the beta variant, bringing the RMSF down to a level comparable to a WT complex. The L452R and E484Q mutation in the B.1.617 strain causes an increase in overall rigidity of the antibody structure, particularly in the distant domain. The change of antibody RMSF for the delta variant case is almost identical to the kappa variant, except the antigen binding domain for BD23 shows a slight increase of flexibility in comparison to the kappa variant.

Change in correlated motion and allosteric information flow

Significant changes in the residue-wise LMI-based cross-correlation can be observed throughout the complex (both RBD and antibody) in presence of the RBD mutations (see ESI Fig. S11 and S13†). We particularly focus here on the cross-correlation coefficient between the RBD and the antibody residues, as this can provide a direct quantitative measure of the degree of dynamical coupling between them. In Fig. 4 we depict the relative RDB–antibody cross-correlation for each mutant system, with respect to the WT complex. One of the key features of the B38-Fab is that Ab residues show large change in correlation, particularly with the mutated RBD residues and the residues adjacent to the mutation. For example, in the N501Y variant, many residues of B38 Ab show a decrease in correlation, by 0.1–0.2, with the residue 501 and its nearby RBD residues. More interestingly, this effect is propagated over a long distance, even to the distant domains of the antibody, indicating that single point mutations in RBD can have a strong long range perturbative effect on the dynamics of the antibody bound to the RBD. A similar phenomenon is observed for N417 and its neighboring residues (410–420) in the K417N–E484K–N501Y mutant, the Q484 residue in the case of the L452R–E484Q mutant, and the R452 and its nearby residues (440–455) in the L452R–E484Q and L452R–T478K mutant. Although this effect is less pronounced in the BD23-Fab, the change in cross-correlation coefficients with the mutated RBD residues is indeed observed when, for example, in complex with the K417N–E484K–N501Y mutant. Unlike B38, this effect is less uniform within the antibody, where certain groups of residues (e.g. near residue 50, 70 and 100 in the ABD), show large changes in correlation with the mutated residues as well as for the whole RBD. This can be visualized as vertical patches of red or blue color, as opposed to horizontal patches for B38.

Fig. 4 Change in the mutual information based cross-correlation between the antibody and the RBD residues, due to the mutations in the variants. The differences are reported as the cross-correlations of the WT species subtracted from that of the mutant species. The locations of specific mutations are labelled as dashed horizontal lines. The continuous vertical lines separate different domains and chains of the protein (for details, see Fig. 1).

Moreover, there is, on average, an increase in the cross-correlation between the RBD and the B38 antibody, going from the alpha variant to the beta variant and vice versa in BD23. Previous computational studies have established that B38 binds more weakly with the spike of the alpha variant compared to the beta variant, but BD23 binds at least twice as strongly to the alpha variant than to the beta variant.^43,45 This observation indicates that the mutual information cross-correlation and the binding affinity between the monoclonal antibodies and the spike epitopes follow a similar trend.

This opposite effect in the correlated motion and the binding affinity suggests that there are allosteric communication paths between the RBD and antibody that follow different routes in the two antibodies. To probe these pathways, we built a protein graph-connectivity network model, treating each residue as the node and the negative logarithm of the cross-correlation between the residues as the weights of the edges of the graph. We computed the betweenness centrality (BC) measure, a metric which computes the relative importance of each node in the information flow pathway. We then computed the residue wise difference of BC in each mutant relative to the WT system. Changes in BC point us to the residues which either gain or lose importance in the allosteric information flow pathway in the protein. In B38-Fab, significant changes in the BC metric are primarily limited to the residues of the antigen binding domain, showcasing a situation where local correlation propagates through the amino acid sequences, leading to long range dynamical changes. Particularly, we see a large increase in the BC value for a few consecutive residues in the region of residues 92–98 in chain A, and a large decrease in the BC for residues 23–30 in chain B, for all the mutant species (Fig. 5). This imbalance is primarily caused by the differential interaction of the 501 residue in its WT (N501) compared to its mutated form (Y501). The change in BC in chain B is reduced by a significant amount when the N501Y mutation is removed, for example in the L452R–E484Q variant. This effect can be a result of the reduced interaction energy of Y501 with the B38 mAb as well as the loss of the hydrogen bond between the RDB N501 residue and the S30 residue of chain B of B38, although, there is an overall increase in the hydrogen bond occupancy for Y501 as the lost bond is replaced by a hydrogen bond with S67. In essence, the allosteric information flow through the B38 antibody follows a re-wiring mechanism. The mutations switch off the information flow pathway along the residues 23–30 in chain B and increases the information traffic along residues 92–98 in chain A. This effect is highly asymmetric and location specific, and leads to a drop or rise of the cross-correlation of all antibody residues with particular mutated residues in RBD.

Unlike the B38, in BD23-Fab a large fraction of the RBD–mAb interface (specifically the region containing residues 97–107 of the ABD) loses centrality, along with some isolated, distant domain residues; this theme is common across all mutants (Fig. 5). Reduction of BC is also observed for residues 40–50 in the ABD and residues 80–90 in the DD, indicating the possibility of long-distance propagation of an allosteric signal. However, the dynamical effect is different from that in the B38, as the perturbation is not localized in the vicinity of any specific mutation. The same set of residues experiences a decrease in the centrality measure, irrespective of the specific location and the nature of the mutation (except for the alpha and delta variants, where some residues show a gain in BC, particularly the ABD 45–50 and DD 85–90). This explains qualitatively the observation of the mutation-specific long distance correlation changes in the case of B38 (horizontal correlation changes in Fig. 4). Contrarily, in BD23, specific residues lose or gain correlation with the majority of the RBD residues and not particularly with the mutated ones (see the vertical correlation changes in Fig. 4).

Taken together, these findings provide fundamental insight into the molecular mechanism of how mutations in the spike protein of SARS-CoV-2 variants impact antibody recognition and the stability of the antigen–antibody complex. Not only do these observations bolster the necessity to develop vaccines against mutated or alternative epitopes, but also focus attention on specific regions of the antibodies that can be mutated or modified to engineer a more robust neutralizing antibody, capable to act against potential future variants of the virus.

Fig. 5 Change of betweenness centrality (BC) of each residue of the RBD–Ab complex due to the mutations in the variants. The RBD is represented by a transparent surface. The mutated residues are shown as cyan spheres. Only residues with a change in BC of greater than 0.005 is shown in colored sphere. Blue color indicates increase in BC and red color indicates decrease in BC.

Dissociation dynamics and interactions at the RBD–antibody interface

The observation of allosteric effects within the antibody structure resulting from mutations in the spike also motivates further investigation into the role of mutations in antibody binding and unbinding processes. Considering the difficulty of simulating the diffusion-assisted binding of the antibody in a biologically-accurate environment, we only investigated the possible alteration of the antibody detachment mechanism due to the mutations in spike variants. This provides an alternate angle to understand the role of antigen mutations in the stability of the antigen–antibody complex from a more dynamical perspective. From multiple steered molecular dynamics simulations (SMD), we elucidated the unbinding mechanism of the monoclonal antibodies from the spike RBD. For detailed inspection, we picked the pathway that corresponds to the lowest amount of non-equilibrium work for each RBD variant and mAb combination, as it would contribute the most to the real unbinding process. For all variants and both the antibodies, the loop consisting of residues 470–490 is found to play an important role in the unbinding process. It is the last point of contact between the RBD and the antibody in the lowest work dissociation pathway for every system, and it functions like an anchor to adhere to the antibody before a complete detachment occurs. A recent computational study revealed that this loop region also functions like an anchor while binding to the hACE2 receptor.⁵⁶ Unsurprisingly, one of the key mutating residue (E484) is situated in this region of the RBD, and, when not mutated, it is often the last RBD residue that forms a hydrogen bond with the dissociating antibody. But in the K417N–E484K–N501Y and L452R–E484Q variants, this residue is mutated and no longer forms the “last” hydrogen bond to anchor the departing antibody. Instead, other residues from the same loop, primarily S477, which has a hydroxyl side-chain, perform this role. We also observed A475 to form the “last” hydrogen bond for the L452R-E484Q RBD with B38 antibody, with its main chain hydrogen bond donor group (see ESI Fig. S24†). Although the 484 residue is not mutated in the delta variant, the last hydrogen bonds between the delta RBD and the departing antibody is formed by the A475, S477 and E486 residues.

We computed the average non-equilibrium work required to detach the antibody for all systems, expressed as −ln〈exp(−W_i)〉, where W_i denotes the accumulated work in the ith SMD run, and 〈〉 is the arithmetic mean. The average work is heavily dominated by the dissociation path that requires the lowest amount of work. In the limit of an infinite number of pulling trajectories or in the limit of infinitely slow pulling, the average non-equilibrium work will become equal to the unbinding free energy. As of the convergence of average work, we want to emphasize that we do not aim to calculate an accurate non-equilibrium work, instead we merely generate samples of work along a handful of paths. So the average work reported in this work is not meant to be used as an estimate of the free energy (as done, for example, when using the Jarzynski identity⁵⁷ or other fluctuation theorems⁵⁸). To establish the validity of our implicit solvent SMD approach, we performed additional steered molecular dynamics simulations in an explicit solvent environment. While it is true that explicit solvent gives more accurate thermodynamics, when pulling a ligand through water the rate of pulling must be accounted for.⁵⁹ Implicit solvents permit instantaneous (“adiabatic”) rearrangement of the “water”, while explicit water imposes a limit to the speed with which the ligand can be pulled (i.e., slower than the structural relaxation of the explicit water molecules). As a consequence, the absolute work will be different due to interaction with the water molecules during the pulling. However, the relative differences are similar, as we show in the additional calculations (see ESI†). Also, the calculation of uncertainty in exponential averaging of non-equilibrium work can be non-trivial due to the presence of finite statistical bias.⁶⁰ Nevertheless, we showed a very primitive mean-square-error based uncertainty analysis in the ESI (Fig. S22 and S23†). Although we can observe different relative uncertainty across different variants, the absolute values of the variance are not particularly meaningful due to the small sample size.

Provided E484 is the key anchoring residue, we expect the non-equilibrium work will depend on how strongly this particular residue or its mutants interact with the antibody. From the non-bonded energy calculation on the equilibrium trajectory it is apparent that the three mutated residues: E484, K484 and Q484, have little difference in their interaction energy with the B38-Fab antibody, primarily because this residue is not in direct contact with the antibody interface (Fig. 6b). Consequently, the four mutants do not show any discernible difference in the work required for the dissociation of the RBD–mAb complex (Fig. 6c). But, the alpha and beta mutants require the least amount of work for antibody detachment in comparison to the WT, an effect that can partially be attributed to the N501Y mutation. From non-bonded energy calculation, we see that the N501Y mutation destabilizes the interaction between RBD and B38 by lowering (less negative) the interaction energy (see ESI Fig. S16†). Non-bonded interaction energies have previously been used to understand the impact of mutations in the receptor binding affinity of the SARS-CoV-2 RBD.⁴² Similar to the binding of hACE2 receptor, we observed that the change in interaction energy can also have a significant effect on the binding of antibodies to the RBD.

Fig. 6 (a) A representative intermediate structure for the B38 antibody unbinding from the RBD from the SMD simulation. Intermediates for all variants are shown in ESI Fig. S24 and S25.† (b) The non-bonded interaction energy between the RBD residue 484 and the B38 antibody from the unbiased simulation. (c) The average non-equilibrium work performed to dissociate the B38 antibody from all variants of the RBD using SMD simulation. (d–f) Are identical as (a–c) but for BD23 antibody.

For the BD23-Fab, an almost opposite effect is observed. The residue E484, being situated in the RBD–mAb interface, plays a key stabilizing role for the complex via hydrogen bond to the S103 and S104 residues of the Ab. Once mutated, the hydrogen bond propensity of residue 484 diminishes, leading to almost 80 kcal mol⁻¹ decrease (less negative) in the non-bonded interaction energy with the mAb, in agreement with a recent computational study.⁴⁵ The E484 forms hydrogen bonds with the S104 (26.37% with backbone and 54.19% with side-chain) and S103 (3.85% with backbone and 8.62% with side-chain) residues of BD23 with high occupancy in the WT and with G102 (31.9%), S103 (4.15%), Y32 (10.7%), Y110 (3.95%) and Y111 (2.73%) residues of the antibody in the N501Y mutant. In the K417N–E484K–N501Y variant, K484 forms hydrogen bonds only with N92 (6.86%) and Y110 (1.18%) with relatively low occupancy. This effect is clearly reflected in the non-equilibrium work profile, where, going from N501Y mutant to K417N–E484K–N501Y and L452R–E484Q variants, we see a significant drop in the work required for the detachment of the antibody. Conversely, for the delta (L452R–T478K) variant, where the 484 residue is not mutated, the E484 forms hydrogen bonds with >50% occupancy with the S103 and Y109 residue of the BD23 antibody. This increases the amount of non-equilibrium work required for the dissociation of the RBD antibody complex to a strikingly high level. This is reflected also in the high binding free energy of the BD23 antibody with the delta variant RBD.

In contrast, the N501Y mutant, binds more strongly with the BD23-Fab (Fig. 6f) compared to the WT, because the non-bonding interaction energy of Y501 is higher (more negative) than that of the N501 residue (ESI Fig. S17†). This effect can be attributed to the increased probability of hydrogen bonding between the RBD and the antibody by the Y501 residue compared to the N501. In the WT system, N501 participates in one hydrogen bond with the D73 residue of the BD23-Fab with negligible occupancy (0.62%). But in the N501Y mutant, the Y501 participates in two hydrogen bonds with the L72 and the Y80 residues of the antibody with significantly higher occupancy (5.77% and 6.27% respectively). The rupture forces of the spike antibody complexes follow a similar trend as the average non-equilibrium work for the WT and the mutants (see ESI Fig. S20 and S21†).

The relative order of the binding free energies for each RBD antibody pair is also consistent with the average non-equilibrium work necessary for the dissociation of the respective complexes (Fig. 7). In agreement with the earlier work,⁴³ the binding affinity of the BD23 with the beta variant was observed to be weaker than the alpha strain, and vice versa for the B38 Fab monoclonal antibody. This indicates that the arguments presented above, concerning the non-bonded interaction energy and hydrogen bonding propensity, also hold true for the binding free energy of the antigen–antibody complexes. The practical implication of these findings is that antibodies, such as B23, that uses the RBD 470–490 loop as an epitope, are more likely to have significantly different affinity against mutated versions of the virus. Alternatively, fairly consistent neutralization of most variants is possible if antibodies or vaccines are generated using an alternative region of the receptor binding domain as the epitope.

Fig. 7 Binding free energy of the spike protein variants with the B38 and BD23 antibody.

Discussions and conclusions

The simulations of the interaction between the receptor binding domain of the SARS-CoV-2 spike protein and neutralizing monoclonal antibodies reveals the functional role of point mutations in the new variants in the process of immune evasion. Point mutations perturb the microscopic dynamics of the antibodies and destabilize the antigen–antibody complexes. Key findings of our work are the following. First, the structures of the antibodies are largely flexible and can undergo significant conformational transitions over a timescale of several hundreds of nanoseconds. However, the flexibility of certain residues in the antigen–antibody complex are modified disproportionately when the WT RBD is replaced with one from the mutant strains. The mutated RBD residues (and the residues in proximity with the mutations) show the largest change in flexibility. Particularly, in the absence of an antibody, the loop involving residues 470–490 becomes flexible when a Glu to Lys (or Gln) mutation at the 484 position makes it incapable of forming a stabilizing hydrogen-bonding interaction. While local motions in the antibody occurred next to the mutation site in the spike, there were also significant large-scale collective motions throughout the antibody structure. We speculate that these motions can modulate the free energy of binding via contributions to the conformational entropy. Changes in binding can also result from the overall change of stability of the antigen–antibody complex, arising from the noticeable changes in rigidity observed not only in the antigen binding domain, but also in the distant regions. The effect is more prominent for the smaller antibody, BD23 Fab, where the N501Y mutation drastically increases the dynamic motion across both the domains in the antibody structure. The K417N and E484K mutations reverse this effect to a large extent, and, in the B.1.617 variant, both antibodies become structurally more rigid compared to the WT complex. Because the dynamical perturbations in the antibody result from mutations in the antigen, we are presented with a striking example of inter-protein allosteric information propagation, which has the potential to affect the evolutionary driving forces on the virus to adapt against neutralizing interventions.

The linear mutual information-based cross-correlation between the antigen–antibody residue pairs is a useful metric that can capture quantitative information about the allosteric information flow. Cross-correlation values, for both the RBD and the antibody, can change significantly due to mutations in RBD, but there is a stark contrast in the manifestation of this change in correlated dynamics in the two antibodies. Unlike B38-Fab, where correlations with specific mutated RBD residues are modified across the entire antibody, such effects in BD23 Fab are fairly localized albeit universal. Certain groups of residues in the BD23-Fab gain or lose correlation with the entire RBD and not just with the mutated residues. The antibody residues impacted by the mutation are conserved for all mutants, but whether they show a positive or negative change in cross-correlation can depend on the nature of the mutation. We previously observed such an effect within the spike protein,⁶¹ but now we see that the allosteric effect can be propagated into the antibodies that interact with the spike only via non-covalent interactions. One key observation is that the gain or loss in cross-correlation follows the trend of binding affinity fairly well. This phenomenon behooves one to investigate further, for a wide range of protein–protein complexes, the relationship that inter-residue cross-correlations has with, on one hand, binding free energies. On the other hand, simulation-based correlation paths can put in relationship with experimental measures of allosteric signals transduced within and between proteins. For example, using cryo-EM, Frank and collaborators have delineated a way to retrieve functional pathways of biomolecules from single-particle snapshots.⁶² Among the various experimental approaches that can probe allostery at the atomic resolution detail is liquid-state NMR spectroscopy. Using it, one can obtain site-specific dynamical information that can be used to probe allostery experimentally, see for example.⁶³ The dynamical parameters measured by the NMR experiment can also be calculated from the molecular dynamics simulation.⁶⁴ This strategy can serve as an interface with the experimental data to substantiate key findings on allostery.

A protein graph connectivity network provides a deeper understanding of the information flow pathway inside the antibody. With the aid of graph theoretical measures like betweenness or eigenvector centrality, one can quantify the importance of individual resides in the allosteric information propagation. We specifically looked at the change in betweenness centrality for the antibody residues induced by the mutations in the spike antigen. While a group of amino acids in the B38 light chain (proximal to RBD residue 501) experience centrality reduction, it increases in virtually the opposite end of the heavy chain in the RBD–antibody interface. Absence of the N501Y mutation reduces this effect, revealing a antigen variant induced re-wiring mechanism that decreases the traffic of allosteric communication inside antibody in the vicinity of the N501Y mutation and increases near the K417N and E484K/Q mutations. This effect being particularly specific to the location of the mutation in the epitope, it explains the mutation specific collective change in the cross-correlation throughout the B38 antibody. In contrast, all the BD23 residues near the RBD binding interface experience a reduced centrality, irrespective of the specific location of the mutations. Some distant residues, following the trend of interfacial residues, show a change in BC value without the regard of the location of the mutation, explaining the specific RBD–BD23 cross-correlation pattern. These contradictory features point to inherent differences in the microscopic dynamics of the antibody, depending on the nature of mutations in the antigen and on the nature of the antibody. They have repercussions beyond the current pandemic, controlling the fundamental process of antigen–antibody recognition at a molecular level. Such allosteric communications can be present within the hACE2 receptor upon binding to the mutants of the coronavirus spike protein. But, because of the prime role of the glycan moieties in controlling the binding affinity between the spike and the receptor, the demonstration of such effects from inter-residue correlations can be more subtle.

Whether allosteric modulation of monoclonal antibodies can significantly change the effectiveness of the immune response still remains to be fully understood. Nevertheless, we demonstrated that allosteric perturbations do take place within the antibody and such perturbations can be connected to the mutations in the receptor-binding domain of the spike protein. To see if these mechanisms can be generalized to a wider range of pathogens and antibodies, an analysis similar to ours can be applied to other antigen–antibody pairs. We hope to verify this hypothesis in future work.

The process of unbinding of the antibodies from the spike proteins of different variants is also of significant interest, as it provides a more mechanistic insight into the mutation induced modification of the antibody binding propensity. Irrespective of the variant or the antibody, the RBD loop, consisting of residues 470–490, works like an anchor and is always the last region to interact with the antibody in the lowest energy dissociation pathways. The E484K/Q mutation in this region reduces the anchoring capability with the BD23 antibody by the loss of key hydrogen binding interactions. So, the N501Y, K417N–E484K–N501Y, and L452R–E484Q mutants show a decreasing trend in the non-equilibrium work required to dissociate the BD23 antibody, and the recovery of the E484 residue in the delta variant regains the binding affinity. B38, in contrast, does not interact with residues 484 or 452 in the bound state, and mutations in those residues have little to no impact on the interaction energy of B38 with the spike protein. So, all the RBD variants in this study have almost identical binding affinity with the B38 antibody, as apparent from the steered MD results. But the mutants have lower “binding strength” than the WT because of the loss of stabilizing hydrogen bonds as a result of the N501Y mutation. It is interesting to note that the newly emerging highly contagious delta (L452R–T478K), delta plus (K417N–L452R–T478K) and omicron⁶⁵ (S477N, T478K, E484A, Q493K etc.) variants also gained different mutations in this important loop region that is vital in preventing the dissociation of the neutralizing antibodies.

In summary, we quantified how the microscopic dynamics of the antigen–antibody complex is modulated by the mutations in the antigen and the nature of the antibody, with a particular focus on the SARS-CoV-2 spike protein. Although there is significant diversity in the dynamical perturbations observed at the molecular level, there are common underlying themes: modulation of inter-residue correlations, allosteric information propagation, and destabilization of non-bonded interactions dictate the stability of the antigen–antibody complex, which in turn determines the efficacy of monoclonal antibody therapies. Taken together, this study explores the molecular signatures of the antigen-induced allosteric effects within the antibodies for coronavirus. The scope of our findings is, however, broader and future studies in this direction can lead to a deeper understanding of the immune evasion mechanism of mutating viruses to help in designing adaptive therapeutic strategies to induce immune responses against new viral strains via monoclonal antibodies and vaccination.

Methods

System preparation

The coordinates of the antibodies bound to the receptor binding domain (RBD) of the spike protein in wild type SARS-CoV-2, were obtained from the Protein Data Bank (PDB) with the following accession IDs: 7BZ5 (B38)³¹ and 7BYR (BD23).³² Apart from the B38 antibody, the 7BZ5 PDB file contains the RBD domain (only residues 334–526) including one N-acetyl-β-D-glucosamine (NAG) moiety attached to Asn343. This glycan moiety was preserved in our calculation. The 7BYR PDB file includes the entire spike trimer along with the antibody, so we removed all spike residues except for 334–526 of the RDB bound to the antibody. The bare RBD systems were constructed from the 7BZ5 system by deleting the antibody. The mutations in the RBD were added using CHARMM-GUI input generator.⁶⁶ Each system was solvated in a cuboidal water box with 12 Å of padding in each direction, and Na⁺ and Cl⁻ ions were added to neutralize the total charge and maintain an ionic strength of 150 mM. The hydrogen mass repartitioning (HMR) technique⁶⁷ was used to allow for 4 fs time-steps, facilitating longer simulation runs. The topology files for HMR were prepared using the VMD 1.9.4 alpha 51 software.⁶⁸ All systems were modeled using CHARMM36 force field⁶⁹ (with CMAP correction for proteins) and TIP3P parameters⁷⁰ were used for the water and the ions.

Equilibration and unbiased simulation

Each system was minimized using a conjugate gradient algorithm for 10 000 steps followed by an equilibration in NPT ensemble for 250 ps (with 2 fs time step) with harmonic restraint on the alpha carbons of all the protein residues. A 10 ns NPT equilibration was performed (with 4 fs time-step) afterwards for every system with no restraint applied on any atom. The temperature was controlled at 303.15 K using a Langevin thermostat with damping constant 1 ps⁻¹, and a Langevin piston was used to control the pressure at 1 atm.

Production runs were performed for each RBD–antibody complex for 0.5 μs and for each bare RBD system for 0.3 μs. For all systems, (WT and four variants), the cumulative unbiased simulation time was 6.5 μs. Only the last 0.4 μs of each antibody bound system and 0.25 μs of each bare RBD system were used for further analysis. All simulations were performed using NAMD 2.14 simulation package.⁷¹ The trajectories were saved at an interval of 100 ps.

The root mean squared deviation (RMSD) and the residue-wise root mean squared fluctuations (RMSF) were computed using the gmx rms and gmx rmsf tools, respectively, in the GROMACS 2018.1 software.⁷² RMSD for all systems were computed with respect to the first frame of the long production run. Trajectories were aligned only with respect to the α-carbons of the RBD using VMD software.⁶⁸ RMSD was computed for the RBD, antibody, the antigen binding domain (ABD) of the antibody and the whole system. The antigen binding domain is defined as the one of the two monomers of the BD23 protein which interacts with the RBD and the portions of the two monomers directly interacting with RBD (residue 1–108) for the B38 antibody (Fig. 1 in main text). For computing RMSF, the last 400 ns of each trajectory is split into 4 segments of 100 ns each. Residue wise RMSF was computed for each segment and the mean and standard deviations are reported. The RMSF ratio is then computed for each mutant strains as the ratio of RMSF for each residue and that of the corresponding residue in the WT. Principal component analysis was performed using the PyEMMA 2 package.⁷³ The non-bonding interaction energy (including the electrostatic energy and van der Waals energy) between the antibody and the mutated residues of the RBD for different variants were calculated using the NAMD Energy module⁷⁴ of VMD.⁶⁸ Hydrogen bonds between the RBD and antibodies were also analyzed using VMD. Binding free energy between each RBD–antibody pair was computed using the Molecular Mechanics-Generalized Born Surface Area (MM-GBSA) approach⁷⁵ using the last 400 ns of the equilibrium trajectories. The NAMD2.14 package was used to perform the MM-GBSA calculation.⁷⁶

Mutual information-based cross-correlation and protein graph connectivity network

Linear mutual information (LMI) based cross-correlation matrices were computed for all residue pairs in each system. The theoretical details can be found in ref. 77–79. Although the LMI-based cross-correlation is primarily used to quantify the allosteric effect of a small molecule ligand on the receptor protein,^80,81 it can also be used to decipher the allosteric regulation within a protein or protein complex.⁶¹

The linear mutual information (LMI) metric (I_ij) is calculated as

I_ij = H_i + H_j − H_ij (1)

where H is a Shannon type entropy function:x_i, x_j are the vectors containing the Cartesian coordinates of the C_α atoms, p(x_i) and p(x_j) are their marginal distributions and p(x_i,x_j) is their joint distribution.^79,81 A Pearson like correlation coefficient C_ij is defined from I_ij in the following manner:

C_ij = (1 − e^−(2/d)I_ij)^−1/2 (3)

where d is the dimensionality of the space. Unlike the conventional correlation coefficient, C_ij defined by eqn (3) can only take values between zero and one.

The residue-wise cross-correlation matrices were averaged over the four 100 ns segments prepared from the last 400 ns of each unbiased production run. The changes of the cross-correlation across the mutants were used to decipher the role of the mutations in the large scale dynamics. The uncertainty associated with the cross-correlations computed from the four 100 ns segments was found to be quite low, which allowed for relative comparison between different variants (ESI Fig. S12†). We also constructed a protein graph connectivity network model based on the cross-correlation matrices. The α-carbon of each residue is considered as the nodes of the graph and the cross-correlation between residues as the weight of the edges. Two residues were considered connected only if the heavy atoms of those residues were within 5 Å of each other during at least 75% of the entire production run. The betweenness centrality (BC) of each residue in the network was computed and compared between various mutant strains. The BC quantifies the information flow between nodes or edges of a network. As the α-carbon of each residue is a node in the protein-graph-connectivity network, the value of the BC denotes the functional relevance of each residue in the dynamics of allosteric information flow through the protein chains. If a node, i, functions as a bridge between two other nodes along the shortest path joining those two, the BC of the node i is

where g_st is the total number of geodesics a.k.a. the shortest paths joining nodes s and t, out of which nⁱ_st paths include the node i.^61,80,81 The change in BC, due to RBD mutations, is used to quantify the dynamical role of mutations in modifying the information flow pathway in the RBD-spike complex. The LMI based cross-correlations were calculated and the protein graph connectivity networks were constructed using the bio3D package.⁸²

Steered molecular dynamics simulation

To study the unbinding of the antibodies from the spike RBD, steered molecular dynamics (SMD) simulations⁸³ were performed. For SMD simulation, all systems were modeled with the Generalized Born (GB) implicit solvent model⁸⁴ with solvent accessible surface area (SASA) correction (GBSA). The ion concentration was set at 0.15 M. Instead of applying the pulling force to a specific atom or set of atoms, we defined a collective variable as the center of mass distance between the RBD and the antibody. A fictitious particle was attached to the pulling coordinate with a harmonic spring with a force constant of 1 kcal mol⁻¹ Å⁻¹, to allow for sufficient fluctuations. This fictitious particle is pulled at a constant speed of 4 Å ns⁻¹. The external biasing potential applied to the system at time t can be given by:where k is the force constant of the harmonic spring, r^3N is the 3N dimensional atomic coordinates, s is the collective variable, and v is the pulling speed. Five SMD trajectories of 10 ns length were generated for each system. Each trajectory was saved at an interval of 10 ps, and the pulling force and the accumulated work were recorded every 0.1 ps. The starting structures for the SMD trajectories were sampled from the unbiased production runs at the following points: 300 ns, 350 ns, 400 ns, 450 ns and 500 ns. All SMD simulations were started without further equilibration. The colvars⁸⁵ module in conjunction with NAMD 2.14 (ref. 71) was used for SMD simulations. The SMD input files are provided in the GitHub repository and the SMD trajectories and the colvars output files are provided in the Zenodo repository (see Data availability statement).

We also performed additional SMD simulations in an explicit solvent environment to compare our results. The methodological details are provided in the ESI.†

Data availability

All the trajectory data are available from the Zenodo repository: https://doi.org/10.5281/zenodo.5676393. The codes and analysis scripts are available from the GitHub repository: https://github.com/dhimanray/COVID-variant-antibody.git.

Author contributions

DR designed research with inputs from IA. DR and RNQ performed unbiased molecular dynamics simulations and analyzed results. DR performed SMD simulations and analyzed results. DR constructed the protein-graph connectivity network model using mutual information based cross correlation. DR and IA wrote the paper. All authors contributed to discussing the results and reviewing the manuscript.

Conflicts of interest

The authors declare no competing financial interest.

Acknowledgements

The authors thank Ly Le, Trevor Gokey, Sharon Stone, and Moises Romero for their suggestions and feedback on the manuscript. This work was supported partially by the National Science Foundation (NSF) via grant MCB 2028443. The authors thank the Triton Shared Computing Cluster (TSCC) in San Diego Supercomputer Center (SDSC), and University of California Irvine High Performance Computing (HPC) facility for providing the computational resources.

References

1. P. Zhou, et al., A pneumonia outbreak associated with a new coronavirus of probable bat origin, Nature, 2020, 579(7798), 270–273.
2. COVID Live Update: 254,302,532 Cases and 5,118,809 Deaths from the Coronavirus – Worldometer, https://www.worldometers.info/coronavirus/#countries.
3. C.-L. Hsieh, et al., Structure-based design of prefusion-stabilized SARS-CoV-2 spikes, Science, 2020, 369, 1501–1505.
4. E. Andreano, G. Piccini, D. Licastro, L. Casalino, N. V. Johnson, S. Dal Monego, E. Pantano, N. Manganaro, A. Manenti and R. Manna, SARS-CoV-2 escape from a highly neutralizing COVID-19 convalescent plasma, Proc. Natl. Acad. Sci. U. S. A., 2021, 118(36), e2103154118.
5. S. A. Serapian, F. Marchetti, A. Triveri, G. Morra, M. Meli, E. Moroni, G. A. Sautto, A. Rasola and G. Colombo, The Answer Lies in the Energy: How Simple Atomistic Molecular Dynamics Simulations May Hold the Key to Epitope Prediction on the Fully Glycosylated SARS-CoV-2 Spike Protein, J. Phys. Chem. Lett., 2020, 11, 8084–8093.
6. T. N. Starr, N. Czudnochowski, Z. Liu, F. Zatta, Y. J. Park, A. Addetia, D. Pinto, M. Beltramello, P. Hernandez, A. J. Greaney and R. Marzi, SARS-CoV-2 RBD antibodies that maximize breadth and resistance to escape, Nature, 2021, 597(7874), 97–102.
7. A. Morris, W. McCorkindale, T. C. Moonshot Consortium, N. Drayman, J. D. Chodera, S. Tay, N. London and A. A. Lee, Discovery of SARS-CoV-2 main protease inhibitors using a synthesis-directed de novo design model, Chem. Commun., 2021, 57, 5909–5912.
8. W. Tai, L. He, X. Zhang, J. Pu, D. Voronin, S. Jiang, Y. Zhou and L. Du, Characterization of the receptor-binding domain (RBD) of 2019 novel coronavirus: implication for development of RBD protein as a viral attachment inhibitor and vaccine, Cell. Mol. Immunol., 2020, 17(6), 613–620.
9. X. Ou, et al., Characterization of spike glycoprotein of SARS-CoV-2 on virus entry and its immune cross-reactivity with SARS-CoV, Nat. Commun., 2020, 11(1), 1–12.
10. T. F. Rogers, et al., Isolation of potent SARS-CoV-2 neutralizing antibodies and protection from disease in a small animal model, Science, 2020, 369, 956–963.
11. W. Li, et al., Rapid identification of a human antibody with high prophylactic and therapeutic efficacy in three animal models of SARS-CoV-2 infection, Proc. Natl. Acad. Sci., 2020, 117, 29832–29838.
12. S. J. Zost, et al., Potently neutralizing and protective human antibodies against SARS-CoV-2, Nature, 2020, 584, 443–449.
13. S. Jiang; X. Zhang and L. Du, Therapeutic antibodies and fusion inhibitors targeting the spike protein of SARS-CoV-2, 2020,  DOI:10.1080/14728222.2020.1820482.
14. R. Chowdhury, V. S. S. Adury, A. Vijay, R. K. Singh and A. Mukherjee, Atomistic De-novo Inhibitor Generation-Guided Drug Repurposing for SARS-CoV-2 Spike Protein with Free-Energy Validation by Well-Tempered Metadynamics, Chem.–Asian J., 2021, 16, 1634–1642.
15. M. Smith and J. C. Smith, Repurposing Therapeutics for COVID-19: Supercomputer-Based Docking to the SARS-CoV-2 Viral Spike Protein and Viral Spike Protein-Human ACE2 Interface, 2020.
16. Y. Han and P. Král, Computational Design of ACE2-Based Peptide Inhibitors of SARS-CoV-2, ACS Nano, 2020, 14, 5143–5147.
17. F. P. Polack, et al., Safety and Efficacy of the BNT162b2 mRNA Covid-19 Vaccine, 2020, vol. 383, pp. 2603–2615,  DOI:10.1056/NEJMoa2034577.
18. L. R. Baden, et al., Efficacy and Safety of the mRNA-1273 SARS-CoV-2 Vaccine, 2020, vol. 384, pp. 403–416,  DOI:10.1056/NEJMoa2035389.
19. M. Voysey, et al., Safety and efficacy of the ChAdOx1 nCoV-19 vaccine (AZD1222) against SARS-CoV-2: an interim analysis of four randomised controlled trials in Brazil, South Africa, and the UK, The Lancet, 2021, 397, 99–111.
20. G. Köhler and C. Milstein, Continuous cultures of fused cells secreting antibody of predefined specificity, Nature, 1975, 256(5517), 495–497.
21. I. Correia, Stability of IgG isotypes in serum, mAbs, 2010, 2, 221–232.
22. S. Sun, P. Akkapeddi, M. C. Marques, N. Martínez-Sáez, V. M. Torres, C. Cordeiro, O. Boutureira and G. J. Bernardes, One-pot stapling of interchain disulfides of antibodies using an isobutylene motif, Org. Biomol. Chem., 2019, 17, 2005–2012.
23. A. McAuley, J. Jacob, C. G. Kolvenbach, K. Westland, H. J. Lee, S. R. Brych, D. Rehder, G. R. Kleemann, D. N. Brems and M. Matsumura, Contributions of a disulfide bond to the structure, stability, and dimerization of human IgG1 antibody CH3 domain, Protein Sci., 2008, 17, 95–106.
24. P. J. Carter, Potent antibody therapeutics by design, Nat. Rev. Immunol., 2006, 6(5), 343–357.
25. L. A. Rabia, A. A. Desai, H. S. Jhajj and P. M. Tessier, Understanding and overcoming trade-offs between antibody affinity, specificity, stability and solubility, Biochem. Eng. J., 2018, 137, 365–374.
26. C. W. Davis, et al., Longitudinal Analysis of the Human B Cell Response to Ebola Virus Infection, Cell, 2019, 177, 1566–1582.
27. P. C. Taylor, A. C. Adams, M. M. Hufford, I. d. l. Torre, K. Winthrop and R. L. Gottlieb, Neutralizing monoclonal antibodies for treatment of COVID-19, Nat. Rev. Immunol., 2021, 21(6), 382–393.
28. C. O. Barnes, et al., Structures of Human Antibodies Bound to SARS-CoV-2 Spike Reveal Common Epitopes and Recurrent Features of Antibodies, Cell, 2020, 182, 828–842.
29. L. Liu, et al., Potent neutralizing antibodies against multiple epitopes on SARS-CoV-2 spike, Nature, 2020, 584(7821), 450–456.
30. J. Hansen, et al., Studies in humanized mice and convalescent humans yield a SARS-CoV-2 antibody cocktail, Science, 2020, 369, 1010–1014.
31. Y. Wu, et al., A noncompeting pair of human neutralizing antibodies block COVID-19 virus binding to its receptor ACE2, Science, 2020, 368, 1274–1278.
32. Y. Cao, et al., Potent Neutralizing Antibodies against SARS-CoV-2 Identified by High-Throughput Single-Cell Sequencing of Convalescent Patients' B Cells, Cell, 2020, 182, 73–84.
33. D. Zhou, et al., Structural basis for the neutralization of SARS-CoV-2 by an antibody from a convalescent patient, Nat. Struct. Mol. Biol., 2020, 27(10), 950–958.
34. Z. Lv, et al., Structural basis for neutralization of SARS-CoV-2 and SARS-CoV by a potent therapeutic antibody, Science, 2020, 369, 1505–1509.
35. J. Huo, et al., Neutralizing nanobodies bind SARS-CoV-2 spike RBD and block interaction with ACE2, Nat. Struct. Mol. Biol., 2020, 27(9), 846–854.
36. K. Leung, M. H. Shum, G. M. Leung, T. T. Lam and J. T. Wu, Early transmissibility assessment of the N501Y mutant strains of SARS-CoV-2 in the United Kingdom, October to November 2020, Eurosurveillance, 2021, 26, 2002106.
37. M. S. Graham, et al., Changes in symptomatology, reinfection, and transmissibility associated with the SARS-CoV-2 variant B.1.1.7: an ecological study, The Lancet Public Health, 2021, 6, e335–e345.
38. A. J. Greaney, A. N. Loes, K. H. Crawford, T. N. Starr, K. D. Malone, H. Y. Chu and J. D. Bloom, Comprehensive mapping of mutations in the SARS-CoV-2 receptor-binding domain that affect recognition by polyclonal human plasma antibodies, Cell Host Microbe, 2021, 29, 463–476.
39. Public Health England, SARS-CoV-2 variants of concern and variants under investigation, 2021.
40. B. A. Cohn, P. M. Cirillo, C. C. Murphy, N. Y. Krigbaum and A. W. Wallace, SARS-CoV-2 vaccine protection and deaths among US veterans during 2021, Science, 2021, 375(6578), 331–336.
41. Y. Goldberg, M. Mandel, Y. M. Bar-On, O. Bodenheimer, L. Freedman, E. J. Haas, R. Milo, S. Alroy-Preis, N. Ash and A. Huppert, Waning Immunity after the BNT162b2 Vaccine in Israel, 2021,  DOI:10.1056/NEJMoa2114228.
42. M. Koehler, A. Ray, R. A. Moreira, B. Juniku, A. B. Poma and D. Alsteens, Molecular insights into receptor binding energetics and neutralization of SARS-CoV-2 variants, Nat. Commun., 2021, 12(1), 1–13.
43. M. H. Cheng, J. M. Krieger, B. Kaynak, M. Arditi and I. Bahar, Impact of South African 501.V2 Variant on SARS-CoV-2 Spike Infectivity and Neutralization: A Structure-based Computational Assessment, bioRxiv, 2021 DOI:10.1101/2021.01.10.426143.
44. G. M. Verkhivker, S. Agajanian, D. Y. Oztas and G. Gupta, Comparative Perturbation-Based Modeling of the SARS-CoV-2 Spike Protein Binding with Host Receptor and Neutralizing Antibodies: Structurally Adaptable Allosteric Communication Hotspots Define Spike Sites Targeted by Global Circulating Mutations, Biochemistry, 2021, 60(19), 1459–1484.
45. W. B. Wang, Y. Liang, Y. Q. Jin, J. Zhang, J. G. Su and Q. M. Li, E484K mutation in SARS-CoV-2 RBD enhances binding affinity with hACE2 but reduces interactions with neutralizing antibodies and nanobodies: Binding free energy calculation studies, J. Mol. Graphics Modell., 2021, 109, 108035.
46. A. Spinello, A. Saltalamacchia, J. Borišek and A. Magistrato, Allosteric Cross-Talk among Spike's Receptor-Binding Domain Mutations of the SARS-CoV-2 South African Variant Triggers an Effective Hijacking of Human Cell Receptor, J. Phys. Chem. Lett., 2021, 12, 5987–5993.
47. M. Golcuk, A. Hacisuleyman, B. Erman, A. Yildiz and M. Gur, Binding mechanism of neutralizing Nanobodies targeting SARS-CoV-2 Spike Glycoprotein, J. Chem. Inf. Model., 2021, 61(10), 5152–5160.
48. A. Pavlova, Z. Zhang, A. Acharya, D. L. Lynch, Y. T. Pang, Z. Mou, J. M. Parks, C. Chipot and J. C. Gumbart, Machine Learning Reveals the Critical Interactions for SARS-CoV-2 Spike Protein Binding to ACE2, J. Phys. Chem. Lett., 2021, 12, 5494–5502.
49. M. L. Mugnai, C. Templeton, R. Elber and D. Thirumalai, Role of Long-range Allosteric Communication in Determining the Stability and Disassembly of SARS-COV-2 in Complex with ACE2, bioRxiv, 2020 DOI:10.1101/2020.11.30.405340.
50. M. H. Cheng, J. M. Krieger, A. Banerjee, Y. Xiang, B. Kaynak, Y. Shi, M. Arditi and I. Bahar, Impact of new variants on SARS-CoV-2 infectivity and neutralization: A molecular assessment of the alterations in the spike-host protein interactions, iScience, 2022, 25, 103939.
51. A. Uyar and A. Dickson, Perturbation of ACE2 structural ensembles by SARS-CoV-2 spike protein binding, J. Chem. Theory Comput., 2021, 17(9), 5896–5906.
52. L. Casalino, Z. Gaieb, J. A. Goldsmith, C. K. Hjorth, A. C. Dommer, A. M. Harbison, C. A. Fogarty, E. P. Barros, B. C. Taylor, J. S. Mclellan, E. Fadda and R. E. Amaro, Beyond shielding: The roles of glycans in the SARS-CoV-2 spike protein, ACS Cent. Sci., 2020, 6, 1722–1734.
53. T. Sztain, et al., A glycan gate controls opening of the SARS-CoV-2 spike protein, Nat. Chem., 2021, 13(10), 963–968.
54. K. Nguyen, S. Chakraborty, R. A. Mansbach, B. Korber and S. Gnanakaran, Exploring the Role of Glycans in the Interaction of SARS-CoV-2 RBD and Human Receptor ACE2, Viruses, 2021, 13, 927.
55. A. M. Harbison, C. A. Fogarty, T. K. Phung, A. Satheesan, B. L. Schulz and E. Fadda, Fine-tuning the spike: role of the nature and topology of the glycan shield in the structure and dynamics of the SARS-CoV-2 S, Chem. Sci., 2022, 13, 386–395.
56. Y. Cong, Y. Feng, H. Ni, F. Zhi, Y. Miao, B. Fang, L. Zhang and J. Z. Zhang, Anchor-Locker Binding Mechanism of the Coronavirus Spike Protein to Human ACE2: Insights from Computational Analysis, J. Chem. Inf. Model., 2021, 61, 3529–3542.
57. C. Jarzynski, Equilibrium free-energy differences from nonequilibrium measurements: A master-equation approach, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1997, 56, 5018.
58. G. E. Crooks, Entropy production fluctuation theorem and the nonequilibrium work relation for free energy differences, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1999, 60, 2721.
59. H. A. Pham, D. T. Truong and M. S. Li, Dependence of Work on the Pulling Speed in Mechanical Ligand Unbinding, J. Phys. Chem. B, 2021, 125, 8325–8330.
60. J. Gore, F. Ritort and C. Bustamante, Bias and error in estimates of equilibrium free-energy differences from nonequilibrium measurements, Proc. Natl. Acad. Sci. U. S. A., 2003, 100, 12564–12569.
61. D. Ray, L. Le and I. Andricioaei, Distant residues modulate conformational opening in SARS-CoV-2 spike protein, Proc. Natl. Acad. Sci., 2021, 118(43), e2100943118.
62. A. Dashti, G. Mashayekhi, M. Shekhar, D. Ben Hail, S. Salah, P. Schwander, A. des Georges, A. Singharoy, J. Frank and A. Ourmazd, Retrieving functional pathways of biomolecules from single-particle snapshots, Nat. Commun., 2020, 11(1), 1–14.
63. Y. Toyama and L. E. Kay, Probing allosteric interactions in homo-oligomeric molecular machines using solution NMR spectroscopy, Proc. Natl. Acad. Sci. U. S. A., 2021, 118(50), e2116325118.
64. C. Musselman, Q. Zhang, H. Al-Hashimi and I. Andricioaei, Referencing Strategy for the Direct Comparison of Nuclear Magnetic Resonance and Molecular Dynamics Motional Parameters in RNA, J. Phys. Chem. B, 2009, 114, 929–939.
65. A. Gowrisankar, T. M. C. Priyanka and S. Banerjee, Omicron: a mysterious variant of concern, Eur. Phys. J. Plus, 2022, 137(1), 1–8.
66. S. Jo, T. Kim, V. G. Iyer and W. Im, CHARMM-GUI: A web-based graphical user interface for CHARMM, J. Comput. Chem., 2008, 29, 1859–1865.
67. C. W. Hopkins, S. Le Grand, R. C. Walker and A. E. Roitberg, Long-Time-Step Molecular Dynamics through Hydrogen Mass Repartitioning, J. Chem. Theory Comput., 2015, 11, 1864–1874.
68. W. Humphrey, A. Dalke and K. Schulten, VMD: Visual molecular dynamics, J. Mol. Graphics, 1996, 14, 33–38.
69. J. Huang and A. D. Mackerell, CHARMM36 all-atom additive protein force field: Validation based on comparison to NMR data, J. Comput. Chem., 2013, 34, 2135–2145.
70. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, Comparison of simple potential functions for simulating liquid water, J. Chem. Phys., 1983, 79, 926–935.
71. J. C. Phillips, et al., Scalable molecular dynamics on CPU and GPU architectures with NAMD, J. Chem. Phys., 2020, 153, 044130.
72. M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess and E. Lindah, Gromacs: High performance molecular simulations through multi-level parallelism from laptops to supercomputers, SoftwareX, 2015, 1–2, 19–25.
73. M. K. Scherer, B. Trendelkamp-Schroer, F. Paul, G. Pérez-Hernández, M. Hoffmann, N. Plattner, C. Wehmeyer, J.-H. Prinz and F. Noé, PyEMMA 2: A Software Package for Estimation, Validation, and Analysis of Markov Models, J. Chem. Theory Comput., 2015, 11, 5525–5542.
74. NAMD Energy Plugin, Version 1.4, accessed August 2021, https://www.ks.uiuc.edu/Research/vmd/plugins/namdenergy/.
75. N. Homeyer and H. Gohlke, Free Energy Calculations by the Molecular Mechanics Poisson-Boltzmann Surface Area Method, Mol. Inf., 2012, 31, 114–122.
76. A. Vergara-Jaque, J. Comer, L. Monsalve, F. D. González-Nilo and C. Sandoval, Computationally Efficient Methodology for Atomic-Level Characterization of Dendrimer–Drug Complexes: A Comparison of Amine- and Acetyl-Terminated PAMAM, J. Phys. Chem. B, 2013, 117, 6801–6813.
77. A. Kraskov, H. Stögbauer and P. Grassberger, Estimating mutual information, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2004, 69, 066138.
78. O. F. Lange and H. Grubmüller, Generalized correlation for biomolecular dynamics, Proteins: Struct., Funct., Bioinf., 2006, 62, 1053–1061.
79. S. Bowerman and J. Wereszczynski, Methods in Enzymology, Academic Press Inc., 2016, vol. 578, pp. 429–447.
80. I. Rivalta, M. M. Sultan, N.-S. Lee, G. A. Manley, J. P. Loria and V. S. Batista, Allosteric pathways in imidazole glycerol phosphate synthase, Proc. Natl. Acad. Sci., 2012, 109, E1428–E1436.
81. C. F. A. Negre, U. N. Morzan, H. P. Hendrickson, R. Pal, G. P. Lisi, J. P. Loria, I. Rivalta, J. Ho and V. S. Batista, Eigenvector centrality for characterization of protein allosteric pathways, Proc. Natl. Acad. Sci., 2018, 115, E12201–E12208.
82. B. J. Grant, A. P. C. Rodrigues, K. M. ElSawy, J. A. McCammon and L. S. D. Caves, Bio3d: an R package for the comparative analysis of protein structures, Bioinformatics, 2006, 22, 2695–2696.
83. S. Park, F. Khalili-Araghi, E. Tajkhorshid and K. Schulten, Free energy calculation from steered molecular dynamics simulations using Jarzynski's equality, J. Chem. Phys., 2003, 119, 3559–3566.
84. A. Onufriev, D. Bashford and D. A. Case, Exploring protein native states and large-scale conformational changes with a modified generalized born model, Proteins: Struct., Funct., Bioinf., 2004, 55, 383–394.
85. G. Fiorin, M. L. Klein and J. Hénin, Using collective variables to drive molecular dynamics simulations, Mol. Phys., 2013, 111, 3345–3362.

---

Footnotes

† Electronic supplementary information (ESI) available: Table S1: list of all simulations performed in this study. Protocol and results for the steered MD simulation in explicit solvent. Fig. S1–S27: results of RMSD, RMSF, PCA, mutual information, betweenness centrality, hydrogen-bonding analysis, and interaction energy measurements of the RBD–antibody systems and RBD only systems. Results of steered MD simulation and the antibody dissociation mechanism for all variants. See https://doi.org/10.1039/d2sc00534d

‡ Current address: Atomistic Simulations, Italian Institute of Technology, Via Enrico Melen 83, 16152, Genova, Italy, E-mail: dhiman.ray@iit.it.

---