Electrostatic Interactions Explain the Higher Binding Affinity of the CR3022 Antibody for SARS-CoV-2 than the 4A8 Antibody.

A structural understanding of the mechanism by which antibodies bind SARS-CoV-2 at the atomic level is highly desirable as it can tell the development of more effective antibodies to treat Covid-19. Here, we use steered molecular dynamics (SMD) and coarse-grained simulations to estimate the binding affinity of the monoclonal antibodies CR3022 and 4A8 to the SARS-CoV-2 receptor-binding domain (RBD) and SARS-CoV-2 N-terminal domain (NTD). Consistent with experiments, our SMD and coarse-grained simulations both indicate that CR3022 has a higher affinity for SARS-CoV-2 RBD than 4A8 for the NTD, and the coarse-grained simulations indicate the former binds three times stronger to its respective epitope. This finding shows that CR3022 is a candidate for Covid-19 therapy and is likely a better choice than 4A8. Energetic decomposition of the interaction energies between these two complexes reveals that electrostatic interactions explain the difference in the observed binding affinity between the two complexes. This result could lead to a new approach for developing anti-Covid-19 antibodies in which good candidates must contain charged amino acids in the area of contact with the virus.

Introduction

The first outbreak of coronavirus disease 2019 was known in Wuhan, China, in December 2019; then, it became a global pandemic in March 2020 and was named Covid-19.[1] Covid-19 is caused by a novel coronavirus, a severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2).[2] As of 21 March 2021, Covid-19 has resulted in a total of more than 123 million infections and more than 2.7 million deaths (https://coronavirus.jhu.edu/map.html).

Drugs, vaccines, and antibodies can be used to combat Covid-19. However, no new medication has been developed at this time, although several older drugs have been reported to be effective. For example, FDA-approved remdesivir[3] and dexamethasone[4] improve the conditions of severe patients, but they may weaken the immune system.[5] Currently, vaccines developed by various companies such as Pfizer-BioNTech, Moderna, and AstraZeneca are being widely used, but there are cases of resistance and their side effects have not been fully studied. More importantly, Johnson&Johnson (J&J) and Novavax vaccines may not be effective against the South Africa B.1.351 variant of SARS-CoV-2.[6] Antibodies isolated from the plasma of recovered SARS-CoV-2 patients have been proven to effectively treat new patients.[7] However, the amount of plasma available will be insufficient for the growing number of cases, which requires the production of antibodies on a larger scale.

Coronaviruses are spherical in shape with protruding molecules from the viral surface called spike (S) proteins (Figure A,B). The S protein decorates the surface of coronavirus and plays a pivotal role in viral replication by binding to human angiotensin-converting enzyme 2 (ACE2).[8] Antibodies can bind with the S protein, preventing the virus from entering cells (Figure C). The S protein is cleaved into the N-terminal S1 subunit and C-terminal S2 subunit by host proteases and changes conformation from the prefusion to the postfusion state[9],[10] (Figure A). The S1 and S2 subunits comprise an extracellular domain and a single transmembrane helix that function to mediate receptor binding and membrane fusion, respectively.[9,11] Importantly, S1 contains the N-terminal domain (NTD) and the receptor-binding domain (RBD), which are critical in determining tissue tropism and host range.[12],[13] The NTD may recognize specific sugar moieties upon initial attachment[14,15] and plays an important role in the pre- to postfusion transition of the S protein.[16,17] RBD binding to human cells is a critical step, allowing coronaviruses to enter cells to cause infection.[18,19] Since most of the antibodies bind to either NTD or RBD (Figure A), understanding the interactions of antibodies with these regions of SARS-CoV-2 at the atomic level is important for Covid-19 therapies and vaccinations. There are many antibodies that target SARS-CoV-2,[7] but in this study, we focus on two antibodies CR3022 and 4A8 because they hold promise for treating Covid-19 (see below); our computational results may offer insight into the controversial experimental results for these two antibodies, and they bind to different regions of the S protein, allowing more targets to be explored.

Figure 1

(A) Schematic description of the S protein of SARS-CoV-2, which consists of subunits S1 and S2. Monoclonal antibody (mAb) can bind to RBD, NTD, and FP (fusion peptide). (B) S protein of SARS-CoV-2 binds to human ACE2 before its entry to cells. (C) Antibody binds to the S protein, preventing the virus from entering cells.

CR3022, a neutralizing antibody that targets RBD of old SARS-CoV, was previously isolated from a convalescent SARS patient.[20] Recent studies indicated that CR3022 can also bind to RBD of SARS-CoV-2[21] (Figure A), suggesting a potential opportunity to uncover a cross-reactive epitope. Yuan et al. showed that CR3022 can neutralize SARS-CoV but not SARS-CoV-2 RBD at a maximum concentration of 400 μg/mL.[21] Namely, CR3022 binds to SARS-CoV RBD (KD = 1 nM) with a much higher affinity than it does to SARS-CoV-2 RBD (KD = 115 nM) (Table ), implying that CR3022 could not be a candidate for the treatment of SARS-CoV-2.[21] In contrast to Yuan et al., Tian and colleagues found that CR3022 binds efficiently with SARS-CoV-2 RBD (KD = 6.3 nM)[22] (Table ), suggesting that CR3022 alone or in combination with other neutralizing antibodies has potential for prevention and treatment of Covid-19.

Table 1

KD (nM) of CR3022–RBD and 4A8–NTD Complexes Obtained by Experiments and Simulations

complex	experiment results (KD)	our simulation results (KD)
4A8–NTD (PDB ID: 7C2L)	92.7 nM (Chi et al.)[23]	9.1 nM
CR3022–RBD (PDB ID: 6W41)	6.3 nM (Tian et al.)[22]	3.0 nM
	115.0 nM (Yuan et al.)[21]

In contrast to CR3022, 4A8 is a monoclonal antibody that targets SARS-CoV-2 NTD (Figure A) and does not bind RBD.[23] Chi et al. reported that 4A8 is a good candidate for the treatment of Covid-19, as it has a strong neutralizing capacity against both authentic and pseudotyped SARS-CoV-2 NTD (KD = 92.7 nM)[23] (Table ). The epitope of 4A8 on SARS-CoV-2 S protein NTD was determined by cryoelectron microscopy, and its structure in complex with the S protein was obtained with an overall resolution of 3.1 Å and a local resolution of 3.3 Å at the 4A8–SARS-CoV-2 NTD interface.[23] These findings indicate that 4A8 is also a promising therapeutic antibody against SARS-CoV-2 infection.

This work has two goals: (i) to understand the molecular mechanism of CR3022 and 4A8 binding to the S protein and their ability to treat Covid-19 and (ii) to shed light on the controversy between the two experimental groups.[21,22] Using all-atom steered molecular dynamics and coarse-grained simulations, we showed that CR3022 binds strongly to SARS-CoV-2 RBD, especially compared to the binding affinity of 4A8 to SARS-CoV-2 NTD. Our results are consistent with the experimental data of Tian et al.[22] and Chi et al.[23] but not with the data of Yuan et al.[21] (Table ). The difference between the KD values obtained by Tian et al.[22] and Yuan et al.[21] for CR3022 is likely due to the different experimental conditions they used (see below for more details). Therefore, our models perform better under the in vitro conditions adopted by Tian et al.[22]

One of our most important results is that electrostatic interactions are more dominant than van der Waals interactions in antibody binding to the S protein. This may lead to a new strategy for antibody design, according to which the binding region of potential therapeutic agents should be rich in charged residues.

Materials and Methods

Preparation of Input Protein Structures

The structures of CR3022–SARS-CoV-2 RBD (CR3022–RBD) and 4A8–SARS-CoV-2 NTD (4A8–NTD) were extracted from the Protein Data Bank with PDB ID 6W41(21) (for CR3022–RBD) and 7L2C[23] (for 4A8–NTD) (Figure A1, A2). The 4A8–S protein complex (4A8–S protein) was taken from PDB (ID: 7C2L, Figure B2), while the CR3022–S protein complex (CR3022–S protein) was constructed by the ClusPro server,[24,25] where CR3022 was docked to SARS-CoV-2 RBD of the S protein in the open state (was taken from 7L2C without 4A8). The docking result of CR3022 binding to SARS-CoV-2 RBD was then made by a structural alignment with CR3022–RBD (ID: 6W41), and the root-mean-square deviation (RMSD) was 0.1 nm. The missing residues were added by the Modeler package.[26] The structure of CR3022–S protein is shown in Figure B1.

Figure 2

(A1) Structure of the CR3022–RBD complex, retrieved from PDB with ID 6W41. RBD is orange, while green and lemon describe CR3022. (A2) 4A8–NTD complex was extracted from the PDB structure with ID 7L2C. NTD is magenta, and blue and purple-blue refer to 4A8. (B1) Structure of the CR3022–S protein complex, which was obtained by docking CR3022 to the PDB structure 7L2C without 4A8. (B2) PDB structure of the 4A8–S protein complex (PDB ID: 7C2L). The pulling direction in SMD simulations is shown with a spring. The plot was made by the PyMOL 2.0 package.[43]

All-Atom Molecular Dynamics Simulations

The simulation process of the complexes was performed by CHARMM36[27] and AMBER99SB-DISP[28] force fields implemented in the GROMACS 2016 package[29] at 310 K and an isotropic pressure of 1 bar, which was obtained using the v-rescale and Parrinello–Rahman algorithms.[30,31] The Tip3p water model[32] was used in all simulation systems. Bond lengths were constrained by the linear constraint solver (LINCS) algorithm,[33] which allows us to use a time step of 2 fs. The electrostatic and van der Waals interactions were used to depict nonbonded interactions, with the nonbonded interaction pair list being updated every 10 fs using a cutoff of 1.4 nm. The particle mesh Ewald algorithm[34] was used to treat the long-range electrostatic interactions. Periodic boundary conditions were applied in all directions. From these structures, the energy of the system was minimized by the steepest-descent algorithm; then, a short 2 ns MD simulation was performed in the NVT ensemble, which was followed by 3 ns of NPT simulation. Next, a 100 ns production MD simulation was run with an integration time step of 2 fs and the leap-frog algorithm.[35] The “gmx_mpi cluster” tool in the GROMACS 2016 package was used to collect a set of five trajectories for each system to perform steered molecular dynamics simulations.[36−39]

Steered Molecular Dynamics

To investigate CR3022 and 4A8 binding to the S protein, we used SMD, which is as useful as other computationally demanding MD methods in accessing relative binding affinities of ligands.[40−42] This method was also helpful to analyze the interaction between SARS-CoV RBD and human ACE2.[39]

We carried out SMD simulations to pull CR3022 and 4A8 from their RBD and NTD binding regions using the CHARMM36 force field. For each complex, five different trajectories were run at pulling speeds of v = 0.5, 1.5, and 5 nm/ns. To check the robustness of the results to a change in the force field, additional simulations were conducted using the Amber99sb-disp force field.

Rectangular boxes with dimensions of 10 × 6 × 23 nm3 and 7 × 7 × 25 nm3 were used for CR3022–RBD and 4A8–NTD, respectively. However, for the much larger CR3022–S protein and 4A8–S protein complexes, boxes of 18.4 × 21.4 × 37 and 26 × 27.4 × 37 nm3 dimensions were used to allow enough room to pull CR3022 and 4A8 from the binding region. All complexes were immersed in a 0.15 M sodium chloride salt solution to neutralize the total charge.

A spring is attached to a dummy atom on one side and on the other side to the center of mass (CoM) of the antibody (Figure ). The dummy atom is then pulled from its initial position along the line connecting the antibody CoM and CoM of RBD or NTD at a constant speed v. The complexes were rotated so that the unbinding pathway CR3022 and 4A8 is along the z-axis (Figure ), which is displayed using the PyMOL 2.0 package.[43] The pulling force is calculated according to the following equationwhere k is the stiffness of the spring connecting the dummy atom and the antibody CoM, n⃗ is the normal direction of pulling, and r⃗ and are the positions of the system at time t and initial time, respectively. Spring constant k was set to 600 kJ/(mol nm2) (≈ 1020 pN/nm), which is a typical value used in atomic force microscopy (AFM) experiments.[44]

Using the force–displacement profile obtained from the SMD simulation, the pulling work (W) performed by the antibody was estimated using the trapezoidal rulewhere N is the number of simulation steps and F and x are the forces experienced by the target and position at step i.

To estimate the nonequilibrium binding free energy (ΔG), we used Jarzynski’s equality[45] extended to the case when the external force grows at a constant speed v(46)where ⟨...⟩ is the average over N trajectories, z is the time-dependent displacement, and W is the nonequilibrium work at time t, i.e., W = W(t), where W is defined by eq .

In general, using eq , we can extract the equilibrium free energy when the number of simulations is large enough. However, in this study, when the pulling is not slow enough and the number of SMD runs is limited, we can only evaluate the nonequilibrium binding and unbinding energy barriers separating the transition state (TS) from the bound state at t0 and the unbound state at tend, respectively.[40]

Definition of Hydrogen Bond and Nonbonded Contact

A hydrogen bond (HB) is formed when the distance between donor D and acceptor A is less than 0.35 nm, the H-A distance is less than 0.27 nm, and the D-H-A angle is larger than 135°. A nonbonded contact (NBC) between two residues of an antibody and a protein is formed when the shortest distance between their atoms is within 0.39 nm. The two-dimensional (2D) contact networks of HBs and NBCs of CR3022–RBD and 4A8–NTD were constructed using the LIGPLOT package.[47]

Coarse-Grained Simulations

Coarse-Grained Model

The potential energy of the system in this model is given by the following expression:[48]where the terms in order represent the potential energy contributions from bonds, dihedrals, bond angles, electrostatics, native contacts, and non-native interactions, respectively.[49−51] The functional forms of the first three terms have been described in detail previously.[49−51]

The Debye–Huckel theory was employed to model the electrostatic interactions with Debye screening length lD = 1 nm and a dielectric constant of 78.5.[52] Lysine and arginine residues were assigned a charge of +1e, glutamate and aspartate were assigned a charge of −1e, and all other residues were assigned a charge of zero. The contribution from attractive native interactions was computed using the 12-10-6 potential of Karanicolas and Brooks.[53] Collision diameters σ between the Cα interaction sites were set equal to their distance in the crystal structure divided by 21/6. The value of εNC, which sets the depth of the energy minimum for a native contact, was calculated to be εNC = n εHB + ηε. Here, εHB, and ε represent energy contributions arising from hydrogen bonding and van der Waals between residues i and j from the all-atom structure of the protein, respectively. The number of hydrogen bonds n formed between residues i and j is defined using STRIDE software,[54] and εHB is set equal to 0.75 kcal/mol. Intraprotein and interprotein Lennard–Jones (LJ) contacts are defined using a cutoff distance of 0.45 nm between any two heavy atoms of a pair of residues. The value of ε is based on the Betancourt–Thirumalai pairwise statistical potential.[55] While the other terms are transferable among proteins, the LJ well depths for native contacts ε are scaled by a factor η to reproduce the stability of the modeled structures.[48] η values for intraprotein and interprotein interactions of native contacts will be discussed below. All non-native interactions are treated by the final term in the summation using ε = 0.000132 kcal/mol and σ values used as previously described.[56]

Parameterizing the LJ Well Depths for Protein Stability

We applied a previously published training set and parameter tuning procedure to select realistic intra- and interprotein energy scales for native contacts.[56] Sets of ten 1 μs simulations were run with values of η = 1.442, 1.759, 2.480 and 1.235, 1.507, 2.124 for domain and interface of antibodies, respectively, while 1.114, 1.359, 1.916 for SARS-CoV-2 RBD domain and 1.442, 1.759, 2.480 for SARS-CoV-2 NTD domain. The smallest η values were chosen that results in a model that is folded ≥98% of the time in each simulation, and a given conformation was considered to be folded when its fraction of native contacts is greater than 0.69. To assess how strong 4A8 binds to the SARS-CoV-2 NTD domain as compared to CR3022 to the SARS-CoV-2 RBD domain, we plan to employ the same set of parameters for two complexes. For this reason, we set η = 1.442 in simulations instead of 1.359 for the RBD domain, as listed in Table S2. The interaction energy scale for contacts between antibodies and SARS-CoV-2 domains was assigned by an η value of 1.4 to reproduce experimental KD values on the order of nm.

Replica Exchange Umbrella Sampling (REX-US) Simulations

Coarse-grained protein simulations were carried out using Chemistry at Harvard Macromolecular Mechanics (CHARMM) Software, version c35b5.[57] The distance between the centers of mass of the interface residues of the antibody and virus domains is defined as the reaction coordinate. The initial structure of the complex is aligned along the z-axis of the local coordinate system, and the virus domain is translated by 0.05 nm increments along the z dimension to generate a total of 200 umbrella windows. For both complexes, the largest CoM distance for sampling is around 11 nm. A harmonic restraint with a force constant of 70 kcal/mol·Å2 was applied to restrain the relative distance between antibody and virus domain to the target umbrella distance. For each umbrella window, Langevin dynamics simulations were then run at 310 K using a frictional coefficient of 0.050 ps–1, an integration time step of 0.015 ps, and the SHAKE algorithm applied to virtual bonds between coarse-grained particles. Exchanges between neighboring windows were attempted every 5000 integration time steps (75 ps). In total, 10 000 exchanges (750 ns of simulation time) were run with the acceptance ratios between neighboring umbrellas between 0.4 and 0.75. The first 1000 attempted exchanges were discarded to allow for equilibration, and the remaining 9000 exchanges were used for analysis.

Determining Dissociation Constant KD from REX-US Simulations

We can consider CR3022–RBD and 4A8–NTD complexes as two-body systems and define [A], [B], and [AB] as the respective concentrations of the free monomers and the dimer. For example, with CR3022–RBD, [A] = [CR3022], [B] = [SARS-CoV-2 RBD], and [AB] = [CR3022–RBD]. The simulation results are interpreted under the assumption of a two-state binding model. Pb is the probability of the system being in the bound states, with Pu = 1 – Pb defined as the probability of being in the unbound state. In the unbound state, [A] ≡ [B]. The dissociation constant can be calculated as a function of Pb, Pu, and [A] as[39]where the free monomer concentration is defined asIn eq , C0 = 1660 is the standard concentration used to normalize [A] to the units of molarity. V(r*) is the simulation volume in which we found free monomers in the unbound state. As simulations have radial symmetry, r* is the maximum distance between unbound monomers found during the simulation. Pb is calculated from numerical integration of the potential of mean force (PMF) asHere, G1D(r) is the one-dimensional (1D) PMF constructed from the REX simulations using the weighted histogram analysis method (WHAM) equations,[58]rb is the distance threshold separating bound and unbound states, β = 1/kBT, and rb is the value of r at which the 1D PMF reaches maximum.

Results and Discussion

CR3022–RBD Is More Stable than 4A8–NTD: Analysis Based on PDB Structures

Using PDB structures with missing residues rebuilt as described in the Section , we obtained the CR3022–RBD and 4A8–NTD interfaces shown in Figure S1A1,A2. RBD (chain C) formed contact with both chains A and B of CR3022, while NTD (chain A) interacted with chain B but not with chain C of 4A8. There were more than 42 residues in the CR3022–RBD binding region, while only 27 residues were present at the 4A8–NTD interface.

The network of hydrogen bonds and nonbonded contacts of CR3022–RBD is richer than 4A8–NTD (Figure S1B1,B2). There are 10 and 6 HBs for CR3022–RBD and 4A8–NTD, respectively, while CR3022–RBD has 21 nonbonded contacts, which is more than 14 contacts between 4A8 and NTD. Thus, CR3022 is likely to associate with RBD more strongly than 4A8 with NTD, and this will be confirmed by molecular simulations.

SMD Results

CR3022 Binds to RBD More Strongly than 4A8 to NTD: CHARMM36 FF

We first discuss the results obtained by using the CHARMM36 force field. The force–time profiles of two complexes, obtained at v = 0.5 nm/ns (Figure ) and 1.5 and 5 nm/ns (Figure S2), show that CR3022 binds to RBD more strongly than 4A8 to NTD as the corresponding rupture force Fmax is higher. For v = 0.5 nm/ns, Fmax = 1665.2 ± 121.3 and 638.2 ± 57.1 pN for CR3022–RBD and 4A8–NTD, respectively. As expected, the rupture force increases with increasing pulling speed (Table ).

Figure 3

(Left) Pulling force, pulling work, and nonequilibrium energy profiles of CR3022–RBD and 4A8–NTD complexes averaged from five independent SMD runs at v = 0.5 nm/ns. (Right) Pulling force, pulling work, and nonequilibrium energy profiles of CR3022–S protein and 4A8–S protein complexes averaged from five independent SMD runs at v = 0.5 nm/ns. The results were obtained by the CHARMM36 force field.

Table 2

Rupture Force (Fmax), Unbinding Time (tmax), Work of the External Force (W), Nonequilibrium Binding (ΔGbind), and Unbinding (ΔGunbind) Free Energy Barriers Obtained from the Five Independent SMD Trajectories of Four Complexes at Pulling Speeds v = 0.5, 1.5, and 5 nm/nsa

	Fmax (pN)		tmax (ps)		W (kcal/mol)		ΔGbind(kcal/mol)		ΔGunbind (kcal/mol)
pulling speed v (nm/ns)	CR3022–RBD	4A8–NTD	CR3022–RBD	4A8–NTD	CR3022–RBD	4A8–NTD	CR3022–RBD	4A8–NTD	CR3022–RBD	4A8–NTD
0.5	1665.2 ± 121.3	638.2 ± 57.1	6094.0 ± 218.3	2922.0 ± 170.9	438.2 ± 5.9	152.2 ± 4.0	313.5 ± 4.0	79.6 ± 4.1	313.1 ± 4.5	78.4 ± 1.9
1.5	1843.5 ± 146.4	1001.4 ± 85.3	2015.8 ± 131.2	1265.8 ± 154.8	607.2 ± 5.6	305.9 ± 5.7	393.8 ± 5.1	131.1 ± 5.5	387.7 ± 2.6	129.0 ± 5.9
5	2437.5 ± 155.1	1354.8 ± 108.7	726.3 ± 54.8	463.2 ± 57.5	1076.4 ± 7.7	647.7 ± 6.9	498.5 ± 5.6	216.9 ± 7.8	478.9 ± 7.9	189.9 ± 6.7

The errors represent standard deviations. The results were obtained by the CHARMM36 force field.

Because CR3022 binds to RBD more strongly than 4A8 to NTD, the rupture time tmax to reach Fmax of the CR3022–RBD complex is also larger than that of 4A8–NTD (Table ), which is consistent with the results obtained for protein–ligand systems.[39] The rupture time decreases with an increase of v.

The nonequilibrium work W was shown to be a better value for characterizing the relative binding affinity than Fmax.[59] It rapidly increased until CR3022 and 4A8 come out from the binding region and reached a stable value when the antibody ceased to interact with the spike protein (Figures B and S2). At v = 0.5 nm/ns, we obtained W = 438.2 ± 5.9 and 152.2 ± 4.0 kcal/mol for CR3022–RBD and 4A8–NTD, respectively (Table ), which implies that, in agreement with the experiments of Tian et al.[22] and Chi et al.,[23] the former complex is more stable than the latter. This also agrees with the data obtained for higher pulling speeds v = 1.5 and 5 nm/ns (Table ).

Using eq , we obtained the time dependence of the nonequilibrium binding free energy ΔG for two complexes at v = 0.5 nm/ns (Figure C). The maximum corresponds to the transition state with ΔGTS. Since at the beginning the state was bound, we have ΔGbound = ΔG(t0) ≈ 0 kcal/mol, while the unbound state occurs at the end of simulation[40] and ΔGunbound = ΔG(tend) ≈ 0 kcal/mol. The binding and unbinding free energy barriers, which are defined as ΔΔGbind = ΔGTS – ΔGunbound and ΔΔGunbind = ΔGTS – ΔGbound, are nearly equal as ΔGunbound ≈ ΔGunbound ≈ 0.

From Figure C, we obtained ΔΔGunbind = 313.1 ± 4.5 and 78.4 ± 1.9 kcal/mol for CR3022–RBD and 4A8–NTD, respectively (see also Table ), providing additional evidence that CR3022 binds to RBD more tightly than 4A8 to NTD. This conclusion is also valid for other pulling speeds (Table and Figure S2).

CR3022 Binds to RBD More Strongly than 4A8 Binds to NTD: AMBER99SB-DISP FF

To test the robustness of our results against force fields, we additionally performed simulations with the AMBER99SB-DISP force field. The results are shown in Figures S4–S6 and Table S1. Fmax, W, and ΔΔGunbind obtained for three pulling speeds also support the view that the affinity of binding of CR3022 to RBD is higher than that of 4A8 to NTD.

CR3022 Binds to RBD More Strongly than 4A8 to NTD: Effect of the Entire S Protein Structure

So far, we have considered the interaction of an antibody with either RBD or NTD neglecting the rest of the entire S protein. In this section, we ask whether the remainder of the S protein influences our main conclusion that CR3022 binds to the target more strongly than 4A8. To answer this question, we performed SMD simulations for the CR3022–S protein and 4A8–S protein complexes, as shown in Figure B1,B2.

As can be seen from Figures and S6, CR3022 binds to the S protein tighter than 4A8, having higher values of the rupture force, rupture time, pulling work, and unbinding free energy (see also Table ). However, the interaction of the antibody with the entire S protein is weaker compared to NTD and RBD. For example, at v = 0.5 nm/ns, Fmax = 490.8 ± 36.3 and 241.4 ± 21.1 pN for the CR3022–S protein and 4A8–S protein, respectively, while Fmax = 1665.2 ± 121.3 and 638.2 ± 57.1 pN for CR3022–RBD and 4A8–NTD, respectively (Table ). The same trend was obtained for W and ΔΔGunbind at all pulling speeds, while tmax becomes larger due to weaker interactions (Table ).

Thus, taking into account the entire structure of protein S changes the absolute binding affinity but leaves the relative binding affinity unchanged. This suggests that CR3022 is a better candidate for the treatment of Covid-19 than 4A8.

Binding of CR3022 and 4A8 to the S Protein Is Driven by Electrostatic Interactions

The time dependence of vdW, electrostatic, and total (vdW + electrostatic) interaction energies of the CR3022–RBD, 4A8–NTD, CR3022–S protein, and 4A8–S protein complexes, obtained at v = 0.5, 1.5, and 5 nm/ns, is shown in Figures and S7. There is a small difference in vdW interaction energies of CR3022–RBD and 4A8–NTD, but a much more pronounced difference is observed for the electrostatic interactions. The same is true for CR3022–S protein and 4A8–S protein complexes. It is important to note that for all complexes, the energy of electrostatic interactions (Eelec) is significantly lower than the vdW energy (EvdW), which means that their stability is primarily determined by electrostatic interactions.

Figure 4

(A1) Total interaction energy (sum of electrostatic and vdW) for the CR3022–RBD and 4A8–NTD complexes. (A2) Same as in (A1) but for the larger CR3022–S protein and 4A8–S protein complexes. (B1) Electrostatic (Eelec) and vdW (EvdW) interaction energies for the CR3022–RBD and 4A8–NTD complexes. (B2) Same as in (B1) but for larger CR3022–S protein and 4A8–S protein complexes. The results were obtained from five independent SMD runs at v = 0.5 nm/ns using the CHARMM36 force field.

We calculated the mean interaction energy in the bound state by averaging over the time window [0, tmax]. Note that tmax depends on the system, v, force field, and SMD runs (Tables and S1). At v = 0.5 nm/ns for CR3022-RDB, we obtained Eelec = −299.6 ± 1.4 kcal/mol, which is clearly lower than EvdW = −100.6 ± 0.9 kcal/mol (Table ). A similar result was obtained for CR3022 interacting with the entire S protein with Eelec = −299.7 ± 1.9 kcal/mol, which is clearly lower than EvdW = −96.2 ± 1.4 kcal/mol (Table ). The difference between the electrostatic and vdW interactions is much more pronounced in the 4A8 case, namely, Eelec = −817.1 ± 2.4 kcal/mol and EvdW = −58.3 ± 0.6 kcal/mol for 4A8–NTD, whereas Eelec = −1001.9 ± 2.3 kcal/mol and EvdW = −61.4 ± 0.3 kcal/mol for 4A8–S protein (Table ). The dominant role of the electrostatic interactions remains true for other pulling speeds.

Table 3

Interaction Energies Obtained by Averaging Five Trajectories of Four Complexes in the Time Window [0 – tmax] at Pulling Speeds v = 0.5, 1.5, and 5 nm/nsa

	v = 0.5 nm/s		v = 1.5 nm/s		v = 5 nm/s
interaction energy (kcal/mol)	CR3022–RBD	4A8–NTD	CR3022–RBD	4A8–NTD	CR3022–RBD	4A8–NTD
electrostatic	–299.6 ± 1.4	–817.1 ± 2.4	–296.3 ± 1.2	–820.9 ± 2.3	–287.1 ± 2.5	–830.3 ± 2.2
vdW	–100.6 ± 0.9	–58.3 ± 0.6	–98.4 ± 0.5	–58.9 ± 0.9	–95.6 ± 1.2	–58.6 ± 0.9
total	–400.2 ± 2.3	–875.4 ± 3.0	–394.7 ± 1.7	–879.8 ± 3.2	–382.7 ± 3.7	–888.9 ± 3.1

The errors represent standard deviations.

The importance of electrostatic interactions has also been recognized for the entry of SARS-CoV-2 into cells through binding of the S protein to human ACE2.[39,60]

Charged Residues at the Interface Are Important for Binding Affinity

To understand the role of each residue at the interface in binding of the antibody to the S protein, we calculated the average per-residue interaction energy in the [0, tmax] time window at a pulling speed v = 0.5 nm/ns. For both 4A8–NTD and 4A8–S protein complexes, Lys147(A), Lys150(A), and Arg246(A) of the S protein and Glu31(B), Glu54(B), Asp55(B), and Glu72(B) of the 4A8 antibody yielded a significant contribution to the interaction energy (Figure , top). In the case of CR3022–RBD and CR3022–S protein complexes (Figure , bottom), residues Asp55(A), Glu57(A), Asp107(A), and Glu61(B) of the CR3022 antibody and Lys378(C) and Lys386(C) of the S protein dominate. Importantly, all of the most prominent residues are charged, implying that electrostatic interactions play a dominant role in stabilizing the four complexes studied.

Figure 5

(Top) Total interaction energy of the residues at the binding region (see Figure S2) of the 4A8–NTD and 4A8–S protein complexes. (Bottom) Same as on the top but for the CR3022–RBD and CR3022–S protein complexes. Black and red refer to the residues of SARS-CoV-2 and antibody, respectively. The results were obtained in the time window [0, tmax] at pulling speed v = 0.5 nm/ns. The CHARMM36 force field was used.

To show that the most important charged residues govern the binding affinity, we carried out simulations where these residues have been replaced by neutral Alanine (Figure S8). For v = 0.5 nm/ns and the CHARMM36 force field, these mutations reduced Fmax from 1665.2 ± 121.3 to 768.9 ± 62.6 pN for CR3022–RBD and from 638.2 ± 57.1 to 403.1 ± 51.5 pN for 4A8–NTD, which supports our hypothesis.

Estimation of Dissociation Constant of CR3022–RBD and 4A8–NTD Complexes Using Coarse-Grained Simulations

As described in the Section , we employed dissociation constant KD to evaluate the binding affinity of antibodies to SARS-CoV-2 RBD and SARS-CoV-2 NTD domains. First Pb is calculated from 1D PMF (eq ) at different r* values. Figure A shows the most stable state locates near the native states with the CoM distance ≈ 1 nm (CR3022–RBD) and ≈0.86 nm (4A8–NTD). The barrier of 1D PMF separating the bound and unbound regimes occurs at ≈3 nm for both complexes so we decided to choose rb = 3 nm for the numerical calculation of numerator of eq . Pu = 1 – Pb, and free monomer concentration [A] is calculated using eq . Finally, KD is obtained from eq . Figure B plots KD curves as a function of r*. As expected, KD increases and converges at a large radius r*, which means KD physically should not depend on r*. From our simulations, we need to define a cutoff r* corresponding to a total limit volume to define the probability of finding the system in the free monomer state. We determined r* ≈ 11 nm from which there is no longer interaction between the antibody and virus and used this value to calculate KD.[39] Our results showed that both antibodies tightly bind to virus domains with KD = 9.1 nM for 4A8–NTD and KD = 3 nM for CR3022–RBD (Table ). Thus, from our simulations it can be seen that CR3022 binds to the SARS-CoV-2 RBD domain more strongly than 4A8 binds to SARS-CoV-2 NTD, but the difference is not as great as in SMD (Table ), because the difference in KD is only about three times. However, it can somehow explain the discrepancy between the results obtained in different experiments. Tian’s group[22] measured the strong binding affinity of the CR3022–RBD complex with a dissociation constant KD = 6.3 nM. Comparing this result with KD = 92.7 nM obtained for 4A8–NTD by Chi et al.,[23] one can conclude that CR3022 binds strongly to SARS-CoV-2 RBD than 4A8 to SARS-CoV-2 NTD. Thus, our results obtained by both all-atom and CG simulations are consistent with these two groups.

Figure 6

(Left) One-dimensional potential of mean force (1D PMF) of CR3022–RBD (black curve) and 4A8–NTD (red curve). Results were obtained by applying the WHAM analysis for 750 ns REX-US simulations. (Right) KD curves as a function of r* corresponding to the change in the total free monomer concentration from eq . Pb and KD were determined at r* = 11 nm.

It should be noted that for 4A8–NTD our value of KD is about an order of magnitude smaller than that of Chi et al.[23] (Table ). This level of agreement between simulations and experiments is reasonable when we consider the relationship between KD and the binding free energy ΔGbind = −kBT  ln(KD), where KD is measured in M. At room temperature, kBT ≈ 0.592 kcal/mol, meaning that a difference in KD of one order of magnitude results only in a difference in ΔGbind of 1.4 kcal/mol, which is on the order of the calculation error.

Experimentally, Yuan’s group[21] pointed out that CR3022 has a weak neutralizing effect for SARS-COV-2, as it has a low binding affinity to RBD with KD = 115 nM (Table ). This result is in conflict with that of Tian et al.,[22] who reported a lower value of KD.

One advantage of our computational study is that we studied two complexes using the same model, while experiments conducted by different groups were carried out under different conditions, making it difficult to directly compare experimental results. Our simulations showed that CR3022 binds to RBD more strongly than 4A8 to NTD, which is consistent with Tian et al.[22] and Chi et al.[23] but not with Yuan et al.[21] (Table ). From our computational point of view, the fact that the result of Yuan et al.[21] contradicts that of Tian et al.[22] is explained by the different conditions used in their in vitro experiments, namely, CR3022 was expressed in mammalian cells in Yuan et al.[21] but in Escherichia coli in Tian et al.[22] Moreover, SARS-CoV-2 RBD was obtained from insect cells in Yuan et al.[21] but from mammalian cells in Tian et al.[22] Since KD obtained in our CG simulations is close to that of Tian et al.,[22] our model can be expected to reasonably capture the conditions used in this group’s experiment. A modification of our models to mimic the experimental conditions in Yuan et al.[21] is nontrivial and requires further study.

Conclusions

In this work, we applied all-atom SMD and coarse-grained simulations to study the binding affinity of CR3022 and 4A8 antibodies to the S protein of SARS-CoV-2. SMD simulations showed that CR3022 displays a higher binding propensity to RBD than 4A8 to NTD, which is consistent with the result obtained by coarse-grained REX-US simulations that the dissociation constant KD of CR3022–RBD is approximately three times smaller than that of 4A8–NTD. Our results are in good agreement with the experimental data of Tian et al.[22] and Chi et al.,[23] but they are in contrast to the experimental results of Yuan et al.[21] The contribution of electrostatic interactions to the stability of four complexes, including CR3022–RBD, 4A8–NTD, CR3022–S protein, and 4A8–S protein, is more significant compared to vdW interactions. In terms of binding capability, CR3022 is a better candidate for Covid-19 treatment than 4A8.

Since the RBD and NTD binding sites contain charged residues, electrostatic interactions are likely to play an important role not only in CR3022 and 4A8 binding but also in other antibodies. Our preliminary results on the binding of antibodies REGN10933 and REGN10987,[61] as well as nanobodies H11-H4,[62] support this hypothesis, but more systems need to be examined to arrive at a firm conclusion.

Our prediction that electrostatic interactions play a key role in the binding of antibodies to SARS-CoV-2 could open up a new strategy for developing effective antibodies against Covid-19. For example, good candidates should contain many charged amino acids in the region that binds to the spike protein. Moreover, since the important residues of the spike protein are positively charged (Lys and Arg, Figure ), potential antibodies must have negatively charged residues (Asp and Glu). From a methodological point of view, it is important to emphasize that coarse-grained models in combination with REX-US provide a reasonable tool for estimating the dissociation constant of two proteins.