Impact of new variants on SARS-CoV-2 infectivity and neutralization: A molecular assessment of the alterations in the spike-host protein interactions.

The emergence of SARS-CoV-2 variants necessitates rational assessment of their impact on the recognition and neutralization of the virus by the host cell. We present a comparative analysis of the interactions of Alpha, Beta, Gamma, and Delta variants with cognate molecules (ACE2 and/or furin), neutralizing nanobodies (Nbs), and monoclonal antibodies (mAbs) using in silico methods, in addition to Nb-binding assays. Our study elucidates the molecular origin of the ability of Beta and Delta variants to evade selected antibodies, such as REGN10933, LY-CoV555, B38, C105, or H11-H4, while being insensitive to others including REGN10987. Experiments confirm that nanobody Nb20 retains neutralizing activity against the Delta variant. The substitutions T478K and L452R in the Delta variant enhance associations with ACE2, whereas P681R promotes recognition by proteases, thus facilitating viral entry. The Ab-specific responses of variants highlight how full-atomic structure and dynamics analyses are required for assessing the response to newly emerging variants.

Introduction

The high transmissibility and mutation rates of SARS-CoV-2, together with the lack of a robust preexisting immunity by hosts (Petersen et al., 2020; Tortorici et al., 2020) and slow rate of immunization through vaccinations, have caused COVID-19 cases to surge to more than 240 million worldwide by the end of October 2021. Although several antibodies (Abs) proved to neutralize the wild-type (WT) virus, their efficacies against emerging variants ought to be carefully monitored, especially against those noted to possibly evade the host immune response (Andreano et al., 2021; Kemp et al., 2021; McCarthy et al., 2021; Meng et al., 2021; Planas et al., 2021).

The receptor binding domain (RBD; residues 331–524) of SARS-CoV-2 spike (S) glycoprotein has been a focal point in many studies due to its role in binding to the human receptor ACE2 as well as recognition by the immune system (Belizario, 2021). Many mutations in the emerging variants are located in the RBD, and not surprisingly, the RBD is targeted by the majority of mAbs from COVID-19 patients and mAb therapies under investigation. Several RBD residues underwent mutations arisen independently in multiple lineages (Plante et al., 2021). These include N501Y in the Alpha, Beta, Gamma, and Omicron variants; K417 mutated to N in the Beta, Delta, and Omicron variants and to T in the Gamma variant; and E484 mutated to K in the Beta and Gamma variants and to A in the Omicron variant.

SARS-CoV-2 variants have been studied in deep scanning mutagenesis (Dong et al., 2021; Greaney et al., 2021a, 2021c; Starr et al., 2020, 2021a, 2021b), in vitro evolution (Liu et al., 2021d; Muecksch et al., 2021; Wang et al., 2021c; Weisblum et al., 2020; Zahradnik et al., 2021), cryo-EM (Cai et al., 2021; Dejnirattisai et al., 2021b; Du et al., 2021; Gobeil et al., 2021; Han et al., 2021; Liu et al., 2021a; McCallum et al., 2021a, 2021b; Pymm et al., 2021; Supasa et al., 2021; Wang et al., 2021a; Zhang et al., 2021a; Zhu et al., 2021), computational studies (Brown et al., 2021; Cheng et al., 2021a; Fratev, 2021; Khan et al., 2021; Villoutreix et al., 2021), and Ab-binding experiments (see below). There have been comprehensive attempts to map out different classes of Abs and relate them to mutations (Barnes et al., 2020a, 2020b; Cerutti et al., 2021; Dejnirattisai et al., 2021a; Greaney et al., 2021b; Liu et al., 2020; Rogers et al., 2020; Wang et al., 2021b; Wheatley et al., 2021; Xiang et al., 2020; Yuan et al., 2021), but a systematic comparative analysis of the behavior of selected variants using structure-based simulations and analytical methods is lacking.

In addition to the RBD, another site of critical importance is the S1/S2 cleavage site P681RRAR↑S of the spike (Hoffmann et al., 2020), which is implicated in multiple roles. The cleavage of the spike monomers at this site is an essential step priming the subunit S2 to its fusogenic activity for viral entry. Of interest is to understand how mutations at this site affect the interaction with host cell proteins such as furins involved in the proteolytic cleavage. Furthermore, this region belongs to a superantigenic site shown in earlier work to trigger an adaptive immune response that could escalate into hyperinflammation and cytokine storm typical of the multisystem inflammatory syndrome (MIS) observed in children and severely infected adults (Cheng et al., 2020, 2021b). Of interest is mutations at P681, to R in the Delta variant (Sheikh et al., 2021) and to H in Alpha and Omicron, which lies within this superantigen-like motif and has been proposed to be a key site in SARS-CoV-2 infection and pathogenesis (Johnson et al., 2021).

In the present study, we present results from structure-based in silico modeling and full-atomic simulations to assess the interactions of the variants listed in Table 1 with cognate proteins essential to viral entry and their ability to escape neutralizing monoclonal antibodies (mAbs) and nanobodies (Nbs) listed in Table 2. We focus on the Alpha (Kent, UK), UK2 (Bristol, UK), Beta (South Africa), and Delta (India) variants. Our study indicates higher affinity of all these variants to bind to ACE2, compared with WT spike. The study further provides a molecular rationale for experimentally observed interactions in the neutralizing potential of specific mAbs and Nbs. Examination of the effect of L452R shared by the Delta and Epsilon (Californian) variants shows that this mutation weakens interaction with, and thereby reduces the efficacy of, selected mAbs (e.g. Eli Lilly LY-CoV555) (Starr et al., 2021b), whereas an ultrapotent inhalable nanobody Nb20 (Xiang et al., 2020) shows a comparable, if not slightly increased, affinity to the Delta variant. Finally, the Delta variant exhibits high affinity to ACE2 and furin, thus facilitating viral entry.

Table 1

SARS-CoV-2 variants and corresponding mutations at the spike proteina

Variant	Lineage	First observed in	Notable mutations in the spike
Alpha	B.1.1.7	Kent, UK	N501Y, P681H, 69-70del, and D614G
UK2	B1.1.7 with E484K	Bristol, UK	N501Y, P681H, 69-70del, E484K and D614G
Beta	B.1.351	South Africa	N501Y, E484K, K417N, and D614G
Gamma	B.1.1.28, P.1	Brazil	N501Y, E484K, K417T, and D614G
Delta	B.1.617.2	India	L452R, T478K, P681R, and D614G

Countries of origin are listed, but these variants have now been observed in multiple countries (Kemp et al., 2021; Plante et al., 2021; Tegally et al., 2020; Zhang et al., 2021c).

Table 2

Complexes of SARS-CoV-2 spike (or RBD) with Abs and host proteins analyzed in the present study

mAb/Nb (PDB id)a	Residues at the mAb/N-binding epitope on SARS-CoV-2 spikeb	Refc	Change in neutralization effect				Epitope classf
			Current ΔΔGbinding prediction (in kcal/mol)d		Literaturee
			Beta	Delta	Beta	Delta
REGN10933 casirivimab (6XDG)	K417, Y453, L455,F456, E484, G485, F486,N487, C488, Y489, Q493	(Hansen et al., 2020)	1.93 (↓)	1.06 (↓)	(Wang et al., 2021b) ↓	(Tada et al., 2021) ↓(Mlcochova et al., 2021) (−)	I1,E
C105 (6XCN, 6XCM)	D405,T415, G416, K417, Y421, Y453, F456, R457, K458, N460, Y473, A475, G476, F486, N487, N501, G502, Y505	(Barnes et al., 2020b)	0.42 (↓)	0.10 (↓)	(Greaney et al., 2021b) ↓
LY-CoV016 etesevimab (7C01)	R403, D405, E406, R408, Q409, T415, G416, K417, D420, Y421, L455, F456, R457, K458, N460, Y473, Q474, A475, G476, S477, F486, N487, Y489, Q493, Y495, G502, Y505	(Shi et al., 2020)	−1.86 (↑)	0.06 (−)	(Arora et al., 2021; Wang et al., 2021b)↓	(Connor et al., 2021; Mlcochova et al., 2021) (−)	II1,E/F
					(Chen et al., 2021a) (−)
B38 (7BZ5)	R403, D405, Q409, T415, G416, K417, D420, Y421, Y453, L455, F456, R457, K458, N460, Y473, A475, G476, S477, F486, N487, Y489, F490, Q493, Y495, G496, Q498, T500, N501, G502 Y505	(Wu et al., 2020)	0.61 (↓)	1.13 (↓)	(Chen et al., 2021b) ↓
BD23 (7BYR)	G446, Y449, E484, G485, F486, Y489, F490, L492, Q493, S494, G496, Q498 N501, Y505, N165 glycan	(Cao et al., 2020)	−0.94 (↑)	0.33 (↓)	(Ray et al., 2021)∗ ↑	(Ray et al., 2021)∗ ↑	III2,F
H11-H4 (6ZBP, 6ZHD, 6ZH9)	Y449, N450, L452, L455, F456, T470, G482, V483, E484, Y489, F490, L492, Q493, S494	(Huo et al., 2020)	0.32 (↓)	0.61 (↓)	(Huo et al., 2021) ↓	(Golcuk et al., 2021)∗ ↓
LY-CoV555 bamlanivimab (7KMG)	Y351, Y449, N450, L452, L455, T470, N481, G482, V483, E484, G485, F486, C488, Y489, F490, L492, Q493, S494	(Jones et al., 2021)	0.29 (↓)	0.31 (↓)	(Hoffmann et al., 2021; Wang et al., 2021b)↓	(Mlcochova et al., 2021; Planas et al., 2021)↓
Nb20 (7JVB)	Y449, N450, L452, Y453, L455, V483, E484, Y489, F490, L492, Q493, S494	(Xiang et al., 2020)	0.28 (↓)	−0.15 (↑)	(Sun et al., 2021)↓
REGN10987 imdevimab (6XDG)	R346, N440, L441, K444, V445, G446, N448, Y449, Q498	(Hansen et al., 2020)	_ _				IV3, H
EY6A (6ZDH)	Y369, F374, S375, T376, F377, K378, C379, Y380, G381, V382, S383, P384, T385, K386, D389, L390, F392, P412, G413, D427, D428, F429, T430	(Zhou et al., 2020)	_ _				V4, G
H014 (7CAI, 7CAC, 7CAB, 7CAK,7CAH)	Y369, A372,S373, F374, S375, T376, F377, K378, C379, Y380, V382, S383, P384, T385, D405, V407, R408, A411, P412, Q414, N437, V503	(Lv et al., 2020)	_ _
CR3022 (6ZLR, 6W7Y, 6ZH9, 7A5R)	Y369, N367, S375, T376, F377, K378, C379, Y380, G381, V382, S383, P384, T385, K386, L390, D428, F429, T430, F515, L517	(Huo et al., 2020)	_ _

Protein Data Bank (PDB) ids of the structure resolved for S (or RBD) complexed with the indicated antibodies (top 12 rows) or with the human receptor ACE (second row from bottom) or furin (bottom row). Of the 12 listed Abs, 10 are mAbs and 2 (H11-H4 and Nb20) are nanobodies (Nbs). Those mAbs and/or Nbs written in boldface contain mutated residues at their interface with S. For furin (last row), the Zenodo id for the model complexed with S is given (https://zenodo.org/record/4667694#.YVTJTprMJPZ).

Residues on S protein making close (<4.0 Å atom-atom) contacts with the mAbs, Nbs, ACE2, or furin. Mutation sites observed in the investigated variants are in boldface.

Reference for structural data.

PRODIGY energies Boltzmann-averaged over 63 conformations corresponding to top three ANM modes, Max RMSD = 5 Å, using stringent Energy Minimization (OpenMM tolerance = 1 kJ/mol). In each case the arrows (up/down) indicate the increase/decrease in neutralization. The cases with no significant effect are indicated by dashes in parentheses.

References to prior studies (experiments, except for those indicated by an asterisk, which are computational) and/or corresponding qualitative assessments indicated by arrows (up: increased neutralization; down; decreased neutralization; (−) no effect).

Three classifications are provided based on our observations first (excluding the NTD) and structural studies (Barnes et al., 2020a) and competition epitope mapping (Liu et al., 2020).

Results and discussion

Interactions of the Alpha, Beta, Gamma, UK2, and Delta variants with ACE2 and Furin

N501Y shared by Alpha, UK2, Beta, Gamma, and Omicron variants induces an increase in ACE2 binding

N501Y was the hallmark of the Alpha variant, shared with UK2, Beta, and Gamma variants. To examine the effect of N501Y in these contexts, we generated structural models for three mutant RBDs complexed with ACE2: N501Y (Alpha); double mutant E484K + N501Y (UK2; previous variant of concern VOC-202102–02); and K417N, E484K, and N501Y (Beta) using the approach described in STAR Methods—and performed a series of molecular dynamics (MD) simulations of the WT RBD-ACE2 complex (run 1 in Table S1) and these mutants (runs 2, 3, and 5).

We first verified that the Alpha RBD-ACE2 model (Cheng et al., 2021a) we constructed by in silico mutagenesis shows a root-mean-square deviation (RMSD) of only 0.9 Å from the more recently resolved cryo-EM structure (PDB: 7MJN) (Zhu et al., 2021) (Figure 1A), and our in silico Beta RBD-ACE2 model (Cheng et al., 2021a) aligned against the cryo-EM structure of Gamma RBD-ACE2 (PDB: 7NXC) (Dejnirattisai et al., 2021b) shows an RMSD of 1.0 Å only (Figure 1B), in support of the models used in simulations.

Figure 1

Change in interactions with ACE2 between WT RBD and Southern African, Beta variant 501.V2

(A) Alignment of computed Alpha RBD-ACE2 model (tan) with the cryo-EM resolved Alpha RBD-ACE2 structure (green: PDB: 7MJN) (Zhu et al., 2021). Right panel compares the side chain orientations in the model (thin sticks) and cryo-EM structure (thick sticks).

(B) Alignment of the Beta RBD-ACE2 model (tan) with the cryo-EM resolved Gamma RBD-ACE2 (green; PDB: 7NXC) (Dejnirattisai et al., 2021b). Positions of the three mutated amino acids (K417, N501, and E484 in the WT RBD) are displayed in yellow (model) and blue (cryo-EM) van der Waals (vdW) spheres.

(C and D) Comparison of the interactions associated with K417, E484, and N501 in the WT RBD (C) with those in the Beta variant triple mutant K417N, D501Y and E484K (D). The variant with substitution N501Y enables more contacts, resulting in higher binding affinity compared with WT. Two potential salt bridges (with ACE2 D30 and K31) formed by K417 and E484, respectively, in the WT (C) are broken in the Beta variant. Instead, the K484-E75 salt bridge intermittently forms. Here and in all ACE2 (and all nonspike) residues in italic.

(E–H) Detailed time evolution of relevant interfacial salt bridges (measured by intermolecular O-N distances at 0.2 ns intervals) formed between ACE2 and WT S (E–F) and the Beta variant (G–H). The salt bridge K417-D30 is persistently maintained in the WT, whereas E484-K31 shows very low (0.06) occupancy. In contrast, K484 forms a salt bridge with E75 (H) and thus promotes binding to ACE2. Results are based on three independent runs of 200 ns each carried for each of the complexes (runs 1 and 3; see Table S1). The mean values and standard deviations showed in F and H were calculated based on three runs.

Simulations of the WT RBD-ACE2 complex (PDB: 6LZG (Wang et al., 2020a); run 1) showed that N501 forms hydrogen bonds with ACE2 residues Y41 and K353 (Figure 1C; nonspike residues written in italic in the text and figures to distinguish them from spike residues). During the entire course of the three mutant sets of simulations, N501Y formed stable cation-π interactions with ACE2 K353 and engaged in aromatic-aromatic interactions with Y41 (Figure 1D). These associations are also captured by cryo-EM structural studies (Dejnirattisai et al., 2021b; Zhu et al., 2021). The same ACE2 residues have also been reported to regulate binding to SARS-CoV S with K353H, K353A, K353D, and Y41A disrupting interactions (Li et al., 2005).

Binding free energy histograms for the complexes formed between ACE2 and the RBD (mutants and WT) evaluated using the Molecular Mechanics/Generalized Born Surface Area (MM/GBSA) method (Genheden and Ryde, 2015) (Figure 2A) applied to MD snapshots showed that van der Waals interactions played a dominant role in defining the significantly higher stability of the Alpha variant RBD-ACE2 complex compared with WT RBD-ACE2. This effect persists in UK2, being further strengthened by electrostatic and solvent-accessible surface area (SASA) effects (Figure 2B; see also Figures S1A′-C′). Independent computations performed using PRODIGY (Xue et al., 2016) for evaluating the binding free energy of the MD snapshots (see STAR methods) confirmed the higher affinity of all Alpha, Beta, and UK2 variants for binding ACE2 compared with that of the WT RBD (Figures S1A–S1C), UK2 yielding the strongest binding affinity. Therefore, shared substitution N501Y promotes ACE2 binding to potentially increase infectivity. This is also indicated by binding assays (Cai et al., 2021; Supasa et al., 2021; Zhu et al., 2021), deep mutational scanning (Starr et al., 2020), and free energy perturbations (Fratev, 2021).

Figure 2

Histograms of different types of interaction energies contributing to the binding of the complexes formed between ACE2 and the WT RBD (blue bars) or (A) Alpha, (B) UK2, and (C) Delta variant RBD (orange bars)

Results are obtained by applying the MM/GBSA method to 800 evenly collected snapshots between 20 ns and 100 ns from each trajectory, generated in triplicates for each complex (runs 1–2 and 4–5; see Table S1). Contributions from electrostatic and solvation energies and vdW interactions are shown in the top and bottom panels, respectively. In each case, the change in binding energy ΔΔG (relative to the WT) contributed by the indicated subgroup of interactions is shown.

K417N mutation abolishes a salt bridge otherwise formed between ACE2 and the RBD, whereas E484K mutation favors a new interfacial salt bridge

K417, mutated to N and T in the Beta and Gamma variants, is specific to SARS-CoV-2 S, compared with SARS-CoV. The salt bridge between K417 (RBD) and D30 (ACE2) (Figure 1C) assists in RBD-ACE2 binding (Wang et al., 2020b). Simulations showed that this salt bridge is highly stable (Figure 1E), and its occurrence probability (also called salt bridge (SB) occupancy) of 0.59 ± 0.2 (averaged over three runs) is substantially higher than other (intermittently formed) interfacial salt bridges (Figure 1F). Interestingly, the local conformational changes induced by the mutation N501Y further stabilize the K417-D30 salt bridge (to increase SB occupancy to 0.81 ± 0.1; Figures S2A and S2B). The loss of this salt bridge by the mutations K417N and K417T would therefore decrease the binding affinity of the respective Beta and Gamma variants to ACE2 (compare Figures S2C and S2D). This has been observed for the Delta+ variants AY.1 and AY.2 where the K417N mutation reduced ACE2 binding affinity by at least 3-fold relative to the Delta variant (McCallum et al., 2021b).

The effect of E484K (shared by UK2, Beta, and Gamma; Table 1) is more complex. An interfacial salt bridge between E484 and ACE2 K31 was observed to form intermittently in the complex with the WT spike, albeit with low occupancy (Figures 1E and 1F). MD simulations with the variants (runs 3, 5 and 6 in Table S1) showed that a new interfacial salt bridge, K484-E75, forms instead for relatively extended durations (up to 35% of the trajectory in the Beta variant; Figures 1G and 1H). This salt bridge was also noted in a recent MD simulation study of E484K-associated mutants (Nelson et al., 2021). Simulations also revealed that the UK2 variant could form up to three interfacial salt bridges (K417-D30, K484-E75, and K458-E23; Figures S2A–S2C).

Notably, E484 is in a loop region that undergoes significant conformational fluctuations (Figures S2C and S2D), as also noted earlier (Bhattarai et al., 2021). This loop contains a potential disulfide bridge, C480-C488. Simulations performed with this disulfide (Table S1) showed a decrease in RMSD from the crystal structure (PDB: 6LZG) from 2.8 ± 1.0 Å in the absence of Cys-Cys link (runs 1 to 6) to 2.0 ± 0.5 Å in its presence (runs 7–11) due to suppression of loop mobility. Concomitantly, the occurrence of the salt bridge K484-E75 was prevented, drawing attention to the potential role of this Cys-Cys crosslink in the binding of ACE2. Yet, cumulative results from all runs confirmed that the ACE2 binding enthalpy ΔHMM/GBSA of the Beta variant (as well as the Alpha variant) remained significantly higher than that of the WT RBD (Figures S1E and S1F), consistent with deep mutational scanning (Starr et al., 2020) and in vitro evolution experiments (Zahradnik et al., 2021).

Delta RBD is distinguished by strongest association with ACE2 among all variants studied

The residues L452 and T478 (mutated in the Delta variant) do not make direct contacts with ACE2 in the resolved structures (see Table 2). Yet, given that simulations may induce conformational changes and expose newly formed interactions in the mutants not evident from examination of static structures or homology models (as observed in Figures 1C–1H), we carried out multiple MD runs and evaluated the energetics of conformers sampled in those simulations. MM/GBSA computations indicated the significantly higher affinity of the Delta RBD to bind ACE2 compared with other variants (Figure S1G), the major contribution again originating from van der Waals interactions (Figure 2C). Examination of the Delta RBD further revealed a change in the surface distribution of charges in comparison to the WT RBD, the two mutations T478K and L452R rendering it more basic (positively charged), which would overall increase its attraction to the predominantly acidic ACE2 interfacial surface (Figure 3A). A significant increase in binding affinity may partly explain the high virulence of the Delta variant.

Figure 3

L452R, T478K, and P681R in the Delta variant promote ACE2 and furin binding

(A) Electrostatic complementarity between RBD and ACE2 is strengthened in the Delta variant by substitution of positive charges at L452R and T478K.

(B) Molecular modeling of the Delta variant spike interactions with furin (cyan). The S subunit with the RBD in the open conformation is shown in magenta, and the two closed subunits are in green and light blue. The closeup view of the highly attractive interactions between the spike and furin shows the involvement of the new basic group (P681R) in strong interactions with the polyacidic binding site on furin. Basic residues are blue, acidic residues are red, and S368 (furin) and S686 (spike) are shown in yellow, highlighting the proximity of S368 to R685↑S686.

Overall evaluation of relative binding affinities of Alpha, Beta, Gamma, Delta, and UK2 variants to ACE2 and comparison with experimental data

In summary, out of the three mutations that define the Beta or Gamma variants, N501Y promotes the association with ACE2 through amino-aromatic and aromatic-aromatic interactions; E484K further contributes by electrostatic attractions that are otherwise “mute,” such that the double mutant E484K + N501Y (UK2) exhibits a higher affinity. Notably, the third mutation K417N or K417T in the Beta or Gamma variant reduces the ACE2-binding affinity and may help moderate the infectivity of these variants. In contrast, the mutations T478K and L452K significantly increase the affinity of the Delta variant RBD to bind ACE2.

Table 3 provides a summary of the computed and experimental dissociation constants Kd for the complexes that the WT, Alpha, UK2 (double mutant), Beta, and Delta RBDs form with ACE2. Notably, all variants exhibit higher affinity to bind ACE2 compared with WT RBD. The first three rows represent the computational results obtained here with two different methods (Methods I and II both use PRODIGY, the former applied to MD snapshots and the latter applied to HADDOCK-refined conformers [Cheng et al., 2021a]; see STAR methods) and their average. Method I yielded the histograms presented in Figures S1A–S1D. The experimentally measured dissociation constants show broad variations between different studies, e.g. the dissociation constant for WT RBD-ACE2 varies from 3.1 nM (Yang et al., 2021) to 218.29 nM (Gobeil et al., 2021). Therefore, it is hard to obtain a precise absolute value, but the relative values reported by the same group for different variants (listed in the same row) as well as the averages over the data compiled for the same variant indicate some consistent patterns. Overall, combining the data from independent computations and experiments, the rank order WT < Beta ≈ Gamma < Alpha ≈ Delta < UK2 emerges in increasing ACE2-binding affinity. The Delta and Alpha variants' affinities, as well as that of UK2, are stronger than that of the WT by a factor of approximately two based on our computations.

Table 3

Spike-ACE2 dissociation constants (Kd) from experiments and computations

	Kd (nM)					Ref
	WT	Alpha	UK2	Beta	Delta
Method Ia	23.1 ± 1.8	15.7 ± 1.2	15.0 ± 5.3	11.1 ± 5.0	9.1 ± 2.5	Pres
Method IIb	1.75 ± 0.4	0.76 ± 0.8	0.54 ± 0.2	2.00 ± 0.3	1.5 ± 0.4	pres, (Cheng et al., 2021a)
Comp avg	6.34	3.45	2.84	4.71	3.69
Experiments	12.8 ± 0.6	2.7 ± 0.4		15.1 ± 0.8		(Cai et al., 2021)
	75.1	10.7		4.0		(Supasa et al., 2021), (Zhou et al., 2021)
	3.1	1.3				(Yang et al., 2021)
	27.5 ± 4.8	11.8 ± 0.8		(23.5)	(22)	(Kim et al., 2021)c
	8.58 ± 0.4	6.25 ± 0.1				(Zhu et al., 2021)
	133	22		64		(Collier et al., 2021)
	17 ± 0.6	2.4 ± 0.4	1.40 ± 0.02	5.8 ± 0.8		(Laffeber et al., 2021)
	5.76	0.566				(Liu et al., 2021b)
	9.33	2.27	2.13	5.28		(Tanaka et al., 2021)
	8.3 ± 0.25	0.5 ± 0.01		0.5 ± 0.01		(Tian et al., 2021)
	23.9	2.26		7.81		(Bayarri-Olmos et al., 2021)
	62.6 ± 7.7	5.5 ± 2.4	3.7 ± 2.7	17.4 ± 3.1		(Barton et al., 2021)d
	218.29	45.20		72.63		(Gobeil et al., 2021)e
Exp Average	21.26	3.93	2.23	9.42

Based on 3 × 800 evenly collected MD snapshots taken from triplicate runs 1–4 in Table S1.

Based on HADDOCK refinement and energy minimization of ensemble of structural models.

Values from Figure 4 in the indicated reference are written in parentheses (not included in the average on the last line).

Kd measured by equilibrium binding.

Kd based on all mutations in variant.

P681R renders the polybasic S1/S2 cleavage site more attractive to furin, thus promoting proteolytic cleavage and fusogenicity

Proteolytic cleavage of the spike into S1 and S2 subunits at the site R685↑S686 is an essential step for the acquisition of fusogenic properties by the S2 trimer. Recognition by host cell proteases such as TMPRSS2 and furin is enabled by a polybasic stretch, R682RAR685, that serves as furin-recognition motif. The third mutation characterizing the Delta variant, P681R, is adjacent to this motif. The same residue is also mutated in the Alpha, UK2, and Omicron variants as P681H.

These substitutions render the segment P681RRAR685↑S686 even more basic than the WT cleavage site and might be expected to be more attractive to the acidic epitope on furin or other proteases, thus increasing the probability (or population) of cleaved spikes. To test this hypothesis, we examined the propensity of the Delta variant to bind to furin, compared with WT SARS-CoV-2 S using a previously generated complex of the WT spike with furin (Cheng et al., 2021b). The majority of models generated by ClusPro for the furin-spike complex showed that R681 associates with the acidic residues of furin, complementing other basic residues in the spike polybasic motif. Figure 3B illustrates such a highly stable binding pose between furin and the Delta variant. Three arginines, R681, R682, and R683, form a network of salt bridges with acidic residues from furin to enable the insertion of the catalytic serine S368 (forming the catalytic triad of furin with H194 and D153) near the cleavage site. In particular, the spatial and sequential juxtaposition of furin acidic residues E257 and E258 to R681 is an additional attractor to the segment RRRAR685↑S686, which implies an increased probability of interaction with furin in favor of the ensuing transition of the S2 trimer into its fusogenic form and shedding of the S1 subunits to the extracellular region. Notably, a recent experimental study indicated that P681R in the Delta variant enhanced replication via increased S1/S2 cleavage (Liu et al., 2021c).

It is worth noting the location of the mutation P681R overlaps with a superantigenic region (Cheng et al., 2020; Noval Rivas et al., 2020; Porritt et al., 2021) that has been shown to elicit multi-inflammatory response. The increased S1 shedding and exposure in the Delta variant may enhance such superantigenic effects. Notably, a recent study has shown that increased S1 levels in the extracellular medium correlate with more severe disease (Ogata et al., 2020).

N501 and K417 are essential residues that mediate the global dynamics of the RBD-ACE2 complex

Apart from interfacial interactions, it is of interest to examine the role of mutated residues in the structural mechanics of the spike and the complexes it forms. Gaussian Network Model (GNM) analysis (Bahar et al., 1997) of the global dynamics of the RBD-ACE2 complex revealed that N501 plays a central mechanical role as well as T500 and Y505, participating in a hinge center (minimum in the global mobility profile) as illustrated in Figures 4A and 4D. Further analysis of essential sites based on ESSA (Kaynak et al., 2020) yielded the profiles in Figures 4B and 4C for the RBD-ACE2 complex and isolated RBD, respectively. Here the peaks correspond to key sites/residues that play a dominant role in altering the essential/global dynamics of the structure (see STAR Methods). We present in panels F and G color-coded ribbon diagrams (peaks are colored red and minima, blue), which show that the residues at the interface between ACE2 and RBD as well as those lying in the ACE2 catalytic cleft are essential residues in the global dynamics of the complex.

Figure 4

Essential sites dominating the structural dynamics of RBD complexed with ACE2

(A) Mean-square fluctuations of residues driven by three GNM softest modes. Mutation sites occupying minima (hinge sites) are labeled.

(B) Results from ESSA of RBD-ACE2 complex. K417 and N501, and associated ACE2 D30 and K31, are distinguished as essential residues (peaks), as is L452. Three RBD residues implicated in the human-to-animal passage of the virus, also detected as essential residues, are labeled in cyan.

(C) ESSA profile for isolated RBD (red curve), compared with its counterpart in the complex (black), shows that K417, N501, and Y505 acquire a significant mechanical role upon complexation, whereas L452 retains its essential role.

(D) Position of N501 in the RBD-ACE2 complex. The RBD is color-coded by GNM mobility profile (from red, most constrained, to cyan, most mobile).

(E) Closeup view of the interfacial regions displaying an allosteric pocket predicted by ESSA (red wire) and key residues in the vicinity.

(F and G) Color-coded diagrams illustrating the residues distinguished by ESSA to serve as essential residues (peaks in the ESSA profile for the complex in panel c; colored red in the diagram). Two regions, the catalytic cleft of the ACE2 receptor and the interface between ACE2 and RBD, stand out as essential sites.

In Figure 4B, our analysis highlights two of the residues mutated in the Beta variant (K417 and N501; in red) and two ACE2 residues (Y41 and K353; in blue) reported (Li et al., 2005) to regulate ACE2-spike binding, as essential residues. It is interesting to note that several residues gain functional importance (emerge as peaks) upon complexation (Figure 4C), including K417 and N501 as well as T500 and Y505 in the RBD and D30, K31, Y41, and K353 in ACE2. Moreover, the mutation site L452 shared by Delta and Epsilon variants stands out as an essential site in the RBD. In contrast, the Delta mutation site T478 and Beta mutation site E484 do not exert a substantial effect on the global dynamics.

We also note that in an insightful study of the spike interactions with ACE2 and Regeneron mAbs, the authors predicted N501 and K417 (as well as E406 and N439) as effectors of allosteric communication (Verkhivker et al., 2021). Our ESSA-based allosteric pocket prediction corroborates that this interfacial region (comprising RBD K417, N501, and Y505 and ACE2 D30 and K31; Figure 4E) has an allosteric potential to interfere with the overall conformational dynamics of RBD-containing complexes.

Interactions of variants with monoclonal antibodies and nanobodies

We analyzed a set of 12 Abs, consisting of ten mAbs and two Nbs H11-H4 and Nb20, whose complexes with SARS-CoV-2 S or RBD have been structurally resolved (Table 2, columns 1 and 3). Figure S3A displays their binding poses, generated by structurally aligning the RBD regions of the corresponding S/RBD–Ab complex structures. Table 2 column 2 lists the spike residues that make interfacial contacts with the listed Abs, which reveals five epitope classes, termed here classes I to V (last column). Examination of these residues shows that none of the Abs makes contacts with the residues T478, D614, and P681 mutated in one or more of the examined variants (Table 1). Likewise, the binding epitopes of four mAbs (REGN10987, EY6A, HO14, and CR3022) do not include the mutation sites K417, E484, N501, L452, and T478 (Figure S3B). These four mAbs belong to class 3 (REGN10987) or 4 in the structural nomenclature (Barnes et al., 2020a) (corresponding to classes H and G in the competition epitope mapping) (Liu et al., 2020), which are designated as IV and V here. We focused therefore on the remaining eight Abs, and quantitatively examined the effects on the binding to Beta and Delta variants. The computed changes in binding free energies (see STAR methods) and corresponding data from the literature (experiments and/or prior computations) are listed in columns 4–8, with the up/down arrows in green/red indicating an increase/decrease in binding affinity, compared with the WT RBD. Table S2 lists the computed binding affinities to WT RBD.

The mutations K417N and E484K reduce the binding affinity of REGN10933, C105, B38, LY-CoV555, and H11-H4 against the Beta variant, while increasing that of BD23

Two mAbs adversely impacted by the mutations in the variants are C105 (Figures 5A–5C) and REGN10933 (Figures 5D–5I), which both belong to Class 1 or E and depend heavily on K417. In C105, K417 coordinates a polar network that supports its high affinity to bind the SARS-CoV-2 RBD. K417 forms a salt bridge with C105 residues E96/E99, which also engage RBD Y453 (blue circle in Figure 5A; zoomed-in view in Figure 5B). This network is completely abolished by the substitution K417N (Figure 5C) in the Beta variant. This significant loss in the attractive interaction is partially alleviated by a new favorable interaction between N417 and mAb Y52, and a cation-π interaction that Y501 makes with mAb K31, but the net effect is a reduction in binding affinity by 0.42 kcal/mol relative to the WT RBD (Table 2, column 4). The Gamma variant with the shorter residue T417 cannot form compensatory interactions (not shown), suggesting an even greater escape from C105 and other Abs interacting with this site and potentially explaining its prolonged dominance in Brazil, Chile, and Argentina.

Figure 5

Disruption of salt bridges formed by K417 weakens the association of the mAbs C105 and REGN10933 with Beta and Delta RBD, compared with WT RBD

(A–I) The panels compare the interactions of C105 (A–C) and REGN10933 (D–I) with WT RBD (A, B, D, and G), the Beta variant (C and H), and the Delta variant (E, F, and I). Central salt bridges (blue circles) with C105 (K417-E99/E96) and REGN10933 (K417-D31) are lost due to the substitution K417N as shown for C105 (B and C). Another salt bridge, R403-D31, is formed stably in (D) and (H). Both interactions are disrupted in the Delta variant (E and I) through rearrangement of this region by the L452R mutation. N501 makes few interfacial contacts with either mAb, but is in the vicinity of C105 K31, with which N501Y can form a new cation-π interaction. (G–I) Time evolution of O-N distances of potential salt bridges between REGN10933 and WT RBD (G), the Beta variant (H), and the Delta variant (I). Salt bridges form once the distance is less than 4 Å. MD simulations reveal additional interfacial interactions of REGN10933 including E484-R100 and R403-D31 salt bridges revealed by multiple MD simulation runs. A new salt bridge K478-D92 stably formed in the complex with the Delta variant upon rearrangement of the region around E484 and T478K (green circle in D) is shown in (F) and (I). Results are based on runs 12–14 (see Table S1).

The effect of mutations on the interaction with the Regeneron mAb REGN10933 (Figures 5D–5I) is more complex. In the experimentally resolved wt RBD-REGN10933 complex, K417 forms a salt bridge with mAb D31 (as well as engaging T102), which in the Beta variant is replaced by weaker polar interactions. This implies a weakening in the neutralizing efficacy of REGN10933. MD simulations reveal two additional salt bridges, R403-D31 and E484-R100, that stabilize interactions with the WT (Figures 5D and 5G). This is accompanied by rearrangement of the loop containing E484 (and T478) as noted in previous simulations (Bhattarai et al., 2021). Even though R403-D31 partly compensates for loss of K417-D31 in the Beta variant (Figure 5H), the additional loss of the E484-R100 salt bridge and accompanying changes in interactions overall decrease the binding efficacy of REGN10933 to the Beta variant by 1.93 kcal/mol.

For Eli Lilly LY-CoV555, E484 makes two critical salt bridges with mAb residues R50 and R96 (Figure 6A, red circle). E484K severely disrupts these interactions via charge repulsion. Indeed, during our MD simulations, we observed a dissociation of mAb LY-CoV555 from its initial contact pose with the Beta variant RBD within 10 ns, succeeded by occasional return to bound forms, and eventual drift away after ∼150 ns (Figure 6B). The significant reduction in MM/GBSA interaction energy (Figure 6C) is also consistent with the escape of Beta RBD from mAb LY-CoV555. Independent computations using anisotropic network model (ANM) (Eyal et al., 2006) conformers (see STAR Methods) also confirmed a reduction in binding affinity of 0.29 kcal/mol (Table 2).

Figure 6

Loss of salt bridges and hydrophobic contacts disrupts the association of Eli Lilly LY-CoV555 with Beta, Gamma, Delta, and Epsilon variants

(A) The interactions of LY-CoV555 with WT RBD. E484 forms salt bridges with R50/R96 from LY-CoV555; L452 makes hydrophobic contacts with I54/L55 from LY-CoV555.

(B and C) Time evolution of the RMSDs of LY-CoV555 from its initial conformer (B) and binding free energies (ΔHMM-GBSA; C) in the WT (blue) and Beta variant (green). Notably, E484K in Beta (or Gamma) variant weakens interfacial interactions and LY-CoV555 drifts away from the Beta variant with significant reduction of binding free energy.

(D) Interactions of the LY-CoV555 with Epsilon or Delta variant RBDs with the L452R mutation, which reduces hydrophobic interactions associated with L452.

Figure S4 panels A–B present the results for the Nb H11-H4. The strong attractive interaction between WT RBD E484 and H11-H4 R52 turns into a strong repulsion between R52 and K484, and the free energy change ΔΔGbinding of the Beta variant with respect to WT confirms a weakening in interaction by 0.32 kcal/mol. Although a similar effect is observed in BD23 (panels C–D) where the original favorable interaction between E484 and BD23 R107 is disrupted in the Beta variant, compensating interactions take place near K484, plus Y501 of the Beta RBD engages in close interactions with BD23, such that this mAb exhibits an increase in binding affinity to Beta variant (ΔΔGbinding = −0.94 kcal/mol relative to WT RBD; Table 2), in accord with recent computations (Ray et al., 2021).

In the case of the second Eli Lilly mAb, LY-CoV016, a salt bridge K417-D104 remained stably formed during the entire duration of simulations of the WT complex. The mutation K417N abolished this contact. Hydrogen bonds between N417 and Y33/Y52 partially compensated for the loss of the salt bridge in the Beta variant. PRODIGY calculations applied to the ensemble of conformers generated by MD simulations indicated a 2-fold decrease in binding affinity. On the other hand, ANM sampling of conformers identified favorable interaction poses enabled by collective motions of the complex with the Beta variant (see STAR methods), leading to an increase in binding affinity by 1.86 kcal/mol. It is possible that these “optimized” conformers involving global rearrangements in the structure are not accessible during simulations of 100s of nanoseconds. Experiments, on the other hand, indicate either a decrease in neutralization efficacy (Arora et al., 2021; Wang et al., 2021b) or no effect (Chen et al., 2021a).

The Abs REGN10933, B38, BD23, LY-CoV555, and H11-H4 undergo a decrease in binding affinity to Delta RBD compared with the WT RBD

Similar analysis conducted for the Delta variant led to a reduction in binding affinity of many Abs (Table 2, column 5) consistent with a recent study where 16 out of 26 mAbs showed a marked decrease or complete loss in neutralizing activity to the Delta variant (Mlcochova et al., 2021).

In the case of the complex formed with REGN10933, the L452R mutation is observed to induce a conformational rearrangement leading to the disruption of the K417-D31 and R403-D31 salt bridges (Figures 5E and 5I). The T478K mutation generates a new salt bridge K478-D92 (Figure 5F), which may help in retaining the binding and neutralizing activity of mAb REGN10933 to the Delta variant observed in an extensive study (Mlcochova et al., 2021). On the other hand, a decrease in neutralization activity has been reported in another study (Tada et al., 2021). Our current evaluation indicates a 6-fold decrease in the binding affinity of REGN10933 to Delta variant (ΔΔGbinding = 1.06 kcal/mol; Table 2).

The RBD residue L452 is also engaged in a network of hydrophobic contacts with LY-CoV0555 I54 and L55 in WT S (complemented by L492) (Figure 6A; blue circle). The L452R mutation breaks these favorable hydrophobic interactions (Figure 6D). Even though the Delta variant RBD remained bound to LY-CoV555 due to the salt bridge between E484 and R50/R96 in our 200 ns MD simulations, energy calculations showed a 3-fold reduction in binding affinity of LY-CoV555 due to L452R mutation shared by the Epsilon and Delta variants in line with recent experiments showing reduced neutralizing activity by this mAb (Mlcochova et al., 2021). Notably, a previous experimental study also found that L452R enabled escape from LY-CoV0555 (Starr et al., 2021b). The Kappa variant (B.1.617.1) containing E484Q and L452R (Cherian et al., 2021; Verghese et al., 2021) is less able to escape from Abs (Ferreira et al., 2021; McCallum et al., 2021b), which is likely due to Q484 still being able to make favorable interfacial contacts. Computations repeated for the second Eli Lilly mAb, LY-CoV016, on the other hand, indicated no clear effect on binding of this mAb to Delta variant (Table 2), in agreement with its experimentally observed retained neutralization activity (Mlcochova et al., 2021).

Overall, the results compiled in Table 2 columns 4–5 are in agreement with experiments or prior simulations whenever there are available data. A few cases stand out where the predicted changes in binding affinities differ from previous observations. For example, in the case of BD23, we predict a decrease in binding affinity (by 0.33 kcal/mol) to the Delta variant compared with the WT RBD. This contrasts the favorable interaction observed in a recent study (Ray et al., 2021). Figures S4C and S4D show that the mAb BD23 harbors the hydrophobic residue W105 in the vicinity of L452. This raises the possibility that the W105-R452 cation-π interaction may strengthen the interactions of BD23 with the Delta and Epsilon variants.

Computations and experiments indicate that Nb20 retains its neutralizing efficacy with respect to the Delta variant

Finally, as a proof of concept, we explored the validity of our computational predictions with the help of enzyme-linked immunosorbent assays (ELISA) performed for Nb20, an ultrapotent and inhalable nanobody to block WT SARS-CoV-2 (Xiang et al., 2020), for the WT RBD and a mutant RBD with critical substitutions L452R and T478K of the Delta variant.

In the case of the WT RBD, Nb20 forms a salt bridge E484-R31/R97 and hydrophobic interaction L452-M55 (Figure 7A), both of which remained stable during MD simulations, the E484-R31/R97 salt bridge with an occupancy of 96 ± 4%, suggesting these may play a substantial role in stabilizing the Nb20 complex.

Figure 7

Computational and experimental assessment of the nanobody Nb20 binding to the WT and Delta RBDs

(A) Crystal structure of SARS-CoV-2 RBD with Nb20 (PDB: 7JVB) (Xiang et al., 2020). A key interaction is the salt bridge between E484 and R31 (Nb20).

(B) MD-resolved Delta variant RBD with Nb20.

(C) Comparison of the salt bridge occupancy in the WT (blue) and Delta variant (orange), showing sustained occupancy of E484-R31 in both. The mean values and standard deviations were calculated based on three runs.

(D) ELISA binding of Nb20 to WT RBD and Delta variant. The experiment was repeated three times, and data points are shown as mean ± SD The efficacy of Nb20 to the Delta variant is thus expected to be unchanged compared with WT spike.

For the Delta variant, our modeling suggests that the hydrophobic interaction (L452-M55 in the WT) would not significantly contribute to the overall RBD binding energy as a result of the L452R substitution. However, L452R promoted a new interaction with S52 via hydrogen bonding (Figure 7B). Interestingly, the T478K mutation favored a new salt bridge formation with E42/E44 (Figure 7C). The E484-associated salt bridges remained stable (>95% occupancy), indicating that E484 was essential to binding Nb20. Mutations breaking this salt bridge, e.g. E484K of Beta variant, reduce, if not abolish, its ultrahigh-affinity to the RBD (Sun et al., 2021). The predicted change in interaction energy of −0.15 kcal/mol (Table 2) indeed suggests that the Delta variant would not evade neutralization by Nb20. In agreement with this prediction, Nb20 is observed in binding assays to retain (if not increase) its potency for the Delta RBD (Figure 7D).

Conclusion

In silico models and simulations shed light onto the molecular basis of the altered binding and neutralization properties of the variants

Recently observed variants continue to adapt to binding the host receptor ACE2, which together with the high mutation rate of SARS-CoV-2, points to a risk of still-undetected additional members of the coronaviridae family making the jump from their traditional hosts to humans in the future. Our analysis shows an enhancement in the binding affinity to all variants compared with the WT RBD (Table 3). Computed dissociation constants fall within the broad range observed in experiments. The relative ACE2-binding affinities of variants is indicated by both computations and experiments to yield a rank order of WT < Beta ∼ Gamma < Alpha < Delta < UK2 with increasing affinity.

To assess the impact of the variants on interaction with and hence neutralization by antibodies, we examined 10 mAbs and 2 Nbs (Table 2). We found five Ab binding epitopes (I to V) on the SARS-CoV-2 RBD similar to those identified in other studies (Barnes et al., 2020a; Liu et al., 2020; Xiang et al., 2020), including major classes (I to III) where residues E484 and K417 contribute to the binding and their mutations could be responsible for evading selected Abs. The interactions between either the Beta or Delta spike and many Abs such as C105, B38, REGN10933, LY-CoV555, and H11-H4 were shown to be weakened if not abolished. MD simulations of hundreds of nanoseconds revealed associations and dissociations that stabilized or destabilized the specific complexes, leading to the final assessments. Notably, many structure-based arguments we recently made with regard to the infectivity and neutralization of the Alpha and Beta variants in our recent report (Cheng et al., 2021a) are corroborated by structural data in a recent study (Cai et al., 2021), and simulations may further help provide efficient assessments of the effect of future mutations.

A recent mapping of RBD mutations that escape REGN-COV2 cocktail and Eli Lilly’s LY-CoV016 antibodies (Starr et al., 2021a) using a deep mutational scanning method (Greaney et al., 2021c) revealed several sites whose point mutations mediate escape. Interestingly, the mutation (E406W) escaped the cocktail of both REGN antibodies even though it does not make direct contacts with either mAb in the resolved structure (Table 2), inviting attention to the necessity of considering structural changes induced by substitutions, hence the need to conduct MD simulations as well as network-based analysis of allosteric dynamics to identify critical substitutions. Our analysis indicates the critical mechanical roles of N501, Y505, and E406 (Figure 4).

A substitution at the polybasic cleavage site with superantigen-like properties may contribute to the transmissibility and infectivity of the Alpha and Delta variants

The location of the Delta variant spike mutation P681R is worth special attention. Our simulations demonstrated that the mutation rendered this site more attractive (than the WT) to furin (Figure 3B), which could thus enhance proteolytic cleavage and thereby viral entry. Another important aspect of this specific site P681RRAR685 is that it belongs to a putative superantigenic segment having sequence and structure features highly similar to those of a bacterial toxin, Staphyloccocal enterotoxin B (SEB), potentially causing a cytokine storm or MIS-C (Cheng et al., 2020). Enhanced cleavage of S and shedding of the S1 subunits would further expose the highly reactive C-terminal end R681RRAR685 of the superantigenic segment in the Delta variant, which might further escalate hyperinflammatory responses. Indeed, recent studies show that higher serum levels of circulating free S1 spike protein are associated with more severe disease among adult COVID-19 patients (Ogata et al., 2020), as well as among children with severe MIS-C disease (Porritt et al., 2021). The importance of the furin cleavage site in SARS-CoV-2 infection and pathogenesis is confirmed in several studies (Johnson et al., 2021; Lau et al., 2020).

Given the importance of the furin cleavage site and the shed S1 in infectivity and disease pathology, it is of great interest that two dominant SARS-CoV-2 variants, Alpha and Delta, both have altered furin cleavage sites. In the Alpha variant, which was 50% more contagious than the original Wuhan strain, P681 is substituted by a histidine (P681H); in the Delta variant, which is about 50% more contagious than the Alpha variant, it is substituted by an arginine (P681R) (SPI-M-O, 2021). Both changes make the sequence more basic, potentially leading to a more effective recognition by the polyacidic furin motif. Our simulations have shown that the P681R form binds more efficiently to furin, which would result in increased S1/S2 cleavage as well as cleaved products. Notably, a recent study reported that Delta variant (P681R) S proteins were highly cleaved (>50%; Peacock et al., 2021). A more recent preprint has also shown that P681R mutation enhances the fitness of Delta over the Alpha variant and that increased cleavage leads to increased infection via cell surface entry (Liu et al., 2021c). Finally, because circulating levels of cleaved S1 correlate with disease severity, our simulations are thus in line with reports that the Delta variant may lead to more severe disease (Fisman and Tuite, 2021; Sheikh et al., 2021). Our data, together with other recent studies, suggest that spike mutations that affect the furin cleavage site should be carefully monitored, as they may increase viral fitness, infectivity, and virulence.

A loss in mAb efficacy may be alleviated by combination treatments consisting of mAbs that bind alternative, non-overlapping epitopes

Overall, our analysis indicates the ability of variants to evade mAb or Nb treatment depends on the mutations and on the specific Abs and may be alleviated by combination treatments consisting of mAbs that bind non-overlapping epitopes. Notably, the mutation sites on the Delta variant make no contacts with most of the investigated mAbs, except the mAb LY-CoV555 and the Nbs H11-H4 and Nb20. So, this variant would not be expected to escape all the mAbs, except for these three, which we examined more carefully. Our analysis demonstrated that among them Nb20 would not lose its potency, which was also confirmed by ELISA experiments (Figure 7); whereas LY-CoV555 would exhibit a significant weakening in its neutralization activity.

A potential approach to minimize the impact of mAb escape mutations is to develop additional mAbs with epitopes that are outside of the RBD and/or cross-reactive across coronaviridae, with low tolerance for mutation (Pinto et al., 2020; Wec et al., 2020). We have recently described such a mAb (Cheng et al., 2021b) that binds the putative superantigenic-like motif of SARS-CoV-2 S (Cheng et al., 2020). That particular mAb, 6D3, presents the dual advantage of potentially blocking the S1/S2 cleavage site, thus interfering with proteolytic cleavage that is essential for viral entry, in addition to targeting the superantigenic region that may contribute to MIS in adults, children, and adolescents with severe COVID-19 (Cheng et al., 2021b). Given the potential emergence of immune evasion mutations that maintain virulence or those that confer resistance to immunizations or therapies (Weisblum et al., 2020), in silico assisted genomic/molecular surveillance may provide valuable feedback for accelerating the design of experimental studies in response to the pandemic.

Limitations of the study

In the present study, we have benefited from the wealth of structural data resolved for several complexes of the SARS-CoV-2 spike (or RBD) with ACE2 or various Abs. Although the in silico structural models for mutants often yield good agreement with resolved structures (as illustrated here in Figures 1A and 1B), model selection and refinement often require extensive simulations (docking, energy minimization and/or MD simulations), which lead to multiple models and conformers, and accurate sampling of the complete space of conformations is often a challenge. Even though MD simulations are considered as the ground-truth in molecular computations, they are also known to suffer from sampling inaccuracies, in addition to approximations in the force fields and parameters. PRODIGY, optimized to correlate experimental binding affinities with the number of interfacial contacts and non-interacting surface properties (Kastritis et al., 2014) derived from resolved structures (Vangone and Bonvin, 2015), appears to be one of the best-performing empirical predictors reported to date (Vangone and Bonvin, 2015), which we adopted here, together with methods of different complexities, from HADDOCK refinement with short energy minimization to detailed MM/GBSA evaluation of thousands of conformers sampled in MD trajectories. Current simulations starting from equilibrium conformations stabilized under crystallization conditions may remain in the close vicinity of this highly stable state and not sample the whole conformational space accessible under physiological conditions. Finally, MM/GBSA energies evaluated here represent the enthalpic contribution to binding free energies. Examination of the differences in entropic effects between WT and variant complexes would be valuable. Despite all these limitations, the combination of an empirical machine-learning type approach (PRODIGY energies) with full-atomic free energy calculations (MM/GBSA) helped consolidate the data from simulations and make robust inferences on the relative affinity of different variants to bind ACE2 (Table 3) or on the impact of mutations on the neutralizing activity of various Abs (Table 2). The method of analysis adopted here may be applied to other variants to estimate their ACE2-binding affinity relative to that of the WT spike, as well as their propensity to escape specific Abs. Furthermore, the comparative analysis gives insights into the mechanistic bases of the observed behavior, while providing a framework for predictions on systems not experimentally characterized to date.

STAR★Methods

Key resources table

Resource availability

Lead contact

All requests should be directed to and will be fulfilled by the Lead Contact, Ivet Bahar (bahar@pitt.edu).

Materials availability

This study did not generate new unique reagents.

Method details

In silico mutagenesis

We performed in silico mutagenesis and assessment of structural models for SARS-CoV-2 S RBD variants complexed with ACE2 and Abs. For ACE2, we used the resolved structure of the complex between WT S RBD and ACE2 (PDB: 6LZG) (Wang et al., 2020a). We generated structural models for five mutant RBDs complexed with ACE2 - (1) N501Y (Alpha); (2) E484K; (3) N501Y and E484K; (4) K417N, E484K, and N501Y (Beta); and (5) L452R and T478K (Delta) - using CHARMM-GUI and MD energy minimization with the refinement protocol and default parameters implemented in the webserver HADDOCK 2.4 (van Zundert et al., 2016). Our in silico models yielded excellent agreement with cryo-EM data (Figures 1A and B). Interactions were visualized with PyMOL (DeLano, 2002) using mutagenesis and sculpting wizards for rotamer optimization. Unless otherwise specified, the RBD residues 333-527 and ACE2 residues 19-614 were included in the complexes. Similar modelling and computational protocols were adopted for the complexes with Abs (both heavy and light chains), using as template the PDB structures listed in Table S1.

MD simulation systems and protocols

All initial MD simulation systems were generated using CHARMM-GUI Solution Builder module (Jo et al., 2008). Simulations were performed in triplicates for the lowest binding affinity HADDOCK-refined WT or mutant RBD-ACE2 as described above (runs 1-6 in Table S1). We also performed additional runs (runs 7-11, in triplicate) to examine the effect of N-linked glycans and disulphides. Both groups had explicit water molecules to a distance of 10 Å away from protein edges, and sodium and chloride ions corresponding to 0.15 M NaCl, which resulted in a simulation box of 133×133×133 Å3 for the RBD-ACE2 complex. The boxes for RBD-Ab complexes (runs 12-23, each in duplicate) varied from 105 × 105 × 105 Å3 (RBD-Nb20) to 147 × 147 × 147 Å3 (RBD-LY CoV-555).

All MD simulations were performed using the NAMD (version 2.13) (Phillips et al., 2005), with a well-established protocol (Lee et al., 2016). The CHARMM36 force field with CMAP corrections was used for protein, water, and glycan molecules (Guvench et al., 2011; Huang et al., 2017). Briefly, prior to productive runs, each system was energy-minimized for 20,000 steps, followed by a modified Nosé–Hoover (Hoover, 1985; Nosé, 1984) constant pressure (1 bar) equilibration with temperature varied from 1 K to 310 K within 62,000 steps, during which the constraint on the protein backbone was set to 10 kcal/mol. A subsequent 4 ns Nosé–Hoover (Hoover, 1985; Nosé, 1984) constant pressure and temperature (1 bar, 310 K; NPT) simulation was performed, during which the constraints on the protein backbone were reduced from 10 to 0 kcal/mol. Finally, the unconstrained protein was subjected to NPT simulations. Periodic boundary conditions were employed for all simulations, and the particle mesh Ewald (PME) method (Essmann et al., 1995) was used for long-range electrostatic interactions with the pair list distance (pairlistdist) 13.5 Å. The simulation time step was set to 2 fs with the covalent bonds to hydrogen atoms constrained with the SHAKE algorithm (Ryckaert et al., 1977). The force-based switching function was used for the LJ interactions with switching distance (switchdist) set to 10 Å. We applied Langevin dynamics conditions with a piston period (langevinPistonPeriod) of 200 fs and a piston decay (langevinPistonDecay) of 100 fs as well as Langevin temperature coupling with a friction coefficient (langevinDamping) of 1ps-1. Table S1 lists all simulations, along with the corresponding systems, durations and PDB files and chains used in simulations. VMD (Humphrey et al., 1996) with in-house scripts was used for visualization and trajectory analysis.

MM/GBSA computational protocol and parameters

We adopted the steps described in earlier work (Tanner et al., 2011). First, we took our MD trajectories generated for the protein−protein complex, i.e. RBD-ACE2, in explicit solvent. Second, all solvent molecules, ions, and glycan ligands if any, were removed from each MD snapshot, yielding trajectories for the RBD, ACE2, and RBD-ACE2 complex. Third, for each of these trajectories, the MM/GBSA free energy is calculated using the GBIS module (Tanner et al., 2011) implemented in NAMD. The binding free energy change is estimated aswhere G, G, G are the MM/GBSA energies for the RBD-ACE2 complex, RBD and ACE2, respectively, calculated using (Genheden and Ryde, 2015):

The first three terms are the standard molecular mechanics (MM) energy terms from bonded (bond length, bond angle, and dihedral angle), electrostatic, and van der Waals (vdW) interactions. Gpol and Gnp are the polar and non-polar contributions to the solvation free energies. Gpol is obtained using the Generalized Born (GB) model, the non-polar term Gnp is estimated from a linear relation to the solvent accessible surface area (SASA) (Genheden and Ryde, 2015). Note that the MM/GBSA energy (POTENTIAL term from the output) released by NAMD does not contain the entropic contribution to free energy. Therefore, in a strict sense, the enthalpy change is evaluated, designated as in the manuscript. Binding is usually accompanied by an entropy loss, and the associated free energy increase is comparable in magnitude to the free energy decrease induced by enthalpic effects.

We used the default parameters set in the GBIS module (Tanner et al., 2011). Briefly, CHARMM36 force field with CMAP corrections (Huang et al., 2017) was used to calculate MM energy. For GBSA calculations, the dielectric constant was set 78.5 and ion concentration to 0.3 M. The surface tension was set to be 0.00542 kcal/(mol·Å2). The switching distance (switchdist) for the LJ interactions was 15 Å; the long-range electrostatic cutoff distance (pairlistdist) was 18 Å; and the cutoff (cutoff) for the non-bonded interactions was 16 Å.

Evaluation of binding energetics

We evaluated the binding energetics ACE2-RBD (WT or variant) complexes using three methods: (I) PRODIGY (Xue et al., 2016), for evaluating the binding energy (ΔG) based on a series of snapshot from MD trajectories; (II) PRODIGY, applied to an ensemble of structural models refined by HADDOCK and energy-minimized by short MD, as described earlier (Cheng et al., 2021a); and (III) the MM/GBSA method described in the previous section, for evaluating the enthalpic contribution ΔH to the binding free energies for the same series of MD snapshots as in (I). In methods I and III, ΔG or ΔH histograms for binding to ACE2 were generated based on 800 snapshots evenly collected during the MD simulation time interval 20 < t ≤ 100 ns for each run. Dissociation constants were calculated using K = exp(ΔG/RT) x 109 with RT = 0.6 kcal/mol at T = 300 K, yielding nM units.

The binding energies of the complexes formed by the Abs were evaluated by generating an ensemble of conformers along ANM modes. We considered the collective motions along the first three lowest frequency modes. A total of 21 conformations generated for each mode resulted in a total of 63 conformers for each complex. The softest (first) mode was allowed a maximum deformation (RMSD) of 5.0 Å and the next two modes were scaled proportionately to their eigenvalues. Each of the 63 conformations were stringently optimized using OpenMM (Eastman et al., 2017) for energy minimization with the Amber14 force field (Wang et al., 2004). Binding energies of resulting conformers were evaluated using PRODIGY.

Impact of mutations on the structural mechanics

We evaluated the global dynamics of the RBD-ACE2 complex using the Gaussian network model (GNM) (Bahar et al., 1997) and essential site scanning analysis (ESSA) (Kaynak et al., 2020). Both tools are accessible and described in detail in our open source ProDy API (Bakan et al., 2011, 2014; Zhang et al., 2021b). GNM analysis provides information on the role of individual residues in the global dynamics (Bahar et al., 2010), i.e. collective motions at the low frequency end of the mode spectrum accessible under equilibrium conditions. Briefly, we evaluated the GNM-predicted square fluctuations of residues as driven by three slowest modes of motions and identified the most constrained residues (i.e., minima, referred to as hinge sites). These sites cannot usually tolerate mutations (Haliloglu and Bahar, 2015). ESSA, on the other hand, scans all residues for identifying the so-called essential sites whose perturbation (by local crowding or ligand binding) would elicit a strong change in the frequency distribution of global motions. Integration of ESSA with a pocket detection algorithm predicts allosteric pockets (Kaynak et al., 2020).

Modelling of Delta variant spike-furin complex

The three mutations L452R, T478K, and P681R were generated using CHARMM-GUI (Jo et al., 2008) using the WT Spike model from our previous studies (Cheng et al., 2020, 2021b). A crystal structure of human furin (PDB: 5JMO) (Dahms et al., 2016) was used for docking simulations by ClusPro (Kozakov et al., 2017) with the same protocol as for WT S (Cheng et al., 2021b). Complexes were refined using HADDOCK 2.4.

ELISA (Enzyme-Linked Immunosorbent Assay)

WT SARS-CoV-2 RBD and Delta variant RBD were coated onto 96-well ELISA plates at 3 μg/ml in coating buffer (15 mM sodium carbonate, 35 mM sodium bicarbonate, pH 9.6) at 4°C overnight. The plates were decanted and washed with washing buffer (1x PBS, 0.05 % Tween 20) three times, and then blocked for 2 hours at room temperature in blocking buffer (1x PBS, 0.05 % Tween 20, 5 % milk powder). Nanobody Nb20 was 5x serially diluted in blocking buffers starting from 2.5 μM with 8 different concentrations and incubated with the plate for 2 hours at room temperature. After three washes, anti-T7 tag HRP-conjugated secondary antibodies (Thermo, cat# PA1-31449) were diluted at 1:7500 and incubated at room temperature for 1.5 hours. Upon washing, freshly prepared 3,3’,5,5’-Tetramethylbenzidine substrate was added to the plate for signal development in the dark for 10 minutes. The reaction was quenched with a STOP solution. The plates were measured at wavelengths of 450 nm with background subtraction at 550 nm. Raw data were processed and fitted to the 4PL curve using Graphpad Prism 9.0.

Quantification and statistical analysis

For each system, we performed multiple MD runs and verified the convergence and reproducibility of the trajectories. We evaluated the binding energetics using three methods, as described above. Free energy histograms were generated by using 800 x n snapshots generated for n = 2 or 3 runs performed for each system, and corresponding mean values and standard deviations are reported for each system. Statistical analyses were performed using Matlab. We repeated ELISA experiments three times and we show data points as mean +/- standard deviation.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Antibodies
Nb20	PMID: 33154108	N/A
Anti-T7 tag HRP-conjugated secondary antibodies	Thermo Fisher	Cat# PA1-31449
Chemicals, peptides, and recombinant proteins
ELISA kit	R&D system	DY990-DY999
SARS-CoV-2 RBD wild-type (his-tag)	Acro Biosystem	Cat# SPD-C52H3
SARS-CoV-2 RBD Delta (his-tag)	Acro Biosystem	Cat# SPD-C52Hh
Deposited data
Antibody-Spike complexes	Protein Data Bank (PDB)	see Table 2
SARS-CoV-2 Spike with one chain in up state	(Wrapp et al., 2020)	PDB: 6VSB
SARS-CoV-2 Spike in down state	(Walls et al., 2020)	PDB: 6VXX
SARS-CoV-2 WT RBD complexed with ACE2	(Wang et al., 2020a)	PDB: 6LZG
SARS-CoV-2 N501Y RBD complexed with ACE2	(Zhu et al., 2021)	PDB: 7MJM
SARS-CoV-2 Gamma variant RBD complexed with ACE2	(Dejnirattisai et al., 2021b)	PDB: 7NXC
Furin	(Dahms et al., 2016)	PDB: 5JMO
SARS-CoV-2 WT Spike bound to furin	(Cheng et al., 2021b)	https://zenodo.org/record/4667694#.YGzUgOhKhPY
SARS-CoV-2 Delta Spike bound to furin	This study	https://zenodo.org/record/6026641#.YgQNdN_MI2w
MD trajectories of ACE2 complexed with the WT, Alpha, Beta, Delta or UK2 RBD	This study	https://zenodo.org/record/6026641#.YgQNdN_MI2w
Structual model of REGN10933 bound to the Delta RBD	This study	https://zenodo.org/record/6026641#.YgQNdN_MI2w
Structural model of Nb20 bound to the Delta RBD	This study	https://zenodo.org/record/6026641#.YgQNdN_MI2w
Software and algorithms
ClusPro	(Kozakov et al., 2017)	https://cluspro.bu.edu
SWISS-MODEL	(Waterhouse et al., 2018)	https://swissmodel.expasy.org/interactive
PRODIGY	(Xue et al., 2016)	https://bianca.science.uu.nl//prodigy
PyMOL	(DeLano, 2002)	https://pymol.org/2
HADDOCK 2. 4	(van Zundert et al., 2016)	https://bianca.science.uu.nl/haddock2.4
VMD	(Humphrey et al., 1996)	https://www.ks.uiuc.edu/Research/vmd/vmd-1.9.3
NAMD	(Phillips et al., 2005)	https://www.ks.uiuc.edu/Research/namd
OpenMM	(Eastman et al., 2017)	https://openmm.org
ProDy API 2.0	(Bakan et al., 2011, 2014; Zhang et al., 2021b)	http://prody.csb.pitt.edu
Prism	GraphPad	https://www.graphpad.com/scientific-software/prism
MATLAB	https://www.mathworks.com	Version R2019b