The pH Effects on SARS-CoV and SARS-CoV-2 Spike Proteins in the Process of Binding to hACE2.

COVID-19 has been threatening human health since the late 2019, which has significant impact on human health and economy. Understanding the SARS-CoV-2 and other coronaviruses is important to develop effective treatments for COVID-19 and other coronaviruses-caused diseases. In this work, we applied multi-scale computational approaches to study the electrostatic features of spike (S) proteins for SARS-CoV and SARS-CoV-2. From our results, we found thatSARS-CoV and SARS-CoV-2 have similar charge distributions and electrostatic features when binding with the human angiotensin-converting enzyme 2 (hACE2). The energy pH-dependence calculation srevealed that the complex structures of hACE2 and the S proteins of SARS-CoV/SARS-CoV-2 are stable at pH values ranging from 7.5 to 9. Molecular dynamics simulations were performed using NAMD to investigate the hydrogen bonds between S proteins and hACE2. From the MD simulations it was found that SARS-CoV-2 has four pairsof essential hydrogenbonds (high occupancy, >80%), while SARS-CoV has three pairs, which indicates the SARS-CoV-2 S protein has relatively more robust binding strategy than SARS-CoVS protein.Four key residues forming essential hydrogen bonds from SARS-CoV-2 are identified, which are potential drug targets for COVID-19 treatments. The findings in this study shed lights on the current and future treatments for COVID-19 and other coronaviruses-caused diseases.

Introduction

The ongoing COVID-19 pandemic is changing human society significantly and causing both economic and social consequences all over the world [1]. Coronaviruses are named for their crown-like spikes on their surface, and they are commonly found in many mammal species [2]. Human coronaviruses were firstly identified in the mid-1960s. There are four main sub-groupings of coronaviruses, known as alpha, beta, gamma, and delta [3]. Among all the coronaviruses, there are seven known types of coronaviruses that can infect human beings. People around the world commonly get infected by human coronaviruses 229E, NL63, OC43, and HKU1 [4, 5]. And some coronaviruses that infect animals are able to evolve and infect humans, among which the three recent cases are SARS-CoV-2, SARS-CoV, and MERS-CoV[6]. The SARS-CoV-2 virus is the novel coronavirus that causes coronavirus disease 2019, or COVID-19. Other than COVID-19, coronaviruses have caused several pandemics before, including severe acute respiratory syndrome (SARS) which was caused by SARS-CoV and the Middle East respiratory syndrome (MERS) which was caused by MERS-CoV. To end the current pandemic soon and be prepared for the future similar challenges for human society, it is essential to understand the binding mechanisms of SARS-CoV-2 infecting human cells. This is achievable by studying the stabilities of SARS-CoV-2 at different pH conditions, and identify the key residues that play significant roles in the binding processes.

Coronaviruses contain membrane glycoprotein (M), nucleocapsid protein (N), spike protein (S), envelope protein (E) and an RNA single chain[7]. For all enveloped viruses, one of the most important steps during the binding process is membrane fusion, which allows viruses to get into host cells [8]. For coronaviruses, the fusion protein is the S protein that leads the binding process to attack human cells through the host cell receptor angiotensin-converting enzyme 2 (hACE2) [9]. Human hACE2 (hACE2) is an enzyme located widely in the human body, including the lungs, kidneys, adipose tissue, central nervous system and cardiovascular system [9-11] and it has multiple essential functions such as the regulation of amino acid transport in the kidney controlling the blood pressure, and viral receptors including both SARS-CoV-2 and SARS-CoV [11]. Since it is of extreme importance to human health, there are numerous research groups have been or are currently working on S proteins and hACE2 using various approaches.

The traditional process of the de novo drug design is a challenging task which consumes resources and time significantly. With the fast developments of computing technology, computational methods have been widely used to drug-related research[12], including protein-protein interactions[13, 14], MD simulations[15], coarse-grained models[16], pH dependence of protein-protein interactions[17-20], etc. Our previous studies have applied multi-scale computational methods to study several pathogens [21-25] including the SARS-CoV-2 viruses [26, 27], which revealed some mechanisms of the SARS-CoV-2 S protein. Besides, many other research groups have made successful progress to understand the SARS-CoV-2 using computational methods [28, 29].

In this work, we first calculated the electrostatic potentials on the surface of S proteins from both SARS-CoV and SARS-CoV-2, followed by the electric field line comparison between SARS-CoV and SARS-CoV-2 when they bind to hACE2. We found that the two viruses have similar pH responses: The pH-dependence of folding energies for S protein receptor binding domains (RBDs) demonstrated that both the S protein RBDs of these two viruses are at the most stable status when pH values ranging from 6 to 9. Also, the pH-dependence of binding energies for S protein RBDs and hACE2 RBD showed that the complex structures of the two viruses are at the most stable status at pH values ranging from 7.5 to 10.5. Therefore, SARS-CoV and SARS-CoV-2 survive in a similar pH environment. The pH 7.5 to 9 is the best condition for both SARS-CoV and SARS-CoV-2 to best perform their functions to bind with hACE2. Also, we analyzed the trajectories from 100ns MD simulations using NAMD [30] and identified hydrogen bonds with the involved key residues using VMD [31]. It is shown that for the high-frequency (>80%) hydrogen bonds, SARS-CoV-2 has four pairs while SARS-CoV has three pairs, which indicates that the S protein of SARS-CoV-2 uses more residues to form strong hydrogen bonds. The key residues forming essential hydrogen bonds from SARS-CoV-2 are ARG-121, TYR103, THR182 and TYR171, which are potential drug targets for COVID-19 treatments. Using multiple computational approaches, the findings in this work pave the way for the current and future treatment development of COVID-19 and other coronaviruses-caused diseases.

Methods

Structure Preparation

The complex structures of SARS-CoV/hACE2 and SARS-CoV-2/hACE2 were downloaded from the Protein Data Bank (PDB ID 6ACG [32] and 7AD1 [33], respectively) . Please note that in 7AD1, the mutations that the authors made during their experiments are not on the interface area. Since we only focus on the interface area between S proteins and hACE2, the mutations do not affect our results. In this work, we used the complex structures to study the electrostatic binding interactions and the relative binding energies in different pH environments between S proteins and hACE2 RBDs. For the missing loops in proteins, we used MODELLER [34] to model the structures based on the sequences. To understand the mechanisms of S protein binding to hACE2 at the interface particularly, S protein RBDs were separated from the hACE2 binding domain by a distance of 10Å for the best results and visualization.

Electrostatic Potential Calculation

In order to study the electrostatic features, DelPhi [35, 36] was utilized to calculate the electrostatic potential for the S proteins and hACE2 RBDs. In the framework of continuum electrostatics, DelPhi calculates the electrostatic potential ϕ (in systems comprised of biological macromolecules and water in the presence of mobile ions) by solving the Poisson-Boltzmann equation (PBE): where ϕ(r) is the electrostatic potential, ϵ(r) is the dielectric distribution, ρ(r) is the charge density based on the atomic structures, κ is the Debye-Huckel parameter, kB is the Boltzmann constant, and T is the temperature. Due to the irregular shape of macromolecules, DelPhi uses a finite difference (FD) method to solve the PBE.

Before the DelPhi calculations, the PQR file of each trimer was generated by PDB2PQR [37]. We used AMBER [38]force field for PDB2PQR calculation, and removed water molecules. For the better results, we ensured the new atoms are not rebuilt too close to existing atoms and optimized the hydrogen bonding network.

During DelPhi calculations, the resolution was set as 0.5 grids/Å. The dielectric constants were set as 2.0 for protein and 80.0 for the water environment, respectively. The pH value for the solvent environment was set to be 7.0. The probe radius for generating the molecular surface was 1.4 Å. Salt concentration was set as 0.15 M. The boundary condition for the Poisson Boltzmann equation was set as a dipolar boundary condition. The calculated electrostatic potential on the surface was visualized with Chimera (figure 2). VMD was used to illustrate electric field lines between S protein and hACE2 (figure 3). Finally, the color scale range was set to be from −1.0 to 1.0 kT/e for the best visual presentation. Besides the calculations of electrostatic potentials, we also used DelphiForce [39] to calculate the electrostatic binding forces between each S protein and hACE2 while separating them in the direction of the mass center connection line (figure S2). Besides the net forces between each S protein and hACE2, the X, Y, Z components of the net forces are also calculated and shown in figure S2.

Figure 2.

Electrostatic potential on surfaces of SARS-CoV and SARS-CoV-2 S proteins. (A) Top view of S protein structure; (B-C) Top views of electrostatic potential on surfaces of SARS-CoV and SARS-CoV-2 S protein, respectively; (D) Front view of S protein structure; (E-F) Front views of electrostatic potential on surfaces of SARS-CoV and SARS-CoV-2 S protein, respectively; (G) Bottom view of S protein structure; (H-I) Bottom views of electrostatic potential on surfaces of SARS-CoV and SARS-CoV-2 S protein, respectively. Negatively and positively charged areas are colored in red and blue respectively, with the color scale from −1.0 to 1.0 kT/e.

Figure 3.

Electrostatic filed lines at the interfaces of S protein and hACE2. (A) Electrostatic filed lines between SARS-CoV S protein and hACE2; (B) A closeup view of binding domain between SARS-CoV S protein and hACE2 (C) Electrostatic field lines between SARS-CoV-2 S protein and hACE2; (D) A closeup view of binding domain between SARS-CoV-2 S protein and hACE2. Negatively and positively charged areas are colored in red and blue, respectively. Color scale is −1.0 to 1.0 kT/e. Yellow square areas are the RBD of S proteins at open state to reach the hACE2, cyan square areas are the the RBD of S proteins at closed state.

Relative Folding Energy Calculation

We used DelPhiPKa [40, 41] to calculate pKa values of DNA and UDG, given the pH ranging from 0 to 14 with the pH interval of 0.5. During the calculations, we used AMBER force field, and removed water molecules and HETATM. For the hydrogen of ASP and GLU attached atom, we used OD1 and OE1, respectively. Variance of Gaussian Distribution was set to be 0.7, salt concentration was 0.15, reference dielectric was 8.0, and external dielectric was 80.0.

The net charges of proteins at the unfolded state were calculated using this equation: where the summation is of all the titratable groups, y(i) value is −1 for acidic groups and +1 for basic groups, respectively. As for the folding free energy, we used this equation: where Q (pH) and Q𝑢(pH) stand for the net charge of folded and unfolded state, respectively. R is the universal gas constant taken as . T is the temperature with the value of 300 K.

Please note that the algorithms we applied to calculate the folding energies are for the relative values, that is, at pH=0 the folding energy is 0 and at any other pH values the folding energies are the relative values to the pH=0 condition.

Relative Binding Energy Calculation

For the binding energy calculation, we involved two methods, which are DelPhiPKa and MM/PBSA[42]. To calculate binding energy using DelPhiPKa, the following equation was used: where ΔN(pH) is the the binding free energy at different pH values, Q(pH), Q(pH) , and Q𝑟(pH) are the net charges of complexes of each model. R is the universal gas constant taken as . T is the temperature with the value of 300 K.

Please note that the algorithms we applied to calculate the binding energies are for the relative values, that is, at pH=0 the binding energy is 0 and at any other pH values the binding energies are the relative values to the pH=0 condition.

Molecular Dynamic (MD) Simulations

To simulate the dynamic interactions between S proteins RBD and hACE2 protein, MD simulations [15] were carried out using NAMD [30] with the help of GPUs on Lonestar5 clusters at the Texas Advanced Computing Center (TACC https://www.tacc.utexas.edu/). A 2000-step minimization was performed for each simulation, followed by a 100 million steps, during which 20,000 frames were saved from two 100ns simulations of both SARS-CoV and SARS-CoV-2 separately (1.0 fs per step, 1 frame at each 5000 steps, 100 million steps in total). The RMSDs of the SARS-CoV and SARS-CoV-2 trajectories are about 3.4Å and 1.1 Å, respectively (figure S1). During the MD simulations, we used CHARMM [43] force field, the temperature was set to be 300 K, and the pressure was set to be standard using the Langevin dynamics. For PME, which is set for full-system periodic electrostatics, with the grid size (86, 88, 132) as (x, y, z) value respectively. In those two simulations, atoms that are not located in binding domains were constrained within a margin of 10.0 Å of their natural movement maximum length values. In order to get a more accurate result of the simulation, data of the last 50 ns of simulations were selected and used for data analysis, since the structure of the first 50 ns is not as stable as the last 50 ns of simulations. The simulation processes are visualized in movies 4 and 5, generated by VMD.

To analyze the interaction between S proteins and hACE2, the hydrogen bonds that formed within the distance of 4 Å were extracted from the last 10,000 frames (50 ns) of simulations. The several top-strongest hydrogen bonds in each binding domain were determined by calculating their formation frequency (the frequency in figure of essential hydrogen bonds is shown in figure 7.

Figure 7.

Essential Hydrogen bonds at the interfaces between SARS-CoV/SARS-CoV-2 RBDs and hACE2 RBD with the frequency above 80%.

Results and Discussions

First of all, the electrostatic features of SARS-CoV and SARS-CoV-2 S proteins were investigated, including electrostatic potential and electric field lines. Secondly, the relative binding energies of complex structures and folding energies of S proteins at different pH values were analyzed. Finally, the hydrogen bonds and related key residues in each complex structure were obtained using MD simulations.

S Protein Trimer Structure

The RMSD between the S proteins of SARS-CoV and SARS-CoV-2 is 0.973 Å, showing that the S proteins of SARS-CoV and SARS-CoV-2 are very similar. The S proteins of SARS-CoV and SARS-CoV-2 are both homotrimers. Each monomer contains an RBD which connects the other part of the monomer via a hinge composed by two flexible loops (as shown in the black circle of figure 1A). The RBD is in closed configuration when there is no hACE2 binds to the S protein. When binding to hACE2, the RBD of one monomer flips out as open configuration and it binds to the RBD of hACE2.

Figure 1.

SARS-CoV S protein structure. Only the SARS-CoV S protein structure is illustrated in this figure, because SARS-CoV and SARS-CoV-2 S proteins are very similar (the RMSD between two S protein RBDs is 0.973 Å). (A) The S protein is a homotrimer (orange, blue, pink), of which one chain (pink) flips out when it binds to hACE2 (gray). The hinge connecting the RBD and the other part of S protein is shown in a black circle; (B) The closeup view of binding domains when S protein RBD (pink) binds to hACE2 RBD (gray).

Electrostatic Potential on Surfaces

To study the electrostatic features, DelPhi was utilized to calculate the electrostatic potential on surfaces of the S protein trimer (full structure) and hACE2 RBD. The electrostatic potential distribution on SARS-CoV S protein trimer structure is showed in figure 2BEH and movie 1, which were rendered by Chimera with a color scale from −1.0 to 1.0 kT/e. The charge distribution on SARS-CoV-2 S protein trimer structure is shown in figure 2CFI and movie 2, which were rendered by Chimera with a color scale from −1.0 to 1.0 kT/e as well, for the comparison. Negatively and positively charged areas are colored in red and blue, respectively.

By comparing the electrostatic potential on surfaces of two trimer structures, it is obvious that the charge distribution of SARS-CoV and SARS-CoV-2 S proteins are different. From the top view (figure 2A–C) and the bottom view (figure 2G–I), we noticed that SARS-CoV has slightly more positively charged area (blue), compared to SARS-CoV-2. It indicates that the SARS-CoV may attract the hACE2 more easily, since the hACE2 binding interface is overall negatively charged (movie 3). Such finding supports the previous studies of our research group [26, 27]. The electrostatic distribution differences observed from front views (figure 2D–F) of the S proteins demonstrate that the electrostatic features may have impacts on the stabilities of the trimers. Here it was not investigated several details about the binding stabilities among monomers in an S protein, due to the scope of this work that mainly focusses on the binding between S protein and hACE2. The electrostatic distributions on S protein RBDs show that the SARS-CoV RBD is more positive, which is consistent with the top view (figure 2BC). The bottom of the SARS-CoV (Figure 2EH) has more positive potential than SARS-CoV-2 (figure 2FI).

Electric Filed Lines

Electric field lines surrounding the two complex structures were calculated. To better visualize the field lines between interfaces, the S protein RBDs are separated from hACE2 RBDs by 10Å (figure 3). The field lines confirmed that both the SARS-CoV and SARS-CoV-2 S protein RBDs have attractive forces to hACE2 protein. In the analysis of field lines, the density of field lines indicates the strength of binding force, which means the denser area has the stronger interactions. The electric field lines demonstrate that when hACE2 is away from S protein, all the three S protein monomers provide attractive interactions to the hACE2. This is expected because the S protein RBDs are positively charged while the hACE2 is negatively charged, as shown in figure 2 and movie 3, respectively. When hACE2 binds to S proteins (as shown in figure 1), the hACE2 only binds with one S protein RBD, which is in open state. Combining the information from figure 1 and 3, it demonstrates that all the three S protein RBDs generate attractive forces to hACE2. However, when hACE2 gets closer to S protein, one S protein RBD flips out and binds to the hACE2 tightly, while the other two S protein RBDs stay in closed state. Even though the monomer with flipped-out S protein RBD is the closest to hACE2 and forms most of the salt bridges and hydrogen bonds, the other two monomers also provide dense field lines and show strong attractive interactions between S proteins and hACE2.

pH-Dependence of Relative Folding Energies

The folding energy of SARS-CoV and SARS-CoV-2 complexes were calculated using DelPhiPKa at different pH values ranging from 0 to 14 with an interval of 0.5 (Figure 4). We observed that SARS-CoV and SARS-CoV-2 have the same trend of folding energy with the change of pH values, which is decreasing from 0 to 6, then becoming stable from 6 to 9, and increasing from 10 to 14. Other than the trend, the optimal values locate between 6 to 9 for both of the viruses.

Figure 4.

pH-dependence of the relative folding energy of S protein RBDs of SARS-CoV and SARS-CoV-2.

Please note that the folding energies in figure 4 are relative values because we set the reference energy to be 0 kcal/mol when pH is equal to 0. We did not calculate the absolute values of folding energies since we focused on the pH dependency of the folding energies.

pH-Dependence of Relative Binding Energies

DelPhiPKa was implemented to calculate the binding energies of two complex structures at different pH values. The results are presented in figure 5, where we noticed that the binding free energies of both SARS-CoV and SARS-CoV-2 complexes are stable at the pH values ranging from 7.5 to 10.5, which indicates that both SARS-CoV and SARS-CoV-2 have a slight preference of weakly basic environment. Note that the method implementing DelPhiPKa calculates the relative folding and binding energies rather than absolute energies. The folding/binding energy at pH 0 is set as reference, which is 0 kcal/mol. The relative energy profile is used to study the folding/binding energy dependence on pHs. The absolute binding energies was calculated in later section using MM/PBSA method. Combine the folding and binding energy profiles, it is concluded that the best pH environment for both the SARS-CoV and SARS-CoV-2 is from pH 7.5 to 9. Please note that the binding energies in figure 5 are relative values because we set the reference energy to be 0 kJ/mol when pH is equal to 0. We did not calculate the absolute values of binding energies since we focused on the pH dependency of the binding stability.

Figure 5.

The relative binding energies of complexes at different pH values.

Hydrogen Bonds Analysis

To analyze the hydrogen bonds distributions on both S proteins RBDs and hACE2 RBD, we colored the residues forming hydrogen bonds which are over 50% frequency during the MD simulations in figure 6. It’s obvious that the SARS-CoV S protein has more residues involved in the hydrogen bonds which are over 50%. Accordingly, the hACE also has more residues forming hydrogen bonds (over 50% frequency) with SARS-CoV S protein.

Figure 6.

Hydrogen bonds distributions at the binding interfaces. (A) Hydrogen bonds distribution (blue) on the interface of SARS-CoV RBD (orange); (B) Turn (A) for 90 degree for the top view, which is the interface that faces hACE2; (C) The hydrogen bonds distribution (pink) at the interface of hACE2 (grey) where SARS-CoV binds; (D) Hydrogen bonds distribution (blue) on the interface of SARS-CoV-2 RBD (orange); (E) Turn (D) for 90 degree for the top view, which is the interface that faces hACE2; (F) The hydrogen bonds distribution (pink) at the interface of hACE2 (grey) where SARS-CoV-2 binds.

In order to consider the most essential hydrogen bonds, which are the hydrogen bonds with relatively high frequencies, we took 80% as a cutoff, which means those hydrogen bonds with 80% or higher frequency are considered as the relatively more essential ones. By comparing the figure 7A and 7B, SARS-CoV-2 RBD forms one more essential hydrogen bonds than SARS-CoV RBD when binding to hACE2. The residues involved in forming hydrogen bonds over 50% frequency were colored with their side chains, in which the residues with over 80% frequency hydrogen bonds were labeled and highlighted in grey squares (figure 8CF). From the analyses of figure 6–8, it is revealed that SARS-CoV uses more hydrogen bonds to bind with hACE2. However, more high frequency hydrogen bonds are formed in the SARS-CoV-2/hACE2 complex. The key residues forming essential hydrogen bonds from SARS-CoV-2 are: ARG-121, TYR103, THR182 and TYR171. Such residues have significant contributions to the binding of SARS-CoV-2 and hACE2. Therefore, these residues have higher potential to be targets for future drug design.

Figure 8.

Key residues involved in essential hydrogen bonds at the interfaces between SARS-CoV/SARS-CoV-2 RBDs and hACE2 RBD with the frequency above 80%. (A) SARS-CoV S protein single chain binds to hACE2; (B) A closeup view of (A) at the binding interface; (C) Labelled key residues that form essential hydrogen bonds (frequency over 80%) at the interface; (D) SARS-CoV-2 S protein single chain binds to hACE2; (E) A closeup view of (D) at the binding interface; (F) Labelled key residues that form essential hydrogen bonds (frequency over 80%) at the interface.

Limitation

The limitation for this work is that we used relative folding energy and binding energy to analyze rather than the absolute values. Since our work is focused on the relative stability under the pH effects, the relative energy calculations do not affect our conclusions.

Conclusion

In this work, we applied several computational methods, including MD simulations, DelPhi, DelPhiForce and DelPhiPKa to study the electrostatic features of S proteins for SARS-CoV and SARS-CoV-2. From our results, SARS-CoV and SARS-CoV-2 S protein RBDs both have positively charged interfaces, which provides attractive interactions to hACE2 as hACE2 has negatively charged surface.

Also, we revealed the pH-dependence calculations of relative folding energy for SARS-CoV and SARS-CoV-2 S protein RBDs. The best pH to stabilize SARS-CoV and SARS-CoV-2 S protein RBDs is in the range of 6 to 9. The study on pH dependence of binding energies revealed that the complex structures of hACE2 and S proteins of SARS-CoV/ SARS-CoV-2 are stable from pH 7.5 to 10.5. Therefore, SARS-CoV and SARS-CoV-2 survive in a similar pH environment. The pH 7.5 to 9 is the best condition for both SARS-CoV and SARS-CoV-2 to best perform their functions to bind with hACE2.

Besides, based on 100ns MD simulations, we found that for the essential hydrogen bonds (>80% frequency), SARS-CoV-2 has four pairs while SARS-CoV has three pairs, which indicates the relatively more robust binding strategy of SARS-CoV-2 compared to SARS-CoV. The key residues forming essential hydrogen bonds from SARS-CoV-2 are ARG-121, TYR103, THR182 and TYR171, which are potential drug targets for COVID-19 treatments. By using multiple computational approaches, the findings in this work shed light on the current and future treatments of COVID-19 and other coronaviruses-caused diseases.