key: cord-0718921-0eqzn37q
authors: Hendy, Mohamed; Kaufman, Samuel; Ponga, Mauricio
title: Molecular strategies for antibody binding and escape of SARS-CoV-2 and its mutations
date: 2021-11-05
journal: Sci Rep
DOI: 10.1038/s41598-021-01081-0
sha: b98e1b8c12b29e35df52d703e4bfd5177aae0be8
doc_id: 718921
cord_uid: 0eqzn37q

The COVID19 pandemic, caused by SARS-CoV-2, has infected more than 200 million people worldwide. Due to the rapid spreading of SARS-CoV-2 and its impact, it is paramount to find effective treatments against it. Human neutralizing antibodies are an effective method to fight viral infection. However, the recent discovery of new strains that substantially change the S-protein sequence has raised concern about vaccines and antibodies’ effectiveness. Here, using molecular simulations, we investigated the binding mechanisms between the S-protein and several antibodies. Multiple mutations were included to understand the strategies for antibody escape in new variants. We found that the combination of mutations K417N, E484K, L452R, and T478K produced higher binding energy to ACE2 than the wild type, suggesting higher efficiency to enter host cells. The mutations’ effect depends on the antibody class. While Class I enhances the binding avidity in the presence of N501Y mutation, class II antibodies showed a sharp decline in the binding affinity. Our simulations suggest that Class I antibodies will remain effective against the new strains. In contrast, Class II antibodies will have less affinity to the S-protein, potentially affecting these antibodies’ efficiency.

classified Class IV as involving previously described antibodies that do not block the ACE2 and bind only up RBDs. An additional class of hNAbs against SARS-CoV-2 has been identified, particularly antibodies that bind outside the RBD motif in the S-protein. For instance, Chi et al. 25 have shown that a potent antibody links to the N-terminal domain 26 . Other effective antibodies include the broadly neutralizing antibodies (bnAbs) category such as ADG-2 and S2X259 27, 28 .

Despite detailed information about these antibodies' structural and biophysical aspects, essential questions remain unanswered. For instance, the discovery of rapidly spreading new SARS-CoV-2 lineages with mutations in the S-protein, particularly in its RBD, has raised concerns among the health authorities and scientific community. These variants of concern have motivated the experimental study of the effect of mutations in the S protein on the binding energy of ACE2 receptors and in hNAbs. Even though deep mutational scanning can be used to investigate the effect of the mutations on antibody response, the molecular mechanisms for antibody binding and escape cannot be understand using these techniques 17, 21, [29] [30] [31] [32] . Moreover, it is also unknown if the new strains propagate faster due to an easier attachment to the receptors or other biological factors. Finally, the different hNAbs classes' detailed binding mechanisms remain unavailable. This information is vital because it can enable researchers to understand what type of antibodies are more efficacious to bind to the S-protein and aid effective therapeutics with antibodies against COVID19. Molecular simulations can be used to answer some of these questions and identify key molecular binding mechanisms between SARS-CoV-2's new strains and hNAbs. This information can be used to understand how viral mutations impact vaccine developed immunity and to generate treatments that are robust against virus evolution 17, 29 .

This work investigates the mechanisms for binding and strategies to escape several hNAbs used by SARS-CoV-2 strains using molecular simulations. We pursued the same approach as described by Ponga 6 to compute binding energies and summarized in "Methods". Remarkably, this approach predicted the binding affinity and detachment force between the S-protein and ACE2 receptor in good agreement with experimental measures 1,3,5 . We benchmarked five hNAbs (two in Classes I and II and one in Class III) with the RBD, including three mutations N501Y (MT1), E484K, and N501Y (MT2), and K417N, E484K, and N501Y (MT3). These mutations resemble similarities with the novel strains B1.1.7 33 37 , were also simulated in the RBD to compute the binding energy with the ACE2. A summary of the simulations can be found in Table 1 .

We first focused on understanding the effect of the mutations on the binding to the human ACE2 receptor. Figure 1a shows the binding energy for the WT, as well as for the RBD with mutations MT1, MT2, MT3, and δ . We observed that the WT's binding energy is in the neighborhood of 70.6 ± 3.7 kJ mol −1 . At the same time, MT1 has slightly smaller binding energy values, e.g., 60 ± 3.3 . On the other hand, binding energy values for MT1, MT2 and MT3 of 72 ± 3 kJ mol −1 , 73.3 ± 3 kJ mol −1 , and 83.0 ± 4 kJ mol −1 were obtained, respectively. δ had a binding energy near to 78 ± 3.8 kJ mol −1 . We also noted that the binding energy between WT RBD and ACE2 receptor is close to reported values reported in literature. For instance, using MD simulations, Ortega et al. 38 reported ∼ 15.7 kCal mol −1 ( ∼ 66 kJ mol −1 ), and other works reported ∼ 77 kJ mol −1 when the two proteins were fully detached 6 . The slight difference could be due to differences in the geometry of the RBD and the model used for interatomic interactions. Despite these small differences, the values found in this work for the WT RDB and ACE2 are in close agreement to other studies 16 . The PMF evolution between RBD and ACE2 for WT, MT1-3 and δ can be found in Fig. S1a .

In order to obtain a better understanding of the effect of the mutation in the RBD, we computed the root mean square fluctuation (RMSF) for the WT and mutated versions for the equilibrium configuration ( = 0 nm), and when the proteins were completely detached ( = 3 nm). The RMSF is shown in Fig. 1b ,c for = 0 and 3 nm, respectively. The shaded zones represent residues in contact with ACE2. We observed that there were only Table 1 . Details of the simulations used in this work, including the antibody name and PDB file used in the simulations. Total number of atoms in the simulation includes the proteins, and solvent used in the different simulations. Number of Cl − ions used to balance the charge in the simulations. The Class of the antibodies is also specified. The symbols indicate whether the antibody links to the RBD in the up ( ↑ ), up/down ( ↑↓ ), or up/ down from the side ( ↑↓← ). The last three columns indicate if the antibody was modeled with mutations MT1: N501Y, MT2: E484K, N501Y, and MT3: K417N, E484K, N501Y. For other residues, the RMSF remained lower around 0.1 nm. The agreement in the RMSF is remarkable, despite having used different models. The higher binding energy between RBD containing MT2 and MT3 can be explained by analyzing the residues in contact. Besides having fewer contacts, MT2, MT3 and δ have a considerably larger binding energy than WT and MT1. We observed the introduction of K417N eliminates the salt bridge between RBD (K417) and ACE2 (D30). The introduction of E484K generated another one between K484 and E35 in ACE2. When looked carefully, we observed that the distances between the new bridge are much smaller ( ∼ 0.13 nm for MT3 compared with ∼ 0.3 nm in the MT1, see Tables S1-S4). The significantly lower distance makes stronger bonds, explaining the higher binding affinity between MT2-3 and ACE2.

Binding energy between RBD and hNAbs. Next, we investigated the binding energy of five hNAbs described above with the SARS-CoV-2 RBD. The WT S-protein and the three mutations MT1, 2, and 3 were included in the analysis. Let us first describe the residues in contact between the WT S-protein and the hNAbs shown in Fig. 2a . The red letters correspond to the amino acid residues in contact with the different hNAbs. The underlined letters show the positions of the mutation. Figure 2b shows the location of these residues in the RBD for each antibody with the locations of the selected mutations placed on the RBD. We observed that the mutations were placed on strategic locations of the S-sequence, with the only unaffected binding site unaltered for C135. In particular, K417 generated interactions with three antibodies (B38, C102, C144), E484 interacted with two antibodies in Class II (P2B-2F6 and C144). In contrast, N501 interacted with Class I hNAbs (B38 and C102). The analysis allowed us to strategically make simulations with the mutated RBD with antibodies as detailed in Table 1 . www.nature.com/scientificreports/ Let us now focus on the analysis of the binding energy between Class I hNAbs and the RBD. As shown in Fig. 2c , the binding energy between the WT and B38 increased to ∼ 66 ± 3.3 kJ mol −1 when all contacts were lost at approximately = 1.5 nm (see Fig. S1b ). For MT1 and MT3 the last contacts were lost around = 2.0 nm reaching a maximum PMF value of ∼ 84 ± 4.4 , and ∼ 114 ± 4.4 kJ · mol −1 , respectively. The significant difference in the binding energy between mutated and original RBD is explained in the next section's residues interactions.

Let us now analyze the binding energy for Class I hNAb C102, shown in Fig. 2c . Remarkably, the binding energy between the WT RBD was in the neighborhood of ∼ 62.7 ± 3.55 kJ mol −1 , while for MT1 and MT3 it reached ∼ 85.6 ± 3.8 , and ∼ 63.4 ± 3.8 kJ mol −1 , respectively (see Fig 2b for position of the mutations). The similar binding energy between ACE2 and C102 indicates a similar dissociation constant, illustrating that both cell receptor and antibody have similar probability to link to the S-protein. Like B38, MT1 increased the binding energy for C102. However, the combination in MT3 significantly reduced the binding energy to levels below those found for the WT and ACE2. Thus, while C102 remained active under MT1 (e.g., B.1.1.7 strain), it might lose efficiency against MT3 (e.g., B.1.351 strain). An evolution of the PMF for C102 is shown in Fig. S1c .

We now proceed to quantify the binding energy for Class II hNAbs, in particular, for P2B-2F6 and C144 hNAbs. The binding energies and PMF are shown in Fig. 2d and Fig. S1d ,e, respectively. The binding energy between the WT S-protein and P2B-2F6 was measured to be ∼ 77.7 ± 2.4 kJ mol −1 . The significantly higher binding energy compared with the ACE2 receptor indicates a strong binding configuration-as denoted in the PMF plot shown in Fig. S1d -and illustrates why P2B-2F6 is a potent hNAb. On the other hand, we observed that MT2-which includes mutation E484K-had a deleterious effect on the binding energy for P2B-2F6, reducing it to ∼ 51 ± 2.3 kJ mol −1 . This reduction of the binding energy −due to removing a salt bridge between RBD and P2B-2F6 (see next section)-rendered the antibody much less effective to bind to the S-protein. The 'x' symbol indicates that the binding energy was not computed for that combination since the mutation did not affect the binding residues. www.nature.com/scientificreports/ For C144, another Class II hNAb, we found a similar behavior. The binding energy between WT RBD and C144 was ∼ 76.1 ± 5.25 kJ mol −1 , and the PMF (see Fig. 2d; Fig. S1e ) was characterized by a strong binding configuration. The similar binding level between WT and Class II hNAbs was remarkable. It illustrates why these types of antibodies were effective against WT SARS-CoV-2. However, the impact of the mutations was drastic for C144. We observed that both MT2 and MT3 produced harmful effects on the binding of the antibody and S-protein. The binding energy was reduced to minimum levels of 41 ± 3.9 and 33.5 ± 3.9 kJ mol −1 for MT2 and MT3, respectively. Unfortunately, the severe reduction in the binding energy suggests that C144 will not be effective against the S-protein's mutated version (e.g., against the B.1.351 strain).

Next, we analyzed the behavior for C135 hNAb that attaches outside the RBD motif. The total binding energy for C135 is close to ∼ 80.5 ± 5 kJ mol −1 . In particular, the simulation setup was such that the S-protein slid over C135, generating a shearing mechanism. As a result, we found a well-defined energy indicating a stable binding configuration. The PMF increased almost monotonically with until the total detachment was performed (see Fig. S1f ). This notable binding energy illustrates the power of Class III hNAbs to block the S-protein. Since the residues in contact between RBD and C135 did not involve any of the mutated versions, we omitted the analysis of MT1-3 (see Fig 2b) . Figure 2e shows a binding energy map summarizing all binding energies calculated in this work. Cells with 'x' indicated that the particular combination of RBD and antibody was not simulated since the mutated residues did not interact with the antibody. The color scale in the map is the binding energy difference, taking the RBD/ACE2 receptor as a reference. The actual values are shown in each cell, with white fonts for those combinations where the binding energy is below the reference and with black fonts for those combinations above the reference.

It is now beneficial to compare the molecular mechanisms for binding between the S-protein and hNAbs. Figure 3 shows the antibodies Fabs linked to the RBD, with several insets showing the residues in contact due to Hydrogen (H) bonds and salt bridges. The contacts are shown for locations near the mutation points discussed above. Table S5 shows the amino acid residues in contact between the hNAb B38 and the RBD for the equilibrium configuration ( = 0.1 nm), including a detailed description between the residues in the heavy chain (HC) and light chain (LC), as well as the ones in the RBD with the distances between these residues. For www.nature.com/scientificreports/ = 0.02 − 0.13 nm, we observed that all twenty-five residue interactions were present via H-bonds without salt bridges. Thirteen of these interactions were made by the HC, while the LC made twelve. We observed that K417 interacts with the HC and LC, with residues Y52 and N92. Residue N501 interacted with residue S30 in the LC with a distance of ∼ 0.27 nm (see inset in Fig. 3a) . As the reaction coordinate increased, the number of residues in contact diminished as shown in Tables S6-S8 . Before complete detachment, at = 1.3 nm, residue N487 in the RBD was in contact with S31 in the HC. The analysis indicates a stronger interaction of the RBD with the HC for B38.

Besides using the same geometry for MT1, we observed (see Table S9 -S12) that the number of contacts was lower for the mutated RBD with B38. However, the HC increased the number of contacts (sixteen), while the LC decreased them (five). It is this shift in the contacts that was responsible for the higher PMF. Of remarkable interest, residue Y501 did not play a role in the contact with the B38 in the equilibrium configuration (see the bottom inset in Fig. 3a) , but it did interact with S30 in the LC when = 0.2 − 0.5 nm (see Table S11 ). For the MT3 mutations, we observed that the mutated residues did not interact in the equilibrium position between B38 and the RBD (see Table S13 ).

The interacting residues for the hNAb C102 are detailed in Table S14 for the equilibrium configuration ( ∼ 0.1 nm) and shown in Fig. 3b . The similarities between B38 and C102 are quite remarkable. Out of the fourteen residues in the RBD contacting C102, eleven were shared with B38, and all interactions were due to H-bonds (see Table S14 and insets in Fig. 3b) . The lesser number of contacts found in C102 explains the lower binding energy for C102 compared with B38. A key difference is that as the reaction coordinate increased, the LC lost contact with the RBD (see Table S15 ). Of tremendous interest, one of the last residues in contact was between N487 in the RBD and with S53 in the heavy chain, which is very similar to the last contact for B38 (see Table S16 ). This observation suggests that the engagement between RBD and the hNAbs could have started in residue N487 of the RBD.

The contacts between MT1 RBD and the hNAb C102 are detailed in Table S17 -S20 for several values of . The HC of C102 was the main responsible for the contacts (ten out of thirteen). However, the LC's contacts were maintained until ∼ 0.55 nm, justifying the higher PMF observed in Fig. S1c in comparison with the original RBD. In the analysis, we did not find evidence that Y501 interacted with the antibody, unlike the original RBD. For MT3, we observed that residue N417 interacted with the HC. However, the other mutations did not interact with the antibody (see Table S21 ).

Let us now focus on hNAbs Class II, namely P2B-2F6 and C144, whose interacting residues with the RBD are detailed in Tables S22-S24 and shown in Fig. 3c ,d. We observed a different binding strategy than Class I hNAbs, with fewer contacts but with a much stronger binding energy. For the hNAb P2B-2F6, we identified five residues interacting in the RBD with seven residues in the antibody (five in the HC-Y27, S31, S30, H54, and R112-and two--Y34 and N33-in the LC). A key difference in the binding strategy is that in addition to the H-bonds, we observed a salt bridge between R112 and E484 (see inset in Fig. 3c and Table S22 ). As increased, the number of residues decreased. Y43 and R112 in the P2B-2F6 and E484 were the last bonds to be broken, indicating the large interaction forces between these two residues.

For C144, a similar binding strategy was observed. Seven residues in the RBD were in contact with eight in the HC hNAb. The LC did not interact with the RBD since the HC makes direct physical contact with the RBD (see Fig. 3d ). Another salt bridge between E484 and R100 in the HC was observed (inset in Fig. 3d ). The simulation of Class II antibodies illustrated the role of salt bridges and how they can serve to initiate the binding between the two particles.

Mutations MT2-3, involving a change in E484K, explain this mutation's deleterious effect in the binding between RBD and P2B-2F6 and C144. The direct removal of the salt bridge reduced the binding energy significantly for both cases. While K484 interacted with the HC, it was through a weak H-bond with S55 (r = 0.26 nm). Also, nine residues in the RBD were in contact compared to eleven in the WT, explaining the lower binding energy.

Class III antibodies. We now describe the bonding between C135 Class III hNAbs and the RBD as shown in Fig. 4 . Figure 2e shows the remarkable differences between Class III hNAbs with the other classes. We observe that Class III targets residues T345 and residues between N437-Y451 as illustrated in Fig. 2e . We also observed no overlap between residues in Class III and Class I-II, suggesting that combined strategies for neutralization are possible. C135 uses seven residues (five in the HC and two in the LC) to bind to five residues in the RBD (see Table S25 and inset in Fig. 4 ). All these bonds were H-bonds, without any salt bridges. However, since the mechanism for detaching is due to sliding on the surfaces of the RBD, it generated long-lasting contacts (and even new ones-e.g., G97 with T345 and R346 with E50− as shown in Table S25 ). Thus, as the two molecules were pulled apart, the sliding mechanism led to increased binding energy compared to other hNAbs classes.

We have investigated the mechanisms and strategies for antibody binding and escape in SARS-CoV-2. On the one hand, analysis of the WT and MT1 of the RBD and the ACE2 receptor indicated that the binding energy remains at similar levels and between our simulations' margin error. On the other hand, for MT2, MT3, and δ RBD versions involving mutations K417N, L452R, E484K, T478K, and N501Y, the binding energy increased compared to the WT. This observation suggests that the strains B.1.351 36 and B.1.617.2 37 could have a better binding avidity with the ACE2 receptor and better effectivity to enter host cells than the WT virus. This observation seems to agree with recent works where the effect of the mutations have been analyzed. In particular, several experimental works using mutated RBD have shown a more pronounced effect of the mutation E484 17 16, 17 reported that the binding between SARS2 and ACE2 increases when the following mutations K417N, E484K, L452R, T478K, and N501Y are present, in good agreement with our findings. We also observed that the mutations increased the RMSF near the mutation zone, suggesting more motion of the residues near the mutation. The larger values of the RMSF indicate contacts that are on and off at the interface in the mutated versions. However, the biological effects on the overall higher RMSF near residue #501 are unclear since larger fluctuations leads to more exposure of the residues. The binding mechanisms between RBD and hNAbs are vast, and various classes of antibodies use different strategies. Based on our study, Class II and III hNAbs were the most effective to bind to the WT RBD. Considering that Class II and III link to different RBD areas, therapeutics strategies combining Class II and III could be very effective to stop viral infection as suggested by others 18 . Unfortunately, we observed that Class II hNAbs show a sharp decline in the binding energy to MT2 and MT3 and are due to two main factors. First, the loss of the salt bridge made by residue E484 and Class II hNAbs that use gene VH3-53 as a germline (and the same strategy has been seen for P2B-2F6 that is different gene). Second, the reduced number of contacts between Class II and RBD (see Tables S21-S23 ) makes the binding energy decay to low levels. This renders Class II antibodies inefficient to bind to the RBD. This observation suggests that Class II hNAbs are less effective in binding to the S-proteins within the B.1.351 strain 17, 30, 42, 43, 45, 46 . The results are in agreement with experimental observations. For instance, recent experimental work on mutated RBD with E484K has shown a 100-1000 fold increase in antibody dose for C144 17, 30 . This effect is more drastic in C144 since it uses the HC to make contact with the RBD, while P2B-2F6 uses both the HC and LC (see Fig. 3c,d) . A key contribution of our work is that, in addition to be in close agreement with the aforementioned experimental works, it provides a molecular explanation for the mechanisms of binding between RBD and hNAbs, that is not available through current experimental works 17, 30, 42, 43, [47] [48] [49] [50] .

Simulations between the RBD and Class I hNAbs indicated that the binding energy increases between MT1 (N501Y) and B38 and C102 compared to WT. This observation suggests that Class I hNAbs will remain effective to bind to the S-protein for the B.1.1.7 strain. This result is supported by experimental works where N501Y mutation has been found to have a modest effect on the binding of monoclonal antibodies 16, 30, 51 . The binding energy between B38 and C102 had different behavior when analyzed with MT3. B38 had shown higher binding energy 114 kJ mol −1 . On the other hand, C102 showed binding energy significantly below the value found for the ACE2 receptor and MT3, indicating less effectivity to bind to the mutated S-protein. Once again, this observation suggests that C102 will be less effective in blocking viral infection for viruses with the B.1.351 lineage.

Our all-atom MD simulations also have several limitations. For instance, to reduce the computational cost, only part of the S-protein was modeled in this work. Previous studies have shown minor differences in the PMF plot for one single RBD and the full trimeric S-protein 6 . In our study, we have also ignored glycans that attach to the S-protein and can help in the binding to receptors/hNAbs. For instance, Watanabe et al. 52 reported a highly glycosylated S-protein with 22 glycans. They reported that the vast majority of these glycands were mainly placed outside the RBD, with the exception of two glycans located at residues 331 and 343, and this has been used to investigate their role in shielding and modulating the conformational dynamics of the RBD 53 . As shown in Fig. 2a , the only hNAbs that is close to these residues, is C135, while others are far way from these residues. Thus, it is unlikely that the glycans have a large impact on the PMF reported in this study, as the interactions between molecules are short range. Additionally, the mutated RBD has been modeled, but cryoEM mutated molecular structures can differ from those used in this work. As far as the authors are concerned, there is not cryoEM of mutated RBD available yet. Differences in the geometry could affect the equilibrium position and, thus, change the energy path. While these limitations can change the actual values found in this work, the mechanisms for molecular binding will remain mostly unaffected and our conclusions will be unaffected.

In conclusion, our results indicate that the landscape for molecular bonding and binding between antibodies and RBD is vast, and several mechanisms are available. Of tremendous interest, we found that Class III hNAbs −that link in the up and down states outside the RBD motif− have the largest binding energy to the wild type (WT) RBD ( ∼ 80 kJ mol −1 ). Class II antibodies, that link in the up and down states in the WT RBD motif, followed ( ∼ 75 − 80 kJ mol −1 ). Class I hNAbs produced slightly smaller binding energies ( ∼ 60 − 70 kJ mol −1 ). www.nature.com/scientificreports/ These binding energy levels are well separated between different classes and all hNAbs produced similar or higher binding energies than the RBD/ACE2 receptor ( 66.9 ± 3.7 kJ mol −1 ) 6 . We also investigated the effect of mutated S-proteins. We found that mutation N501Y marginally enhanced the binding energy between RBD and ACE2, while mutations MT2 (E484K/N501Y) MT3 (K417N/E484K/ N501Y), and δ (L452R/T478K) enhanced the binding avidity between RBD and ACE2. This result could help to partially explain the rapid propagation of the B.1.351, B.1.617.2 and δ strains in different parts of the world. We also investigated the mutations' effect on the binding affinity to several hNAbs, including B38, C102 (Class I), P2B-2F6, and C144 (Class II). The impact of the mutations is antibody-dependent. However, in general, we found less affinity between mutated S-protein and these antibodies, except for B38. Of particular interest, Class II antibodies showed a sharp decrease in their binding energy when mutations MT2-3 were simulated. This result indicates a much lower efficiency to stop the viral infection. Future directions of this work include the computational modeling of hNAbs with better binding affinity to the mutated RBD using data-based approaches and the estimation of antibody dose using mechanistic models.

Molecular dynamics simulations, umbrella sampling and characterization. Molecular dynamics simulations were performed with GROMACS 54-56 simulation package. Molecular geometries were obtained from different protein database (PDB) files, including the 6LZG (RBD and ACE2) 57 , 7BWJ (P2B-2F6) 22 , 7BZ5 (B38) 20 , 7K8M (C102), 7K8Z (C135), 7K90 (C144) 18, 19 . For the antibodies given in files 7BWJ, 7BZ5, 7K8M, the database's molecular geometry includes the SARS-CoV-2 geometry of the RBD in solution with the antibodies. For 7K8Z, 7K90, the PDB files included the full trimeric proteins of the SARS-CoV-2. To make fair comparisons and reduce the simulation runtime, we used a portion of the S-glycoprotein, limited to the RBD between residues 337−505 (NCBI Reference: YP _ 009724390.1). Since some of the PDB files' sequences were incomplete, we used an RBD from 7K8M and aligned the proteins using Pymol. The use of the same RBD ensures the same geometry for all antibodies with a root mean square displacement (RMSD) per atom less than 0.075 nm.

The antibody geometries with the RBD were processed in GROMACS, where the proteins' principal axis was aligned along the x−direction, which was then used as the main reaction coordinate for the steered MD simulation. Some of the protein files had missing atoms that were not resolved in the cryoEM identification. We included these atoms using the SwissViewer. Once the proteins were complete, we solvated them in water using the TIP3P model to achieve a density of approximately ρ = 1000 kg m −3 . In order to allow for sufficient space for the steered simulations, we generated computational cells with more than 1 nm between the proteins and end of the cells and sufficient space on top to perform the pulling simulations. The computational cell had dimensions of ∼ 20 × 7.55 × 7.25 (nm). After adding the solvent, the system had a non-zero charge, and ions (Cl − ) were added as needed to equilibrate in all samples. All interatomic forces were computed with the CHARMM force-field 58 . The computational cells had around ∼106,000 atoms, including proteins and solvent. The details of the antibody simulations is shown in Table 1 .

After solvating the system with water, we minimized the system's energy using a non-linear conjugate gradient, with a force tolerance of 1000 kJ mol −1 nm −1 . The relaxed cells were then subjected to a constant temperature ( T = 310.15 K), and pressure ( P = 1 bar) run for about t = 2 ns using a Berendsen thermostat 59 . A timestep of t = 2 fs was used for all simulations. Short-range interactions were treated with a smooth force-switch cutoff of r = 1.2 nm, and long-range electrostatics were treated using the Particle-Mesh-Ewald (PME) formalism, implemented in GROMACS 60 for all simulations. Hydrogen−bonds were restrained with the LINCS algorithm 61 .

To compute the binding affinity between the RBD of the S-proteins and several antibodies, we used a combination of steered MD simulations with umbrella sampling. In all simulations, the antibody (Heavy and Light) chains were fixed in the simulation cell using a position restrain, while the RBD was subject to a pulling force acting on the x− direction. This force was applied using a spring constant of K = 1000 kJ mol −1 nm −2 . The force was applied between the center of mass between the RBD and the antibody. The proteins were then pulled apart from each other at a constant speed of v x = 1 nm ns −1 along the x−direction. The systematic use of the same pulling direction between RBD and hNAbs ensures a direct comparison between different systems, although other paths might yield another PMF evolution.

Configurations were obtained every 20 ps while the PMF was sampled every 0.1 ps, resulting in a total of 10,000 sampling points during the whole simulation. We selected configurations separated by a distance of = 0.05 nm or less in most cases. A convergence test was also performed using 20,000 sampling points for the RBD and ACE2, and for RBD and B38 with WT, MT1, and MT3. The difference in the PMF was well below 3% for all cases justifying the use of less points for other simulations. These configurations were then relaxed using a short NPT ensemble run of 0.1 ns for a posterior sampling of the PMF using the umbrella technique. Umbrella simulations were performed using NPT ensembles for about 10 ns, using a spring constant of K = 1000 kJ mol −1 nm −2 between the center of mass of the molecules. The average maximum displacement from the center of mass for each sampling point was around 0.1 nm, this ensures that the maximum energy applied with the restrain is about 5 kJ mol −1 ( ∼ 2k B T ). This value has been used for computing the binding energy between RBD and ACE2 receptor in other works shielding optimal results 6 . The PMF was computed with the Weighted Histogram Analysis Method (WHAM) 62 using a bootstrap analysis to estimate the uncertainty in the PMF. We used 200 different measures, using 200 binning spaces along with the reaction coordinate and a tolerance of 10 −8 . Instead of using the moment when the two molecules detach physically as Ponga 6 , the binding energy was taken as the value when the PMF plateaus, taking the PMF's average during the last 0.5 nm of the reaction coordinate.

Analysis of the atomistic configurations was performed with the GROMACS cluster analysis tool. We scanned the configurations with a root mean squared displacement between the range of 0.15 − 0.25 nm 63 . The cluster analysis yielded between three to six clusters for the analyzed configurations. In all cases shown, the most 

Structure, function, and antigenicity of the SARS-CoV-2 spike glycoprotein

A novel coronavirus from patients with pneumonia in China

Structure of the SARS-CoV-2 spike receptor-binding domain bound to the ACE2 receptor

Cryo-EM structure of the 2019-nCoV spike in the prefusion conformation

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

Quantifying the adhesive strength between the SARS-CoV-2 S-proteins and human receptor and its effect in therapeutics

The role of furin cleavage site in sars-cov-2 spike protein-mediated membrane fusion in the presence or absence of trypsin

Inhibition of sars-cov-2 (previously 2019-ncov) infection by a highly potent pan-coronavirus fusion inhibitor targeting its spike protein that harbors a high capacity to mediate membrane fusion

Distinct conformational states of SARS-CoV-2 spike protein

Characterization of spike glycoprotein of SARS-CoV-2 on virus entry and its immune cross-reactivity with SARS-CoV

SARS-CoV-2 cell entry depends on ACE2 and TMPRSS2 and is blocked by a clinically proven protease inhibitor

The mechanistic overview of SARS-CoV-2 using angiotensin-converting enzyme 2 to enter the cell for replication: Possible treatment options related to the renin-angiotensin system

Coronavirus membrane fusion mechanism offers a potential target for antiviral development

Battle at the entrance gate: CIITA as a weapon to prevent the internalization of SARS-CoV-2 and Ebola viruses

Enhancing host cell infection by SARS-CoV-2

Prospective mapping of viral mutations that escape antibodies used to treat covid-19

Comprehensive mapping of mutations in the sars-cov-2 receptor-binding domain that affect recognition by polyclonal human plasma antibodies

SARS-CoV-2 neutralizing antibody structures inform therapeutic strategies

Structures of human antibodies bound to SARS-CoV-2 spike reveal common epitopes and recurrent features of antibodies

A noncompeting pair of human neutralizing antibodies block COVID-19 virus binding to its receptor ACE2

Prospective mapping of viral mutations that escape antibodies used to treat covid-19

Human neutralizing antibodies elicited by SARS-CoV-2 infection

Cross-neutralization of SARS-CoV-2 by a human monoclonal SARS-CoV antibody

Studies in humanized mice and convalescent humans yield a sars-cov-2 antibody cocktail

A neutralizing human antibody binds to the N-terminal domain of the Spike protein of SARS-CoV-2

N-terminal domain antigenic mapping reveals a site of vulnerability for sars-cov-2

Broad and potent activity against sars-like viruses by an engineered human monoclonal antibody

Structural basis for broad sarbecovirus neutralization by a human monoclonal antibody

Complete mapping of mutations to the sars-cov-2 spike receptor-binding domain that escape antibody recognition

Evolution of antibody immunity to sars-cov-2

Deep mutational scanning of sars-cov-2 receptor binding domain reveals constraints on folding and ace2 binding

Sars-cov-2 rbd antibodies that maximize breadth and resistance to escape

Recurrent emergence and transmission of a sars-cov-2 spike deletion h69/v70. bioRxiv

Transmission of sars-cov-2 lineage b.1.1.7 in England: Insights from linking epidemiological and genetic data

Genomic characterization of a novel sars-cov-2 lineage from Rio de Janeiro

Emergence and rapid spread of a new severe acute respiratory syndrome-related coronavirus 2 (sars-cov-2) lineage with multiple spike mutations in South Africa

Sars-cov-2 b.1.617.2 delta variant replication and immune evasion

Role of changes in sars-cov-2 spike protein in the interaction with the human ace2 receptor: An in silico analysis

Dynamics of the ace2-sars-cov-2/sars-cov spike protein interface reveal unique mechanisms

Sars-cov-2 escape in vitro from a highly neutralizing covid-19 convalescent plasma

Potent neutralizing antibodies against multiple epitopes on sars-cov-2 spike

Escape from neutralizing antibodies by sars-cov-2 spike protein variants

Escape of sars-cov-2 501y.v2 variants from neutralization by convalescent plasma. medRxiv

Resistance of sars-cov-2 variants to neutralization by monoclonal and serum-derived polyclonal antibodies

Antibody evasion by the p. 1 strain of sars-cov-2

Increased resistance of sars-cov-2 variant p. 1 to antibody neutralization

Antibody cocktail to sars-cov-2 spike protein prevents rapid mutational escape seen with individual antibodies

The monoclonal antibody combination regen-cov protects against sars-cov-2 mutational escape in preclinical and human studies

Tackling covid-19 with neutralizing monoclonal antibodies

Genetic and structural basis for recognition of sars-cov-2 spike protein by a two-antibody cocktail

Antibody resistance of sars-cov-2 variants b. 1.351 and b. 1.1. 7

Site-specific glycan analysis of the sars-cov-2 spike

Beyond shielding: The roles of glycans in the sars-cov-2 spike protein

GROMACS: A message-passing parallel molecular dynamics implementation

0: A package for molecular simulation and trajectory analysis

GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers

Structural and functional basis of sars-cov-2 entry by using human ace2

CHARMM36 all-atom additive protein force field: validation based on comparison to NMR data

Molecular dynamics with coupling to an external bath

A smooth particle mesh Ewald method

A parallel linear constraint solver for molecular simulation

& van der Spoel, D. g_wham-a free weighted histogram analysis implementation including robust error and autocorrelation estimates

Peptide folding: when simulation meets experiment

Esbri: A web server for evaluating salt bridges in proteins

We gratefully acknowledge the support from the Natural Sciences and Engineering Research Council of Canada 

M.P. designed the research plan, performed the simulations and analysis. M.H. and S.K. performed the analysis of the residues, H-bonds, and salt bridges. All authors wrote and made comments to the manuscript.

The authors declare no competing interests.

The online version contains supplementary material available at https:// doi. org/ 10. 1038/ s41598-021-01081-0.Correspondence and requests for materials should be addressed to M.P.Reprints and permissions information is available at www.nature.com/reprints.Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.