Critical Differences between the Binding Features of the Spike Proteins of SARS-CoV-2 and SARS-CoV.

The COVID-19 caused by SARS-CoV-2 has spread globally and caused tremendous loss of lives and properties, and it is of utmost urgency to understand its propagation process and to find ways to slow down the epidemic. In this work, we used a coarse-grained model to calculate the binding free energy of SARS-CoV-2 or SARS-CoV to their human receptor ACE2. The investigation of the free energy contribution of the interacting residues indicates that the residues located outside the receptor binding domain are the source of the stronger binding of the novel virus. Thus, the current results suggest that the essential evolution of SARS-CoV-2 happens remotely from the binding domain at the spike protein trimeric body. Such evolution may facilitate the conformational change and the infection process that occurs after the virus is bound to ACE2. By studying the binding pattern between SARS-CoV antibody m396 and SARS-CoV-2, it is found that the remote energetic contribution is missing, which might explain the absence of cross-reactivity of such antibodies.

Introduction

The 2019 coronavirus disease (COVID-19) caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has been spread globally since its first outbreak in Wuhan, China in December 2019.[1−3] Common symptoms of SARS-CoV-2 infected patients involve fever, cough, and fatigue with an estimated death rate of 3%–5%.[4] It has already caused more than 7 million confirmed cases and more than 400,000 deaths all over the world (at the time this work is written). In addition to health concerns, the disease also brought severe economic and social issues.[5] While the situation of Wuhan has been stepped down, Europe and the United States are experiencing a major epidemic. Thus, it is crucial to understand the infection and spreading process and find measures to mitigate the epidemic situation.

The SARS-CoV-2 is a member of the beta-coronavirus genus,[6] which also includes the acute respiratory syndrome coronavirus (SARS-CoV) and the middle-east respiratory syndrome virus (MERS). The SARS-CoV-2 virus appears to be optimized for binding to the human receptor ACE2,[7] and the binding patterns between ACE2 and SARS-CoV-2 or SARS-CoV at the receptor binding domain (RBD) are thought to be almost identical.[8] More specifically, SARS-CoV-2 shares 76%–78% sequence identity with SARS-CoV for the whole protein and 73%–76% for the RBD.[9] The trimeric spike glycoprotein of SARS-CoV-2 is comprised of three S1/S2 units, and the RBD locates at S1. One variation, the S1/S2 cleavage site of SARS-CoV-2 is a unique “RRAR” furin recognition site,[10] while in SARS-CoV it is a single arginine.[11] The three S1/S2 units undergo a hingelike conformational switch between “up” and “down” states. Only at the “up” state, the RBD is exposed and is able to bind to the receptor, while at the “down” state the RBD is hidden and is inaccessible by the receptor.[12] For SARS-CoV, the spike trimer with two “down” and one “up” is the most populated state.[13] This could very likely be the case for SARS-CoV-2, but to our knowledge no experimental statistical measurement has been reported yet. A recent study pointed out the possibility of two spike proteins binding with the same ACE2.[14] After binding to the receptor, the following cascade of events is triggered: the spike protein undergoes a large conformational change, the S1 with the receptor is shed, S2 is transformed to a more stable postfusion state, and finally the viral membrane is fused with the cell membrane.[15,16]

Despite the similarities in structures and binding patterns between the two viruses, SARS-CoV-2 spreads faster than SARS-CoV and this might be due to the stronger binding of the ACE2-SARS-CoV-2 complex.[12] The range of experimental binding affinities of the two ACE2-virus is wide, with reports of 15 nM[12] and 150–185 nM[13] for the ACE2-SARS-CoV-2 and the ACE2-SARS-CoV complexes, respectively, and also reports of 4.7 nM and 31 nM,[8] and 1.2 nM and 5.0 nM[17] for both systems, respectively. In all cases, the ACE2-SARS-CoV-2 complex shows a larger binding affinity than ACE2-SARS-CoV. On the other hand, although the sequences and epitope have been studied extensively, it is still unclear what is the structural/energetic basis for the difference between the two complexes. Moreover, the receptor binding is a crucial step for drug and antibody interference with the infection process. Thus, this work will focus on understanding the detailed differences between the binding features of the two coronaviruses and the human receptor ACE2.

Recent works yielded high-resolution structure of SARS-CoV-2 at its prefusion state,[14] as well as the complex of its RBD domain and ACE2.[8] These emerging structures provide an opportunity to use computational modeling to investigate the underlying mechanism behind the differences in binding strengths of the two ACE2-virus complexes.

However, such a task is very challenging. For example, recent theoretical work[18] analyzed the number of contacts, interface area, and fluctuations and concluded that different viruses have different strategy for binding. However, this work could not obtain the correct order of binding affinity between the ACE2-SARS-CoV-2 and ACE2-SARS-CoV complexes.

Obviously, the main issue is the differences in interaction of free energies between the two types of viruses and the receptor, as this energy is essential for understanding the binding process. Evaluating the binding energy of a very large protein complex is an enormous challenge for fully atomistic models and thus we chose here to use our coarse-grained (CG) model[19−21] to study the energetics of the complexes. Our CG model has been consistently developed and systematically calibrated to account for proper evaluation of electrostatic free energies of proteins and membranes, including, of course, solvation and hydrophobic effects. The model was applied extensively to many systems, calculating protein-folding free energies and related properties.[22−24] Here we use the model to evaluate the binding of the virus to ACE2.

Our analysis of the binding patterns found that the substitutions of residues near the RBD of SARS-CoV in the conversion to SARS-CoV-2 is not the reason for the increase in binding energy. It is found that the major contribution actually comes from the body of the spike protein that is away from the RBD. It is also found that the antibody of SARS-CoV, that did not show cross-reactivity, might be partially due to the fact that the binding interface is partially covered, compare to the situation in the ACE2-SARS-CoV-2 complex.

Methods

In this work, we used Modeler[25] to perform homology modeling in constructing the binding complexes of ACE2-SARS-CoV-2 and m396-SARS-CoV-2. The SARS-CoV-2 structure was taken from the recent cryo-EM study (PDB ID: 6VSB)[12] with an incomplete receptor binding domain. For the binding domain, we used the crystal structure of SARS-CoV-2 RBD that is bound to ACE2 (PDB ID: 6M0J).[8] The binding between m396 and SARS-CoV-2 was modeled using the m396-SARS-CoV structure as a template (PDB ID: 2DD8).[26]

Subsequently, we utilized our CG model to calculate the free energy of each structure and the relevant binding energies. Our CG model is focused on the electrostatic free energy of the protein that involves the solvation energy and the interactions between charged and polar residues. The full CG treatment includes membrane terms (see SI) which are not included in the present case, since the membrane is out of the system studied. The total CG energy is defined as follows[20] (see also SI)

The terms on the right are the side chain van der Waals energy, main chain solvation energy, main chain hydrogen bond energy, side chain electrostatic energy, side chain polar energy, side chain hydrophobic energy, main chain/side chain electrostatic energy, and main chain/side chain van der Waals energy, respectively. The scaling coefficients c1, c2, and c3 are taking values of 0.10, 0.25, and 0.15 in this work.

Before the evaluation of the free energy, we used a Monte Carlo proton transfer (MCPT) method[20] to determine the charge configuration of all ionizable residues in the system. In the MCPT approach, the MC procedure controls proton transfer moves between ionizable residues or between one ionizable residue and the bulk. The acceptance probability of the move is determined by standard Metropolis criteria (see SI). By such calculations, we are able to get the CG free energy of each protein configuration and also the electrostatic contribution of each residue when they are either in the ACE2-virus complex or in the unbound state. Note that the CG already represents the free energy of the system and not the potential energy. All calculations and simulations were carried out using the MOLARIS-XG package.[27,28] For more details, see SI.

Results and Discussion

We started by trying to evaluate the binding energy differences between the two complexes (ACE2-SARS-CoV-2 and ACE2-SARS-CoV) and to determine where the difference is coming from. This study utilized the recently published cryo-EM structure of SARS-CoV-2 (PDB ID: 6VSB)[12] and the crystal structure of its RBD and ACE2 (PDB ID: 6M0J)[8] and performed homology-modeling in order to obtain the structure of the ACE2-SARS-CoV-2 complex. The structure of ACE2-SARS-CoV has been taken from a previous work (PDB ID: 6CS2).[13] After getting the structures, we performed energy minimization and molecular dynamics for structural relaxation. This procedure was followed by the MCPT algorithm[20] (that determines the optimal charge distribution of the ionizable residues) and obtained the CG energies of the ACE2-SARS-CoV-2 and ACE2-SARS-CoV complexes. With the same treatment, we obtained the CG energies of the ACE2, SARS-CoV, and SARS-COV-2 monomers at infinite distance separation. The ACE2-virus binding free energy was then calculated by ΔGbinding = Gcomplex – GACE2 – Gvirus.

Even though the binding mode of the two ACE2-virus complexes (Figure B) were argued to be almost identical,[8] we still obtained binding energy of −70.7 kcal/mol for ACE2-SARS-CoV-2 and −66.4 for ACE2-SARS-CoV, respectively. As expected, the major difference come from electrostatic contributions. At this point, the binding energy difference could either be an effect of the nonconserved residues or the change in structures of the two complexes near or outside the RBDs.

Figure 1

(A) The structure of the ACE2-SARS-CoV-2 complex. ACE2 is in blue and SARS-CoV-2 is in cyan. (B) The overlap of the receptor binding domain of the ACE2-SARS-CoV-2 and ACE2-SARS-CoV complex. The ACE2 region of the SARS-CoV-2 complex is shown in blue and SARS-CoV-2 spike protein is shown in cyan. The ACE2 region of the SARS-CoV complex is shown in red and the SARS-CoV spike protein is shown in orange.

To understand this issue in a more quantitative way, we evaluated the electrostatic contributions of each residue to the total binding free energy of the two complexes (see Figure A). It was found that some residues give positive contributions while others give negative contributions. This finding is consistent with the results that show some interactions at the binding interface strengthen while others weaken the binding.[14] To see whether the difference in binding comes from the RBD, we further classified the free energy contributions based on their distances from the binding site. Thus, we plotted in Figure B the contributions of the residues within given range to the total free energy (this was done according to the distances to the N501 residue of ACE2-SARS-CoV-2 or to T487 of ACE2-SARS-CoV). To our surprise, if we consider the residues within 60 Å of the binding site, it is found that the ACE2-SARS-CoV complex has a stronger binding affinity (its curve is below that of the ACE2-SARS-CoV-2 curve). Apparently the residues between 60 and 120 Å are the residues that switched the trend. This result indicates that the binding to the receptor in a remote position from the binding site possibly leads to a stronger binding in ACE2-SARS-CoV-2 than in ACE2-SARS-CoV. Interestingly, the three S1/S2 cleavage sites also locate within the 60–120 Å range.

Figure 2

(A) Electrostatic free energy changes of each residue before and after formation of the ACE2-SARS-CoV-2 and the ACE2-SARS-CoV complexes. (B) The dependence of the binding energy contributions on their distances to the binding interface (N501 for SARS-CoV-2 and T487 for SARS-CoV). (C) Classification of charged residues based on their energy contributions to the ACE2-SARS-CoV-2 complex or (D) to the ACE2-SARS-CoV complex.

In view of the above conclusion, we tried to examine whether the finding of long-range interactions is coincidental. Thus, we tested the influence of effective dielectric constant in the CG model (using a constant between 60 and 90) and also use a function of the form εeff = 1 + ε′[1 – exp(−μr)], where we used ε′ = 80 and μ = 0.8. It was found that the trend stayed similar to the different dielectric constants. It is still possible that including the effect of the ionic strength would reduce the binding difference at the long distance, but the sign is very unlikely to change.

It should be noted that the value of the overall calculated binding energy is most probably an overestimate. One missing effect is the entropy contribution of the separate parts of the complex (which is equal for the two systems). Another missing effect is the above-mentioned effect of the ionic strength that would reduce the electrostatic interaction.

To illustrate the contributions of each residue, we plot their contributions in Figure C,D. The figure assigns each residue one color by its free energy contribution. Both complexes have residues with relatively large energy changes near the binding site (red and yellow colors) and according to Figure B these contributions are more negative for ACE2-SARS-CoV-2. However, when we move outside the RBD, more residues of this type appear in the thinner spike protein body of SARS-CoV2 compared to the fatter spike protein body of SARS-CoV. It suggests that for both ACE2-virus complexes, some residues at the RBD region strengthen the binding while others weaken it (Figure A). However, the interactions change between 60 to 120 Å at the binding site makes the binding affinity of SARS-CoV-2 stronger (Figure B) and may indicate an effective evolution of the spike protein body of the novel virus. This larger binding affinity might explain why SARS-CoV-2 spreads faster than SARS-CoV.

Because of the structural homology and similar binding patterns, several RBD-directed monoclonal antibodies (m396, S230, 80R) of SARS-CoV have been tested for SARS-CoV-2, but none of them could show cross-reactivity.[8,12] Another antibody (CR3022) that was obtained from a convalescent SARS-CoV patient could bind to SARS-CoV-2 but still could not neutralize the virus even at a concentration as high as 400 μg/mL, and its cross-reactivity was attributed to the high percentage of targeted epitope residues (86%).[29]

In this work we tried to understand the absence of cross-reactivity by studying the binding pattern between SARS-CoV-2 with one of the SARS-CoV antibody, m396. We used homology modeling to generate a structure of the m396-SARS-CoV-2 complex from SARS-CoV-2 (PDB ID: 6VSB)[12] and m396-SARS-CoV (RBD) (PDB ID: 2DD8)[26] structures. For the antibody-virus complexes, we used ΔGbinding = Gcomplex – Gm396 – Gvirus. Figure A shows the overlapped structures of the ACE2-SARS-CoV-2 and m396-SARS-CoV-2 complexes that are aligned toward the virus body. Visually the m396 antibody only covers part of the ACE2-SARS-CoV-2 binding interface. As before, we analyzed the distance-dependent electrostatic energy contributions of the m396-virus complex. As seen from Figure B, the binding of m396 and SARS-CoV-2 does not result in structural/energetic differences that can lead to an increase of interactions in the range between 60 and 120 Å, as was observed in the ACE2-SARS-CoV-2 complex. The energy contribution near the RBD (<60 Å) is also weaker in comparison to the corresponding contribution in the ACE2-virus complex. Overall, m396 shows an ineffective binding pattern that misses part of the epitopes of ACE2 that could trigger the following structural changes, and this might also be the case for other antibodies that did not show cross-reactivity. To mimic the ACE2-SARS-CoV-2 binding pattern, Zhang et al. synthesized a 23-mer peptide fragment of the ACE2 peptidase domain α1 helix.[30] However, the binding affinity was not strong. This might be an effect of the instability of a small helix fragment.

Figure 3

(A) The overlap of m396-SARS-CoV-2 and m396-SARS-CoV complexes. The m396 of the SARS-CoV-2 complex is shown in blue and the SARS-CoV-2 spike protein is shown in cyan. The m396 of the SARS-CoV complex is shown in red and the SARS-CoV spike protein is shown in orange. (B) The dependence of the binding energy contributions on their distances to the binding domain.

Structural analysis shows that there are 14 key residues that participate in the binding between ACE2 and SARS-CoV.[17] Among them, eight are conserved in SARS-CoV-2 and the other six are mutated. The six mutated residues are N439/R426 (SARS-CoV-2/SARS-CoV), L455/Y442, F486/L472, Q493/N478, Q498/Y484, and N501/T487 (Figure A). To understand how the substituted residues would affect binding energy, we performed mutation calculations for the two ACE2-virus complexes. The residues of ACE2-SARS-CoV-2 were mutated to the one in the corresponding position in ACE2-SARS-CoV and vice versa. After introducing the mutation, we performed another relaxation run before the CG MCPT free energy evaluation. As shown in Table , all mutated ACE2-SARS-CoV-2 constructs are less stable compared to the wild-type system, while the mutated ACE2-SARS-CoV shows the opposite results. This is consistent with our previous distance-dependent binding energy contributions near the RBD, where we find that the ACE2 and SARS-CoV complex gives a more favorable binding pattern than in the ACE2-SARS-CoV-2 complex (Figure B). The mutation calculation suggests that the effective binding of the novel virus might happen at the remote part.

Figure 4

Positions of key residues at the ACE2-SARS-CoV-2 binding domain that are substituted compared to the SARS-CoV system.

Table 1

The CG Results for Mutating Residues from Their Sequence in ACE2-SARS-CoV-2 (N439, L455, F486, Q493, Q498, N501) to the Sequence in ACE2-SARS-CoV (R426, Y442, L472, N479, Y484, T487) and Vice Versaa

residue #	ACE2-SARS-CoV-2 mutation	ACE2-SARS-CoV mutation
wild-type	–70.71	–66.45
N439/R426	–70.79	–61.98
L455/Y442	–75.61	–62.20
F486/L472	–72.39	–53.23
Q493/N479	–75.14	–63.79
Q498/Y484	–73.90	–57.10
N501/T487	–73.80	–58.08
All 6	–73.05	–58.81

Energy unit in kcal/mol.

The current work explored the structural/energetic basis of the difference in binding energy between the SARS-CoV-2/ACE2 and the SARS-CoV/ACE2 complexes. It was found that the binding of SARS-CoV-2 is more favorable, not because its RBD has been optimized, but because the SARS-CoV-2 glycoprotein trimer body has been evolved to bind stronger at distanced sites (in fact, if we just consider the RBD then SARS-CoV is more favorable). It is not very clear if this stronger binding is converted to conformational changes, cleavage, and subsequent fusion events. However, if this is the case we have an interesting explanation of the reason why the novel coronavirus spreads faster and easier. The results also suggest to use the novel virus as a template during drug/antibody design with the whole spike protein as a binding template instead of a fraction of the RBD that might neglect the essential changes in the virus body and size effect.