Electrostatic Interactions Are the Primary Determinant of the Binding Affinity of SARS-CoV-2 Spike RBD to ACE2: A Computational Case Study of Omicron Variants

Round 1

Reviewer 1 Report

The study by Sang et al present a molecular description of the SARS-CoV-2 spike protein for different omicron linages (BA.1, BA.2, BA.3, BA.4 and BA.). I do appreciate the effort by the author to provide molecular dynamics simulation data to compare all systems in terms of the structural flexibility and energetics. In this regard, I find the length of 200 ns and only one trajectory insufficient to validate their results. This number should be at least "in cumulative form" to reach 1 microsecond or perform several MD trajectories.

Before supporting further the study I provide major comments below.

(1) Please comment the recent article 10.1016/j.bbrc.2021.12.079, which provide a correlation between binding energy and the presence of larger buried surface area in omicron variants respect to wt.

(2) Figure 3: the analysis of the RMSD is subject to noise and so far I can only observe a reduce flexibility for the omicron variant respect to wt. If the author could add the error bar, using several trajectory or perform a block analysis, so we could observe major differences per residue.

(3) Binding free energies: here I have a comment, author has not reported the error bar in their analysis. The methodology employed (MM- PBSA) can be compare with other entries in the literature. The author have not described properly in the methodology this method. Are the average results depending on the selection of frames (i.e. 100 frames)?  Also, here better explanation of the systematic error is requested.

(4) Figure 4: This figure does not highlight the effect of the electrostatics. I see mean values for all cases are within the error bars. Thus, it is hard to asseverate that omicron favor electrostatics. Other argument by Podgornik et al. regarding electrostatics for omicron are requested to commend. 

(5) Moreover, I believe the authors do not provide at any instance a fairly comparison with experimental data. The dissociation constact (KD) has been calculated for several variants of concern by Koehler et al.  https://0-www-nature-com.brum.beds.ac.uk/articles/s41467-021-27325-1 and they can be commented. Does one expect Omicron to increase the binding affinity? In practice, Delta variant was pretty good in comparison with wt,

(6) Figure 8: report on the characterization of HBs. I believe contact map is better tool for this analysis. I suggest to explore a simply cut-off based or chemical based (i.e. OV+rCSU). See here http://pomalab.ippt.pan.pl/GoContactMap/. Additional information of the entire space of contact holding the interface could support for a better characterization as in ref. (https://0-doi-org.brum.beds.ac.uk/10.1021/acs.jcim.9b00883)

Author Response

Point-by-point response to the reviewer 1

REVIEWER REPORT(S):

The study by Sang et al present a molecular description of the SARS-CoV-2 spike protein for different omicron linages (BA.1, BA.2, BA.3, BA.4 and BA.5). I do appreciate the effort by the author to provide molecular dynamics simulation data to compare all systems in terms of the structural flexibility and energetics. In this regard, I find the length of 200 ns and only one trajectory insufficient to validate their results. This number should be at least "in cumulative form" to reach 1 microsecond or perform several MD trajectories.

Response:

Your comments are constructive and are greatly appreciated. In fact, when designing this study, we also considered using the multiple replica MD strategy, which can enhance the efficiency of conformational sampling and this method was often used in our former study. Considering that this study only focuses on thermodynamics rather than kinetics properties of the RBD-ACE2 binding, we finally decided to perform 200 ns single MD simulation on each complex. However, the simulation time is really too short for MM-PBSA calculation, and lack of statistical significance. We greatly appreciate and accept the invaluable suggestion from you, and we re-performed 10 independent 100-ns MD simulations on each complex, initializing MD runs with different initial atomic velocities taken from a Maxwell distribution of the corresponding temperature. The cumulative simulation time of each system reaches 1 microsecond. All the replicas then were used to perform the MM-PBSA calculation and other analyses, which can increase the statistical significance of the results. The results and conclusions of new trajectories, including MM-PBSA, are consistent with the previous manuscript. As a result, the error bar was added into Figure 3 in the revised manuscript, which can reflect systematic error due to different sampling.

Before supporting further the study I provide major comments below.

(1) Please comment the recent article 10.1016/j.bbrc.2021.12.079, which provide a correlation between binding energy and the presence of larger buried surface area in omicron variants respect to wt.

Response: Our conclusion is consistent with this article which indicated that omicron variants have lager buried surface area than WT, providing a correlation between binding energy and the presence of larger buried surface area. This article has been cited and commented in the discussion section of the revised manuscript.

(2) Figure 3: the analysis of the RMSD is subject to noise and so far I can only observe a reduce flexibility for the omicron variant respect to wt. If the author could add the error bar, using several trajectory or perform a block analysis, so we could observe major differences per residue.

Response: In contrast to MM-PBSA and non-bonded attractive interactions analysis, analysis of RMSF is actually not the most critical in our study. However, we must thank you again for your suggestion to perform multiple replica simulations. One single simulation trajectory may indeed not correctly reflect the real protein flexibility during simulation. In the revised manuscript, the error bar was added into new Figure 3. Figure 3 shows that those residues with larger standard errors are mainly located in the protein external loops, which themselves also have larger flexibility and RMSF values. Conversely, those residues with smaller standard errors themselves have smaller RMSF values ​​and are mainly located in the core region of the protein.

(3) Binding free energies: here I have a comment, author has not reported the error bar in their analysis. The methodology employed (MM- PBSA) can be compare with other entries in the literature. The author have not described properly in the methodology this method. Are the average results depending on the selection of frames (i.e. 100 frames)?  Also, here better explanation of the systematic error is requested.

Response: The stand error has been added into Table 1 in the revised manuscript. A more detailed description of MM-PBSA method and the comparison between MM-PBSA and other binding free energy calculation methods have been added into the materials and methods section in the revised manuscript. In MM-PBSA method, the average results are really depending on the selection of number of frames, and more frames can lead to better statistical significance of the calculation results. Since the MM-PBSA calculation is relatively time-consuming, we still evenly select 100 frames in each replica for calculation, and then take the average value of the 10 replicas (a total of 1000 frames were used for calculation) of each simulation system as the final result. Standard errors for each system are placed in parentheses after the final average value.

(4) Figure 4: This figure does not highlight the effect of the electrostatics. I see mean values for all cases are within the error bars. Thus, it is hard to asseverate that omicron favor electrostatics. Other argument by Podgornik et al. regarding electrostatics for omicron are requested to commend.

Response: Podgornik et al have provided excellent and detailed work in exploring the role of electrostatic interactions, expecially the partial charges, in determining binding affinity of RBD to ACE2. Their work has been commended and cited in our revised manuscript. Our conclusions for binding free energy (BFE) calculation and electrostatic interaction analysis are consistent with Podgornik et al that the enhanced electrostatic attractive interactions are the main determinant of the higher ACE2-binding affinity of Omicron RBD than WT RBD. We appreciate your review that Figure 4 does not highlight the effect of the electrostatics. We believe that the reason why there is no significant difference in the electrostatic interactions between Omicron and WT in Figure 4 is mainly because the gmx energy of gromacs can only provide the calculation results of short-range electrostatic interactions, so it cannot fully represent the overall electrostatic interactions. In order not to confuse readers, we have explained this situation in the revised manuscript.

(5) Moreover, I believe the authors do not provide at any instance a fairly comparison with experimental data. The dissociation constact (KD) has been calculated for several variants of concern by Koehler et al.  https://0-www-nature-com.brum.beds.ac.uk/articles/s41467-021-27325-1 and they can be commented. Does one expect Omicron to increase the binding affinity? In practice, Delta variant was pretty good in comparison with wt

Response: We appreciate your kindly suggestion on comparing with experimental results. Since the emergence of SARS-CoV-2, a large amount of experimental data on the binding affinity of RBD to ACE2 has been reported, and our computational results should indeed be compared with these experimental results. Koehler et al provided thermodynamic and kinetic insight into the binding of RBD and ACE2. There results shown that the RBDs of four VoCs (Alpha-RBD, Beta-RBD, Gamma-RBD, Kappa-RBD ) have higher ACE2-binding affinity (or lower KD) than WT. There work has been commented and discussed in the revised manuscript. Since the Delta and Omicron variants had not been presented and reported prior to the study of Koehler et al, we refer to some articles covering Omicron and Delta. Our computational results are consistent with these experimental and structural studies that Omicron RBD has a relatively higher ACE2-binding affinity than that of WT, but slightly lower than that of Delta（Han et al., Cell., 2022; Lan et al., Cell Research.,2022; Mannar et al., Science., 2022; Mccallum et al., Science., 2022）. Although the ACE2-binding affinity of Omicron RBD is lower than that of Delta, it still successfully replaced Delta as the dominant strain. As stated by Koehler et al, apart from ACE2-binding affinity, several other mechanisms might account for increased variant transmissibility, especially immune evasion. In addition, the experimental studies listed above indeed prove that the immune evasion ability of Omicron is indeed much higher than that of other mutants. In particular, our computational study show that although BA.4/5 have became the dominant strain in the world, their RBDs, however, have the lowest ACE2-binding affinity of all calculated Omicron subvariants. The results are consistant with Tuekprakhon et al and Cao et al, which shown that BA.4/5 evolved stronger immune evasion ability when compared with orther Omicron subvariants. The above discussion and references have been added to the revised manuscript.

(6) Figure 8: report on the characterization of HBs. I believe contact map is better tool for this analysis. I suggest to explore a simply cut-off based or chemical based (i.e. OV+rCSU). See here http://pomalab.ippt.pan.pl/GoContactMap/. Additional information of the entire space of contact holding the interface could support for a better characterization as in ref. (https://0-doi-org.brum.beds.ac.uk/10.1021/acs.jcim.9b00883)

Response: We appreciate your kindly suggestion about contact map. After reviewing the literatures, we found that contact map is a very good method to display the interactions between residues. The above two papers have been cited in the revised manuscript. Since the contact results between RBD and ACE2 provided by the online server did not show a clear trend between different virus strains, we put them in the supplementary material as Table S2-S6.

Reviewer 2 Report

The authors constructed homology models of RBD-ACE2 complexes of four Omicron subvariants (BA.1 BA.2, BA.3, and BA.4/5). They performed MD simulations on these complexes and the wild type of SARS-CoV-2 RBD complex with human ACE2 (hACE2). From these simulations, the authors obtained various structural dynamic properties of the binding free energy (BFE) between hACE2 and RBDwt and subvariants. The results suggest that the RBDs of all the Omicron 21 subvariants (RBDomis) feature decreased global structural fluctuations when compared with RBDwt. A comparison between the BFE components reveals that the enhanced electrostatic attractive interactions are the main determinant of the higher ACE2-binding affinity of RBDomis than RBDwt. The authors also found that the enhanced electrostatic attractive interactions are mainly through gain/loss of the positively/negatively charged residues, and the formation or destruction of the interface HBs and salt bridges can also largely affect the ACE2-binding affinity of RBDs.

The submitted research work is interesting since it points to several structural and energetic factors involved in the stronger binding of RBDomis to hACE2 versus RBDwt to hACE2. However, the authors did not present sufficient information as to how the molecular dynamics simulations were performed, nor did they discuss previous reports on the same subject, and how those findings are related to theirs.

The following questions must be answered before considering this manuscript for publication.

Page 13, line 38 - How does this work explain the recently observed action of a glycan to open the RBD of the Spike protein of SARS-CoV-2?. Sztain, T., Ahn, S. H., Bogetti, A. T., Casalino, L., Goldsmith, J. A., Seitz, E., ... & Amaro, R. E. (2021). A glycan gate controls opening of the SARS-CoV-2 spike protein. Nature chemistry, 13(10), 963-968.

Page14, line 39 – Structure Preparation --- The authors must mention if the RBDwt-ACE2 and RBDomis-ACE2 structures were protonated before performing the MD simulations. ---- The authors must mention if the PDB ID 6M0J structure is in the closed or in the open conformation. --- Are the RBDs glycosylated? If they are, what force field was used to treat the glycosylated amino acid residues? --- For all modeled RBD subvariants show the Ramachandran score of the percentage of amino acid residues in the various complexes of ACE2 that fall into the energetically favored region.

Page 14, line 40 – MD Simulation --- The authors performed a single run on each system. To avoid false positive conclusions, it has been advised to run several replicas of the simulations, so that the conclusions will be more reliable. The authors must perform three replicas, at the very least, since it is very important to show that different simulations converge to the same ensemble average properties. Knapp, B., Ospina, L., & Deane, C. M. (2018). Avoiding false-positive conclusions in molecular simulation: the importance of replicas. Journal of Chemical Theory and Computation, 14(12), 6127-6138. --- The authors must mention the net charge of the RBSwt-ACE2 and RBDomis-ACE2 complexes; and how they neutralized the system previous to the MD simulations. --- How many water molecules were added? What are the dimensions of the unit cell? --- How many minimization steps of steepest descent were performed? How often was the system sampled to perform BFE analysis? What algorithms were used to maintain temperature and pressure? How were the long-range electrostatics treated? How did they constrain all bonds involving hydrogen atoms?

Page 15, line 44 --- The authors must present a Table containing the following components: deltaG(binding), deltaG(complex), deltaG(ACE2) and deltaG(RBD); for all complexes.

Page 15, line 47 --- References --- The authors must discuss previous reports on the same subject by Jawad et al. (2022), Nguyen et al. (2022), and Dutta et al. (2022).

Jawad, B., Adhikari, P., Podgornik, R., & Ching, W. Y. (2022). Binding Interactions between Receptor-Binding Domain of Spike Protein and Human Angiotensin Converting Enzyme-2 in Omicron Variant. The journal of physical chemistry letters, 13, 3915-3921.

Nguyen, H. L., Thai, N. Q., Nguyen, P. H., & Li, M. S. (2022). SARS-CoV-2 Omicron Variant Binds to Human Cells More Strongly than the Wild Type: Evidence from Molecular Dynamics Simulation. The Journal of Physical Chemistry B, 126(25), 4669-4678.

Dutta, S., Panthi, B., & Chandra, A. (2022). All-Atom Simulations of Human ACE2-Spike Protein RBD Complexes for SARS-CoV-2 and Some of its Variants: Nature of Interactions and Free Energy Diagrams for Dissociation of the Protein Complexes. The Journal of Physical Chemistry B, 126(29), 5375-5389.

Author Response

Point-by-point response to the reviewer 2

REVIEWER REPORT(S):

   The authors constructed homology models of RBD-ACE2 complexes of four Omicron subvariants (BA.1 BA.2, BA.3, and BA.4/5). They performed MD simulations on these complexes and the wild type of SARS-CoV-2 RBD complex with human ACE2 (hACE2). From these simulations, the authors obtained various structural dynamic properties of the binding free energy (BFE) between hACE2 and RBDwt and subvariants. The results suggest that the RBDs of all the Omicron 21 subvariants (RBDomis) feature decreased global structural fluctuations when compared with RBDwt. A comparison between the BFE components reveals that the enhanced electrostatic attractive interactions are the main determinant of the higher ACE2-binding affinity of RBDomis than RBDwt. The authors also found that the enhanced electrostatic attractive interactions are mainly through gain/loss of the positively/negatively charged residues, and the formation or destruction of the interface HBs and salt bridges can also largely affect the ACE2-binding affinity of RBDs.

The submitted research work is interesting since it points to several structural and energetic factors involved in the stronger binding of RBDomis to hACE2 versus RBDwt to hACE2. However, the authors did not present sufficient information as to how the molecular dynamics simulations were performed, nor did they discuss previous reports on the same subject, and how those findings are related to theirs.

The following questions must be answered before considering this manuscript for publication.

Responce: We highly appreciate your affirmation of our work. We apologize for insufficiently description of the some methods and results in previous manuscript. In the revised manuscript, the descriptions are more detailed and rigorous. Your kindly provided literatures are cited and discussed in the revised manuscript.

• Page 13, line 38 - How does this work explain the recently observed action of a glycan to open the RBD of the Spike protein of SARS-CoV-2?. Sztain, T., Ahn, S. H., Bogetti, A. T., Casalino, L., Goldsmith, J. A., Seitz, E., ... & Amaro, R. E. (2021). A glycan gate controls opening of the SARS-CoV-2 spike protein. Nature chemistry, 13(10), 963-968.

Responce: We would first like to thank you for providing us such a high level of paper, and the paper has been discussed and cited in the revised manuscript. Although the erection of the RBD is a prerequisite for ACE2 binding, the RBD is an independently folded domain and its erection has only a marginal impact on its overall conformation (Benton et al., Nature., 2020; Zhang et al., Cells., 2022). In this study, the binding affinity of the spike protein to the ACE2 was evaluated using the RBD rather than the spike trimer. Therefore, our study may not directly explain the action of a glycan to open the RBD of the Spike protein of SARS-CoV-2. However, our per-residue BFE calculation shows that residue N343 has no significant positive or negative contribution to the final BFE. The reason for this may be that it is protected by glycosylation, which reduces various interactions with ACE2 residues. In addition, the mutation of N343 site has never appeared in any Omicron variant, which also reflects the high degree of conservation of glycosylated N343.

• Page14, line 39 – Structure Preparation --- The authors must mention if the RBDwt-ACE2 and RBDomis-ACE2 structures were protonated before performing the MD simulations. ---- The authors must mention if the PDB ID 6M0J structure is in the closed or in the open conformation. --- Are the RBDs glycosylated? If they are, what force field was used to treat the glycosylated amino acid residues? --- For all modeled RBD subvariants show the Ramachandran score of the percentage of amino acid residues in the various complexes of ACE2 that fall into the energetically favored region.

Responce: The RBDwt-ACE2 and RBDomis-ACE2 structures were protonated according to the protonation states of all titratable residues at pH 7.4, and we have mentioned this in the revised manuscript. Since the authors of 6M0J structure only expressed the SARS-CoV-2 RBD and the N-terminal peptidase domain of ACE2, their structure do not contain information about whether the spike is in closed or in the open conformation. Nevertheless, this structure was widely used in the studies on RBD-ACE2 interactions. In our study, the RBDs are not glycosylated. The Ramachandran plots of the complexes were put into Figure S3-S7, and the analyses of Ramachandran score were performed in the revised manuscript.

(3) Page 14, line 40 – MD Simulation --- The authors performed a single run on each system. To avoid false positive conclusions, it has been advised to run several replicas of the simulations, so that the conclusions will be more reliable. The authors must perform three replicas, at the very least, since it is very important to show that different simulations converge to the same ensemble average properties. Knapp, B., Ospina, L., & Deane, C. M. (2018). Avoiding false-positive conclusions in molecular simulation: the importance of replicas. Journal of Chemical Theory and Computation, 14(12), 6127-6138. --- The authors must mention the net charge of the RBSwt-ACE2 and RBDomis-ACE2 complexes; and how they neutralized the system previous to the MD simulations. --- How many water molecules were added? What are the dimensions of the unit cell? --- How many minimization steps of steepest descent were performed? How often was the system sampled to perform BFE analysis? What algorithms were used to maintain temperature and pressure? How were the long-range electrostatics treated? How did they constrain all bonds involving hydrogen atoms?

Responce: We greatly appreciate and accept the invaluable suggestion from you, and we re-performed 10 independent 100-ns MD simulations on each complex, initializing MD runs with different initial atomic velocities taken from a Maxwell distribution of the corresponding temperature. The cumulative simulation time of each system reaches 1 microsecond. All the replicas then were used to perform the MM-PBSA calculation and other analyses, which can increase the statistical significance of the results. The results and conclusions of new trajectories, including MM-PBSA, are consistent with the previous manuscript. The paper by Knapp et al has been cited in the revised manuscript. We apologize for the imprecise description of the MD simulation method. The quantity information you mentioned above we have added to the revised manuscript as Table S7. In our MD simulation, the temperature and pressure were maintained by v-rescale thermostat and Parrinello–Rahman barostat, respectively. The long-range electrostatics were treated by the Particle-mesh Ewald (PME) algorithm. The bonds, including hydrogen atoms were constrained by LINCS algorithms with the integration time step of 2 fs.

(4) Page 15, line 44 --- The authors must present a Table containing the following components: deltaG(binding), deltaG(complex), deltaG(ACE2) and deltaG(RBD); for all complexes.

Responce: These components were shown in Table S1 in the revised manuscript.

(5) Page 15, line 47 --- References --- The authors must discuss previous reports on the same subject by Jawad et al. (2022), Nguyen et al. (2022), and Dutta et al. (2022).

(6) Jawad, B., Adhikari, P., Podgornik, R., & Ching, W. Y. (2022). Binding Interactions between Receptor-Binding Domain of Spike Protein and Human Angiotensin Converting Enzyme-2 in Omicron Variant. The journal of physical chemistry letters, 13, 3915-3921.

(7) Nguyen, H. L., Thai, N. Q., Nguyen, P. H., & Li, M. S. (2022). SARS-CoV-2 Omicron Variant Binds to Human Cells More Strongly than the Wild Type: Evidence from Molecular Dynamics Simulation. The Journal of Physical Chemistry B, 126(25), 4669-4678.

(8) Dutta, S., Panthi, B., & Chandra, A. (2022). All-Atom Simulations of Human ACE2-Spike Protein RBD Complexes for SARS-CoV-2 and Some of its Variants: Nature of Interactions and Free Energy Diagrams for Dissociation of the Protein Complexes. The Journal of Physical Chemistry B, 126(29), 5375-5389.

Response: Thanks for providing us the above three papers for our comparison and discussion. The three pepers have been discussed and cited in the revised manuscript. The results and conclusions of our study are consistent with theirs that the ACE2-binding affinity of Omicron RBD is higher than that of WT SARS-CoV-2, and the electrostatic interactions play decisive role in determining the RBD-ACE2 binding affinity. However, unlike they selected only one Omicron subvariant for their studies, we studied and compared the differences in ACE2-binding affinity of 4 subvariant RBD. Our results indicate that although Omicron BA.4/5 have replaced the previous subvariants as the predominant SARS-CoV-2 strain, their RBDs have the lowest ACE2-binding affinity of all. Combined with our conclusion, we infer that the high transmissibility of new Omicron subvariant may primarily originate from the enhanced immune evasion rather than the higher ACE2-binding affinity.

Round 2

Reviewer 1 Report

I am

glad to endorse the article for publications.

Reviewer 2 Report

The authors have improved their original version. Performing several MD replicas now validates their conclusions. Furthermore, the write-up is much better, and readers can read the paper pleasantly.