key: cord-0888843-sldhdptw
authors: Cao, Yiwei; Choi, Yeol Kyo; Frank, Martin; Woo, Hyeonuk; Park, Sang-Jun; Yeom, Min Sun; Seok, Chaok; Im, Wonpil
title: Dynamic Interactions of Fully Glycosylated SARS-CoV-2 Spike Protein with Various Antibodies
date: 2021-09-16
journal: J Chem Theory Comput
DOI: 10.1021/acs.jctc.1c00552
sha: 13244f682899f255807917439140122985dd7c1e
doc_id: 888843
cord_uid: sldhdptw

[Image: see text] The spread of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) presents a public health crisis, and the vaccines that can induce highly potent neutralizing antibodies are essential for ending the pandemic. The spike (S) protein on the viral envelope mediates human angiotensin-converting enzyme 2 binding and thus is the target of a variety of neutralizing antibodies. In this work, we built various S trimer–antibody complex structures on the basis of the fully glycosylated S protein models described in our previous work and performed all-atom molecular dynamics simulations to gain insight into the structural dynamics and interactions between S protein and antibodies. Investigation of the residues critical for S–antibody binding allows us to predict the potential influence of mutations in SARS-CoV-2 variants. Comparison of the glycan conformations between S-only and S–antibody systems reveals the roles of glycans in S–antibody binding. In addition, we explored the antibody binding modes and the influences of antibody on the motion of S protein receptor binding domains. Overall, our analyses provide a better understanding of S–antibody interactions, and the simulation-based S–antibody interaction maps could be used to predict the influences of S mutation on S–antibody interactions, which will be useful for the development of vaccine and antibody-based therapy.

Since the outbreak of coronavirus disease 2019 (COVID- 19) , the spread of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has infected over 226 million people and led to 4.6 million deaths as of September, 2021. Because of the extremely high contagiousness of SARS-CoV-2 and lack of effective antiviral medicines, neutralizing antibodies (NAbs) acquired through active or passive immunization play critical roles in preventing healthy individuals from infection and accelerating recovery of infected persons. Several SARS-CoV-2 vaccines have been authorized for emergency use, and the preliminary data show that they provide a high protection rate. 1−3 SARS-CoV-2 is an enveloped virus with a positive-sense single-stranded RNA genome. 4 The spike (S) trimer anchored in the viral envelope is a glycoprotein mediating the binding to human angiotensin-converting enzyme 2 (ACE2). 5−7 S protein consists of multiple functional domains, including the receptor binding domain (RBD) that is responsible for interacting with ACE2. The RBDs on the top of S trimer are conformationally variable. In the so-called "closed state", the RBDs lay flat with their receptor binding motifs (RBMs) occluded by the RBDs of neighboring protomers. The "open states" are characterized by one or more uplifted RBDs, resulting in exposure of their RBMs ( Figure 1A ). Many NAbs have been isolated from the sera of recovered COVID-19 patients, and most of them target the RBD. The structures of S trimers in complex with various antibodies have been solved by cryogenic electron microscopy (cryo-EM). 8−12 Barnes et al. classified the structures of human NAbs solved by themselves and those described in other published studies into four categories. Class I directly occupies the ACE2 binding site and binds only to open-state RBDs (blue in Figure 1B ). Class II has a different epitope, but it also blocks binding of ACE2. It binds to both open-and closedstate RBDs and has contacts with adjacent RBDs (green). Class III targets the outside of RBDs and binds to both openand closed-state RBDs (orange). Class IV targets the inside of RBDs and binds only to open-state RBDs (magenta).

Although these static structures provide valuable information of S−antibody interactions at near-atomic resolution, they may not include all the information that is necessary for understanding the mechanisms underlying antibody binding. In addition, many cryo-EM structures have missing residues in S RBDs and/or antibodies, and the glycans that can have significant influence on antibody binding are mostly not resolved. Therefore, molecular dynamics (MD) simulations based on well-refined initial structures with all missing portions properly modeled can provide a more comprehensive understanding of S−antibody interactions. As an RNA virus, SARS-CoV-2 has a relatively high mutation rate, which is a big challenge to the efficacy of antibodies and vaccines. In the past one and a half years, multiple variants of SARS-CoV-2 have appeared and are now circulating globally. In our recent study, steered MD simulations combined with microscale thermophoresis experiments have been used to show that different combinations of mutations on RBDs lead to various binding affinities between ACE2 and SARS-CoV-2 variants. 13 The US Centers for Disease Control and Prevention classify the circulating variants into variants of concern including Alpha (first identified in United Kingdom, B. In this work, we present all-atom MD simulations of a fully glycosylated S protein trimer in complex with various antibodies. We have selected antibodies that target different epitopes on the RBD ( Figure 1B ) and modeled the structures of S trimers bound with these antibodies ( Figure 1C ). Thirteen systems consisting of different antibodies and S trimers with different RBD open/closed states have been built and simulated (Table 1) . For convenience, a one-letter symbol "O" or "C" is used to represent an open-or closed-state RBD, and thus antibody C002 bound to the S trimer with open-closed-closed RBDs is denoted as "C002_OCC." To sample as many S−antibody binding interfaces as possible in each simulation, we modeled all three RBDs bound with antibodies if the specific epitopes are not occluded and the modeled antibodies do not result in significant steric hinderance. The only exception is antibody C105 whose epitope is not exposed when the RBD is closed (Figures S1−S4) . Note that the reference PDB structures contain some free RBDs without bound antibodies. Our results provide insight into the dynamics of antibodies, contributions of residues for binding interactions, and roles of glycans in antibody binding, which provides a better understanding of S−antibody interactions.

Modeling of Fully Glycosylated SARS-CoV-2 S Protein−Antibody Complex Structures. As reported in our previous work, 16, 17 we have modeled a fully glycosylated full-length SARS-CoV-2 S protein structure using a combination of computational modeling tools including the GALAXY protein modeling suite 18−20 for building missing residues and In this study, we truncated the heptad repeat linker, transmembrane domain, and cytoplasmic domain built by ab initio structure prediction, and only kept the S1 subunit and part of the S2 subunit. Before generating the S− antibody complex structures, we first removed all glycans to avoid bad contacts with antibodies. We extracted the portion of the RBD−antibody from each cryo-EM structure of the S trimer−antibody complex (Table 1 ) and superimposed it onto our S trimer structure by maximizing the overlap in RBD. For example, when building the C121_OCC system, we used the cryo-EM structure (PDB id: 7K8X) as a reference. PDB 7K8X has two closed-state RBDs both bound with antibody C121 and one open-state unbound RBD. We simply extracted the portions of RBD−antibody from two closed chains in 7K8X and aligned them onto the closed chains in our S trimer structure. Because the open chain in 7K8X is unbound, we used the portion of RBD−antibody extracted from the closed chain to model the binding interface in the open chain with minor clashes. The GALAXY protein modeling suite was then used to relax the structure to relieve the clashes. After the S trimer−antibody complex structure was generated, CHARMM-GUI Glycan Reader & Modeler was used to build 19 N-linked and 1 O-linked glycans onto each protomer using the same glycoforms as those in our previous model (Table  S1 ). Simulation Details. In this study, the CHARMM36(m) force field was used for proteins 29 and carbohydrates. 30−32 The TIP3P water model 33 was utilized along with a 0.15 M KCl solution. The total number of atoms is approximately 1,250,000 (∼400,000 water molecules, ∼1100 K + , and ∼1100 Cl − ), but the exact numbers differ between various systems. The van der Waals interactions were smoothly switched off over 10−12 Å by a force-based switching function, 34 and the long-range electrostatic interactions were calculated using the particle-mesh Ewald method 35 with a mesh size of ∼1 Å. To avoid the overestimation of protein− protein interactions, a force field adjustment was made to enhance protein−water interactions. 36, 37 All simulations were performed using the input files generated by CHARMM-GUI, 38−40 and we used GROMACS 2018.6 41 for both equilibration and production with the LINCS algorithm. 42 The temperature was maintained using a Nose−Hoover temperature coupling method 43,44 with a τ t of 1 ps. For pressure coupling (1 bar), a semi-isotropic Parrinello− Rahman method 45, 46 with a τ p of 5 ps and a compressibility of 4.5 × 10 −5 bar −1 was used. During the equilibration run, NVT (constant particle number, volume, and temperature) dynamics was first applied with a 1 fs time step for 250 ps. Subsequently, the NPT (constant particle number, pressure, and temperature) ensemble was applied with a 1 fs time step (for 125 ps) and with a 2 fs time step (for 1.5 ns). During the equilibration, positional and dihedral restraint potentials were applied, and their force constants were gradually reduced. The production run was performed with a 4 fs time step using the hydrogen mass repartitioning technique 47 without any restraint potential. Each system shown in Table 1 ran for 500 ns production time with about 20 ns/day using 768 CPU cores on NURION in the Korea Institute of Science and Technology Information. For comparison, we also ran 500 ns production of 12 S-only systems (3 for each of CCC, OCC, OOC, and OOO) with the same simulation protocols. All 4 S-only and 13 S−antibody simulation systems and trajectories are available in the CHARMM-GUI COVID-19 archive (https://www. charmm-gui.org/docs/archive/covid19).

Simulation Analysis. Interacting Residue Pairs. We first used the MDTraj python library 48 to identify all heavy atom pairs consisting of one atom from S and the other from the antibody within 5 Å, and then we checked whether they formed favorable interactions such as hydrophobic interactions, π−π stacking interactions, hydrogen bonds, and salt bridges. If two residues contain at least one pair of interacting atoms, we consider them an interacting residue pair.

RBD Motion and Root-Mean-Square Fluctuation (RMSF). We measured two structural features to describe the RBD motion: the RBD-NTD (N-terminal domain) distance (d) defined by the minimum distance between the heavy atoms of RBD (N334 to P527) and those of NTD (C15 to S305), and the RBD orientation angle (θ) defined by three points corresponding to the (i) center of mass (COM) of L452 and L492, (ii) COM of N334, and (iii) COM of S1030. To calculate the RMSF, we first aligned the simulating trajectory onto RBD (N334 to P527) in each protomer and calculated the RMSF of that specific RBD using the MDTraj python library.

Critical Residues for Antibody Binding and Potential Influence of Mutant Variants. To identify critical residues for S−antibody interactions, we searched for all residue pairs consisting of one residue from S and the other from antibody, which make favorable interactions including hydrophobic interaction, π−π stacking interaction, hydrogen bonds, and salt bridges. We processed 500 snapshots (every 1 ns) from the 500 ns simulation trajectory of each system and calculated the frequency of interacting residue pairs. Figure 2 shows two representative cases of the antibodies bound to the RBDs. The first one is the interface of C105 bound to chain A in C105_OOC. There are 30 residues of RBD_A (RBD of chain A) that interact with the antibody in at least 10% of snapshots, and more than 10 of them keep their interactions with the antibody in at least 75% of snapshots ( Figure 2A ). Two of these critical residues are mutated in COVID-19 variants, which are K417N in the Beta variant, K417T in the Gamma variant, and N501Y in Alpha, Beta, and Gamma variants. K417 has high frequency of interactions with Y33, Y52, and E96 from C105 V H (variable domain of antibody heavy chains). As shown in Figure 2C , the positively charged side chain of K417 can make a salt bridge with E96 or cation−π interaction with Y33 or Y52 depending on the conformation of these four residues, and such interactions are lost with K417N or K417T mutation. This implies that K417N/K417T will have significant effects on the interaction with C105, which is likely to reduce the efficacy of C105. N501 has high frequency of interactions with several residues from C105 V L (variable domain of the antibody light chain). Different from K417, N501 mostly interacts with C105 by the backbone atoms, and the side chain points toward the RBD itself. If a mutation does not significantly change the structure and dynamics of this local backbone, it should not be harmful to the binding of C105. However, N501Y introduces a new side chain with a larger size, and the local backbone is a flexible loop. The backbone conformation may change in order to Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article accommodate the side chain of Tyr. Therefore, it is possible that N501Y can affect the antibody binding. The second representative case is the interface of C119 bound to chain C in C119_OCC. There are more than 10 residues of RBD_C continuously interacting with C119 V H and/or V L in the simulation trajectory ( Figure 2B ). Two critical residues involved in COVID-19 variants, L452R and E484Q, occur simultaneously as a double mutation in the Kappa variant. E484K also appears in Beta and Gamma variants as well as some sequences but not all in the Alpha variant, and L452R also appears in Delta and Epsilon variants. As shown in Figure 2D , the hydrophobic side chain of L452 continuously interacts with the aromatic side chain of Y29 from C119 V L . With the L452R mutation, the positively charged side chain can still interact with the aromatic ring (i.e., cation−π interaction). However, because the side chain of L452 points toward C119 V L and the space between RBD and C119 V L is limited, residue 452 is very likely to clash with the antibody after mutating to Arg that has a longer side chain. Consequently, L452R may reduce the efficacy of all class II antibodies that interact with L452. This is also mentioned in the study by McCallum et al., where the RBD−antibody complex structures were superimposed onto B.1.427/B.1.429 S structures. 49 Another critical residue is E484 which has a salt bridge with K30 from C119 V L . The E484K mutation not only breaks the salt bridge but also forms an unfavorable contact of two positively charged groups. Therefore, E484K may also seriously impair the efficacy of antibody C119.

The interacting residue pairs of all the systems are shown in Figures S5−S10. In summary, antibody C105 is likely to be sensitive to mutation K417N/K417T as the salt bridge and cation−π interaction involving K417 is frequently observed in all cases except chain A of C105_OOO. However, the influence of N501Y is uncertain, and the interaction between N501 and antibody is only observed in chain A of C105_OOC. Antibodies C002, C119, and C121 have frequent interactions with L452 and/or E484, which is observed in all cases except chain A of C121_OCC. Because all variants of current concern contain at least one of two mutations, it remains a question whether these three antibodies are effective enough to the variants. Antibodies S309 and EY6A target different epitopes that do not contain the mutant residues, and hopefully, they would not be affected by the mutations in COVID-19 variants. The maps of interacting pairs derived from our simulation results highlight the critical residues from RBD and their interacting residues from antibodies. By considering the physicochemical properties of wild-type and mutant amino acids as well as their possible interactions with the antibody residues, we are able to predict the potential impacts of S mutations on antibody binding.

Antibodies Can Have Two Binding Interfaces with RBDs from Different Protomers. The class II antibodies (C002, C119, and C121) bind to similar but not identical epitopes, located around the green region shown in Figure 1 , with different interaction patterns. When these antibodies bind to a closed-state RBD (e.g., RBD_C), they can also have a secondary binding interface with RBD (+1), the next RBD in the anticlockwise direction from top view (e.g., RBD_A). To examine the existence of such a secondary binding interface for other systems, we performed the analysis of interacting residue pairs with consideration of both the RBD containing the primary interface and the neighboring RBD that may contain the secondary interface ( Figures S5−S10) . Except for C002, C119, and C121, all other antibodies have only one primary binding interface. C002, C119, and C121 have a secondary binding interface regardless of whether the neighboring RBD is open or closed, but the secondary interface has more extensive interacting residue pairs when the neighboring RBD is open ( Figures S6−S8) . S309 occasionally has minor contacts with RBD (−1), the preceding RBD in the anticlockwise direction, when both RBDs are closed ( Figure S9) . Figure 3 shows two representative cases of the antibodies interacting with two RBDs. The first is C121 bound to closedstate RBD_C in C121_OCC. C121 has a secondary binding interface with the neighboring open-state RBD_A. Figure 3A shows that V H makes the major contribution to RBD_C binding. In addition, there are multiple residues from both C121 V H and V L that continuously interact with the residues from RBD_A. This suggests that C121 binds to S protein through multivalent interactions which can lead to enhancement of binding affinity. As shown in Figure 3C , both complementarity-determining region 3 (CDR3) of V H (red circle) and CDR1/CRD2 of V L (purple circle) are involved in the secondary binding interface. When C121 binds to RBD_B (primary) and RBD_C (secondary) that are both closed, the secondary interface consists of only a couple of residues from CDR3 of V H , but two of them (V106 and L107) continuously interact with RBD_C ( Figure S8A middle) .

The second representative case is C002 bound to closedstate RBD_C (primary) and open-state RBD_A (secondary) in C002_OCC. Although multiple residues from V L are involved in the secondary interface, none of the interacting residue pairs has a frequency over 60%, indicating that the residue contacts in the secondary interface keep varying with the movements of V L and RBD_A ( Figure 3B ). As shown in Figure 3D , the antibody residues involved in the secondary interface (orange circle) are not within any CDR of V L . Therefore, the binding interaction between C002 and RBD_C/RBD_A is not multivalent, and the secondary interface appears to exist simply because the binding pose of C002 in the primary interface and the relative position of two RBDs leads to such fortuitous contacts.

RBD−Antibody Interfaces Are Stable in All Systems. Although there exist cryo-EM structures for each antibody, many of them have only one or two antibodies bound to RBD(s). In this study, we aimed to build the S−antibody complex systems with all three RBDs bound with antibodies, except for C105 whose epitope is not exposed when RBD is closed. Therefore, we needed to check whether binding interfaces are stable considering the mutual influence between multiple antibodies. For each system, we measured the minimum distance between the heavy atoms of V H /V L and those of RBD and calculated the number of interacting residue pairs between V H /V L and RBD (Figures S11−S16). For the antibodies with primary and secondary interfaces, we considered only the primary interfaces. All antibodies stably stay in the binding sites and continuously interact with the Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article epitopes on RBDs (at least during the current simulation time). We observed three different binding modes: (i) V H and V L make similar contributions to RBD binding ( Figure 4A ), (ii) V H makes the major contribution to RBD binding, but V L continuously interacts with RBD ( Figure 4B ), and (iii) V H makes the major contribution to RBD binding and V L does not always make contacts with RBD ( Figure 4C) . Figure 4D −F shows the snapshots extracted from the trajectories of C119, C002, and C121. Both V H and V L of C119 extensively interact with RBD. In contrast, for C002 and C121, V H makes the major contribution for binding while only CDR3 of V L has contacts with one loop in RBD. C121 V L is distant from RBD; hence, CDR3 of V L can interact with RBD only if the RBD loop is in certain conformations. A summary of binding modes in each system is shown in Table 2 . Clearly, S−antibody binding modes are versatile and specific to individual antibodies and RBD's open/closed states. Glycan Flexibility and Its Roles in Antibody Binding. S protein is highly glycosylated, and it has been shown that the glycans can play important roles in protein stability, RBD open-closed state transition, and interactions with ACE2 and antibodies. 50−52 In our previous work, 17 we aligned the PDB structures of RBD−antibody complexes onto the simulation trajectories of fully glycosylated S trimers and investigated the influences of glycans on antibody binding by measuring the extent of clashes between glycans and superimposed antibodies. In this work, we explored the same question by comparing the glycan flexibility in S-only and S−antibody complex systems. Three glycans attached to N165, N331, and N343 are close to the epitopes marked with green and orange in Figure  1 . All of them could possibly influence antibody binding depending on the glycan conformation and open/closed state of RBD. Among these glycans, N343 glycan is attached to the residue within the orange epitope, and N165 glycan continuously interacts with the green epitope; hence, they are more likely to be involved in antibody binding activity. We measured the φ/ψ angles of the first three glycosidic linkages from the reducing terminal (i.e., the linkages in Manβ1− 4GlcNAcβ1−4 GlcNAcβ1-Asn) for N343 in S309_CCC and N165 glycan in C119_CCC, as well as these two glycans in the S-only system with three closed-state RBDs. The dihedral angles of GlcNAcβ1-Asn linkage are defined by φ: (2) is Asn. The dihedral angles of Manβ1− A  C105_OOO  3  2  2  C002  C002_OCC  2  2  2  C002_OOC  2  2  3  C119  C119_CCC  1  1  1  C119_OCC  1  1  1  C119_OOC  1  1  1  C121  C121_OCC  1  2  3  C121_OOC  3  3  3  S309  S309_CCC  2  2  2  S309_OCC  3  2  2  EY6A  EY6A_OOO  1  2  2 a Three different binding modes are observed; (1) both V H and V L make similar contributions to RBD binding, (2) V H makes the major contribution to RBD binding, but V L continuously has contacts with RBD, and (3) V H makes the major contribution to RBD binding and V L can move away from the RBD. (2) is the reducing terminal unit. The illustration is shown in Figure  S17 .

We superimposed the N343 glycan in each snapshot onto the initial structures by aligning the local backbone consisting of the glycosylation site and two neighboring residues, and the resulting structures of S-only and S−antibody systems are shown in Figure 5A . Figure 5B shows the distributions of dihedral angles in N343 glycosidic linkages, and the glycans attached to RBD_A, RBD_B, and RBD_C are shown in different colors. By comparing two rows, it becomes evident that the conformational space explored by N343 glycan in the S-only system is similar to that in the S−antibody system. Only two minor regions that are explored in the S-only systems (labeled by red arrows) are not explored in the S−antibody systems. In the S-only system, N343 glycan interacts with the closed-state RBD on the neighboring chain in most of the simulation time. This interaction keeps N343 glycan pointing toward the left direction ( Figure 5A ), which does not disturb the binding of S309 antibody. On the contrary, the core portion of the glycan can interact with the antibody, serving as part of the epitope. Only very occasionally, it switches to the other direction and shields the binding interface for the antibody. Therefore, N343 glycan facilitates the binding of antibody S309 rather than blocking such interaction, which is consistent with the conclusion from our previous work. 17 The frequently interacting residue pairs between S309 antibody and N343 glycan are shown in Figure S18 . In four out of six cases (considering three chains in S309_CCC and three chains in S309_OCC), we observe that multiple antibody residues frequently interact with N343 glycan, particularly with fucose. This is consistent with the results reported in the work of Pinto et al. showing that V H CDR3 and V L CDR2 sandwich the core fucose of N343 glycan. 15 In the remaining two cases, the fucose points toward the opposite direction of the antibody in the initial structure, and it did not rotate to make contacts with the antibody during the current limited simulation time. Instead, the antibody has extensive interactions with the α1−3 branch of N343 glycan. Figure 5D shows the comparison of dihedral distributions of N165 glycan in S-only and S−antibody systems. There are two major regions (labeled by purple arrows) that are explored in the S-only system but not in the S−antibody system. As shown in the superimposed glycan conformations ( Figure 5C ), when the antibody exists, the glycan is forced to point toward the left direction. Note that we built the S−antibody system by first generating the S trimer−antibody complex and then modeling the glycans, and consequently we observed that the antibody restricted the glycan motion. However, if we compare this with the simulation of the unbound S glycoprotein, it becomes clear that N165 glycan occupied the binding interface for antibody C119 in about half of the simulation time and moved out of Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article that region in the remaining time. These indicate that N165 glycan does not always block the antibody binding, but the existence of this glycan makes it more difficult for the antibody to get access to the epitope. We also performed the same analysis for N343 glycan in C119_CCC, N331 glycan in S309_CCC, N165 glycan in S309_CCC, and N165 glycan in EY6A_OOO. When N343 glycan in C119_CCC was modeled onto the S trimer in the S−antibody complex system, the glycan had different conformations from that in the S-only system, because of the limited space between antibody and RBD ( Figure S19 ). The ψ angle of linkage #1 is sampled around 0°in S−antibody systems, while it is sampled around 180°in S-only systems, showing a huge conformational difference in the entire glycan. Although the glycan can occlude the epitope and block antibody binding when it is in the conformations with the ψ angle of linkage #1 ≈ 180°, it only has slight clashes with the antibody in a very small number of snapshots when the ψ angle of linkage #1 is around 0°. This again suggests that the glycan can disturb the antibody binding, but the epitope becomes accessible to the antibody when the glycan is in certain conformations. The result of N331 glycan in S309_CCC shows that it does not affect the antibody binding ( Figure  S20 ). For N165 glycan in S309_CCC, the glycan in the initial structure had two different conformations reflected by the ψ angle of linkage #1 (one observed in chains A and C, and the other in chain B) in S-only systems, and conformational transition did not occur during the simulation. The glycan in the S−antibody complex system was built with only one conformation, and the dihedral space sampled is similar to that of chains A and C in S-only systems ( Figure S21 ). Therefore, N165 glycan does not affect the binding of S309 antibody either. Finally, N165 glycan is very likely to disturb the binding of EY6A antibody. The glycan occupies the space for EY6A antibody in the S-only systems, and it is trapped in two different regions (one observed in chain A and the other in chains B and C) in S−antibody complex systems ( Figure S22) .

Effects of Antibody Binding on RBD Motion and Flexibility. To investigate whether antibody binding has significant influences on RBD motion, we measured two structural features: the RBD-NTD distance (d) defined by the minimum distance between the RBD and NTD and the RBD orientation angle (θ) defined by three points on the RBD and S-trimer central axis (see Methods). θ reflects the open/closed state of RBD, and d reflects the size of the U-shape pocket that accommodates the neighboring closed-state RBD. It is expected that antibody binding could stabilize the S trimer and make the RBD open/closed state transition more difficult, particularly for the antibodies having two interfaces with different protomers. For all-closed and one-open S trimers, the RBDs are generally stable and show very limited motions even in the S-only systems; hence antibody binding has generally no influence on RBD orientation and motion during the current simulation time (Figure S23A In addition to RBD motion, we also measured the RMSF of RBD in S-only and S−antibody systems ( Figures S25−S28) . Though there exist differences between different protomers, the open-state RBDs generally have a higher RMSF than the closed-state RBDs. The antibodies target both highly flexible regions and stable regions. There are two large peaks around residues 445 and 480, which are two loops on the top of RBD. They become much more stable upon antibody binding, particularly when these regions are within the epitopes of antibodies (i.e., C002, C119, and C121).

In this work, we have presented a modeling and simulation study of fully glycosylated S protein in complex with various antibodies. The frequency of interacting residue pairs between S and antibody reveals the interaction patterns and the critical S protein residues contributing to antibody binding, which enables us to predict possible effects of mutations in SARS-CoV-2 variants. Such analyses also make it possible to identify the antibodies that bind the S trimer through multivalent interactions. Comparison of glycan conformational flexibility in S-only and S−antibody complex systems reveals whether a specific glycan disturbs antibody binding. During the current limited simulation time, we could not observe that antibodies had any significant influence on RBD open-closed state transition, but the RMSF of the RBD was reduced because of antibody binding. Together, this study provides richer insight into the S−antibody interactions that is not available from the static cryo-EM structures, and we hope that our findings are useful for the development of vaccines and antibody-based therapeutic agents.

The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.jctc.1c00552.

Table S1 Glycosylation sites, selected glycan compositions, and glycan sequences in this study, Figure S1 . Illustration of C105_OCC, C105_OOC, and C105_OOO, Figure S2 Illustration of C119_CCC, C119_OCC, and C119_OOC, Figure S3 Illustration of S309_CCC, and S309_OCC, Figure S4 Illustration of EY6A_OOO, Figure S5 The frequency of interacting residue pairs in C105 systems, Figure S6 The frequency of interacting residue pairs in C002 systems, Figure S7 The frequency of interacting residue pairs in C119 systems, Figure S8 The frequency of interacting residue pairs in C121 systems, Figure S9 The frequency of interacting residue pairs in S309 systems, Figure S10 The frequency of interacting residue pairs in EY6A_OOO, Figure S11 RBD−antibody distance and contact residue number in C105 systems, Figure S12 RBD−antibody distance and contact residue number in C002 systems, Figure S13 RBD−antibody distance and contact residue number in C119 systems, Figure S14 RBD−antibody distance and contact residue number in C121 systems, Figure S15 . RBD−antibody distance and contact residue number in S309 systems, Figure S16 RBD−antibody distance and contact residue number in EY6A_OOO, Figure S17 Illustration of φ/ψ angles in the first three glycosidic linkages from the reducing terminus, Figure S18 Frequency of interacting residue Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article pairs between S309 antibody and N343 glycan, Figure  S19 Glycan flexibility in S-only and S−antibody complex systems: N343 glycan in C119_CCC, Figure S20 Glycan flexibility in S-only and S−antibody complex systems: N331 glycan in S309_CCC, Figure S21 Glycan flexibility in S-only and S−antibody complex systems: N165 glycan in S309_CCC, Figure S22 Glycan flexibility in S-only and S−antibody complex systems: N165 glycan in EY6A_OOO, Figure S23 NTD-RBD distance (d) and RBD orientation angle (θ) in S-only and S−antibody complex systems, Figure S24 Representative structures with various NTD-RBD distances (d) and RBD orientation angles (θ), Figure S25 RMSF of RBD in S-only and S−antibody complex systems with all three RBDs closed, Figure 

Safety and Efficacy of the BNT162b2 mRNA Covid-19 Vaccine

Efficacy and Safety of the mRNA-1273 SARS-CoV-2 Vaccine

BNT162b2 mRNA Covid-19 Vaccine in a Nationwide Mass Vaccination Setting

The Molecular Biology of Coronaviruses

Functional assessment of cell entry and receptor usage for SARS-CoV-2 and other lineage B betacoronaviruses

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

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

Structures of Human Antibodies Bound to SARS-CoV-2 Spike Reveal Common Epitopes and Recurrent Features of Antibodies

SARS-CoV-2 neutralizing antibody structures inform therapeutic strategies

Structural basis for the neutralization of SARS-CoV-2 by an antibody from a convalescent patient

Potent neutralizing antibodies against multiple epitopes on SARS-CoV-2 spike

Differential Interactions Between Human ACE2 and Spike RBD of SARS-CoV-2 Variants of Concern

VMD: visual molecular dynamics

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

Developing a Fully Glycosylated Full-Length SARS-CoV-2 Spike Protein Model in a Viral Membrane

Receptor Binding, and Antibody Binding of the Fully Glycosylated Full-Length SARS-CoV-2 Spike Protein in a Viral Membrane

GalaxyTBM: template-based modeling by building a reliable core and refining unreliable local regions

The FALC-Loop web server for protein loop modeling

Symmetric and asymmetric ab initio protein−protein docking web server with improved energy parameters

ISOLDE: a physically realistic environment for model building into low-resolution electron-density maps

Glycan Reader: automated sugar identification and simulation preparation for carbohydrates and glycoproteins

Glycan Reader is improved to recognize most sugar types and chemical modifications in the Protein Data Bank

CHARMM-GUI Glycan Modeler for modeling and simulation of carbohydrates and glycoconjugates

Automated Builder and Database of Protein/Membrane Complexes for Molecular Dynamics Simulations

CHARMM-GUI Membrane Builder toward realistic biological membrane simulations

CHARMM-GUI Membrane Builder for Complex Biological Membrane Simulations with Glycolipids and Lipoglycans

CHARMM36m: an improved force field for folded and intrinsically disordered proteins

Additive empirical force field for hexopyranose monosaccharides

CHARMM Additive All-Atom Force Field for Glycosidic Linkages between Hexopyranoses

CHARMM Additive All-Atom Force Field for Aldopentofuranoses, Methylaldopentofuranosides, and Fructofuranose

Comparison of simple potential functions for simulating liquid water

New spherical-cutoff methods for long-range forces in macromolecular simulation

A smooth particle mesh Ewald method

Clustering and dynamics of crowded proteins near membranes and their influence on membrane bending

Slow-Down in Diffusion in Crowded Protein Solutions Correlates with Transient Cluster Formation

CHARMM-GUI: a webbased graphical user interface for CHARMM

Simulations Using the CHARMM36 Additive Force Field

CHARMM-GUI supports the Amber force fields

GROMACS: Fast, flexible, and free

LINCS: A linear constraint solver for molecular simulations

A molecular dynamics method for simulations in the canonical ensemble

Canonical dynamics: Equilibrium phase-space distributions

Polymorphic transitions in single crystals: A new molecular dynamics method

Constant pressure molecular dynamics for molecular systems

Long-Time-Step Molecular Dynamics through Hydrogen Mass Repartitioning

A Modern Open Library for the Analysis of Molecular Dynamics Trajectories

The Impact of Mutations in SARS-CoV-2 Spike on Viral Infectivity and Antigenicity

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

Beyond Shielding: The Roles of Glycans in the SARS-CoV-2 Spike Protein

The authors declare no competing financial interest.