Ionization of D571 Is Coupled with SARS-CoV-2 Spike Up/Down Equilibrium Revealing the pH-Dependent Allosteric Mechanism of Receptor-Binding Domains.

As a type I viral fusion protein, SARS-CoV-2 spike undergoes a pH-dependent switch to mediate the endosomal positioning of the receptor-binding domain to facilitate viral entry into cells and immune evasion. Gaps in our knowledge concerning the conformational transitions and key intramolecular motivations have hampered the development of effective therapeutics against the virus. To clarify the pH-sensitive elements on spike-gating the receptor-binding domain (RBD) opening and understand the details of the RBD opening transition, we performed microsecond-time scale constant pH molecular dynamics simulations in this study. We identified the deeply buried D571 with a clear pKa shift, suggesting a potential pH sensor, and showed the coupling of ionization of D571 with spike RBD-up/down equilibrium. We also computed the free-energy landscape for RBD opening and identified the crucial interactions that influence RBD dynamics. The atomic-level characterization of the pH-dependent spike activation mechanism provided herein offers new insights for a better understanding of the fundamental mechanisms of SARS-CoV-2 viral entry and infection and hence supports the discovery of novel therapeutics for COVID-19.

Introduction

The outbreaks of coronavirus disease 2019 (COVID-19) caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) have resulted in an unprecedented global health crisis,[1] and to date, very few effective therapeutics are available for the treatment of the fatal contagious disease.[2] Meanwhile, the fast evolution of the virus has caused stronger infection and toxicity,[3−7] posing great challenges for combating against the virus.

The main infection machinery of the virus relies on its envelope spike glycoprotein to mediate the viral entry into host cells.[8] Thus, a particularly promising target in the viral life cycle for therapeutic design is the spike protein, and currently, it has been taken as the main target for antibody responses against the virus.[9−12] The SARS-CoV-2 spike glycoprotein is a homotrimeric type I fusion protein, as shown in Figure A,B, with each protomer consisting of two subunits, S1 and S2, in which the upper S1 subunit is responsible for recognizing and binding the host cell receptor angiotensin-converting enzyme 2 (ACE2) and the lower S2 subunit contains the fusion machinery.[13−20]

Figure 1

(A) Structures of the spike homotrimer with RBDs in all-down and single-up positions. Three protomers are encoded with violet, light blue, and wheat colors, respectively. (B) Color-encoded domain organization of the spike protomer and the location of internal D571. (C) Spike RBD conformational distribution in terms of the opening angle and rotation for the experimentally determined structures summarized from the protein data bank.

Like other coronaviruses, SARS-CoV-2 infection primarily utilizes an endosomal entry pathway with the pH change of the milieu within 4.5–6.0, along which the receptor-binding domain (RBD) in the S1 subunit undergoes up/down structural transitions (see Figure A,B) by sensing environmental pH changes.[16,21−27] The cryo-electron microscopy (cryo-EM) structures have revealed two prevalent prefusion conformations of SARS-CoV-2 spike protein, single-up and all-down conformations, related to the RBD positioning.[28−31] The RBD single-up state is predominantly populated at serological pH (>5.5), enabling the virus to recognize and bind the ACE2 receptor.[21,32] The up RBD is also related to the epitope availability of RBD-directed antibodies.[31,33,34] While in late-endosome–early-lysosome pH (<5.5), the all-RBD-down conformation is prevalent, potentially inducing immune evasion from RBD-up-recognizing antibodies.[21] This pH-sensitive switching mechanism of SARS-CoV-2 spike protein, which facilitates the viral entry into host cells and self-protection against host immune responses in the internal physiological environment, has been widely observed in many other viral fusion proteins, such as vesicular stomatitis virus glycoprotein,[35−37] influenza virus hemagglutinin,[38] rabies virus glycoprotein,[39] Ebola virus envelope glycoprotein,[40] and other enveloped virus surface glycoproteins.[41] Therefore, we have to ask the critical questions: what are the pH-sensitive elements gating the spike RBD opening? Where is the RBD transition pathway? Answering these questions would provide crucial information for a better understanding of the detailed course of viral infection and, more importantly, help to discover new therapeutics for not only COVID-19 but also the other related fatal viruses’ infectious diseases.

At present, the precise molecular mechanism governing RBD positioning is still being delineated.[42,43] To this end, both experimental and theoretical communities have made great effort to build models of SARS-CoV-2 spike protein as well as its complex form with ACE2.[21,42−47] It is well known that current experimental techniques are normally incapable of accessing the full range of large protein conformational changes and only limited to a few local energy-minimum structures. Moreover, it is very difficult to experimentally determine the ionization states of the ionizable residues directly, especially for those buried in the hydrophobic interior of proteins, which can function as pH sensors and trigger the protein structural reorganization of varying amplitudes.[48−52] Theoretical molecular dynamics (MD) simulations have the capability to capture a broader ensemble of snapshots connecting the various conformational states that a protein can visit and have greatly contributed to our better understanding of the dynamics of the spike protein.[43] Casalino et al.’s conventional MD simulations indicated the mechanical role of spike surface glycans in affecting the RBD-up/down balance.[53] Sztain et al. utilized a path-sampling strategy to simulate the RBD opening pathway and also revealed the role of a glycan gate.[43] Zimmerman et al. discovered dramatic conformational changes of spike protein through long-time MD simulations with the adaptive sampling algorithm.[47] Fallon et al. characterized the free-energy landscapes of spike protein using MD simulation with enhanced sampling techniques and showed that RBD opening was modulated via interactions in an allosteric pocket.[42] Other enhanced sampling MD simulations have also been performed to study this conformational transition.[54] However, all of these MD simulations used fixed charge states for the ionizable residues on spike without rigorously considering the molecular response to environmental pH and thus were unable to elucidate the pH-dependent allosteric mechanism.

The constant pH MD (CpHMD) simulation allows sampling of the coupled protein conformations and corresponding protonation states of the titratable residues along the MD trajectory, which can provide an ensemble of conformations with properly weighted ionization states at a specific pH environment.[55,56] This is critically important as no single protonation state reasonably describes the ensemble when the pKa values of the titratable residues are near the environmental pH. In addition, due to the less polar and polarizable environment of the hydrophobic interior, the internal ionizable residues of protein often titrate with anomalous pKa values, which are very shifted from their standard values when they appear at the solvent-accessible protein surface.[57−60] Hence, the ionization states of the internal ionizable residues are often coupled with protein conformational changes, which enable protein pH-sensitive ability,[61−63] and many proteins harness these pH-dependent conformational changes for functional purposes.[64−67]

To provide more insights into RBD positioning in response to the environmental pH along the endosomal entry pathway, here we carry out microsecond-time scale CpHMD simulations on the full-length homotrimeric spike to quantify domain movements and to delineate the ionization of internal titratable residues that mediate RBD positioning. Previously, the RBD opening event has only been captured in the simulations using enhanced sampling strategy[42,43] or millisecond-time scale.[47] At present, we show that the ionization of the internal titratable residue would couple and hence prompt RBD opening transitions, which enable the observation of the pH-dependent RBD conformational changes in a relatively shorter simulation time scale. According to the recent studies,[47] glycosylation has little effect on the RBD-up/down equilibrium, and the difference between glycosylated and unglycosylated spikes was smaller than that between different spike variants. Therefore, the glycans on the spike surface were not considered in this study for computational efficiency. Our present work uncovers the pH-dependent allosteric mechanism of spike RBD that is essential for SARS-CoV-2’s infection process and hence supports the discovery of new therapeutic opportunities for COVID-19.

Methods

Spike Proteins in PDB

We first accessed the COVID-19/SARS-CoV-2 resources in the Protein Data Bank (PDB, https://www.rcsb.org/) and summarized all of the structures in the category termed “spike proteins and receptor-binding domains”, where there was a total of 775 related structures as of February 15, 2022. Then, the structures that have missing domains were deleted, and only the full-length spike protein structures were kept. Finally, we obtained 341 three-dimensional spike structures, as summarized in the Supporting Information, and all of them were homotrimers. Based on these structures, we plotted the experimentally observed RBD conformational distribution in terms of the internal angle and dihedral (defined in the latter section) measuring the RBD opening and rotation, as shown in Figure C.

Model Setup

The single-up and all-down three-dimensional structures of SARS-CoV-2 spike with PDB ids of 7dk3 and 7df3,[34] respectively, were used for the subsequent atomic simulations. The glycans and other small molecules were removed from the PDB files, and only the homotrimeric protein structure was kept. The missing residues were constructed using the Modeller program.[68] The internal D571 on spike was chosen as titratable. Two models were prepared according to the selection of titratable D571 on the homotrimer, namely, one model with only one titratable D571 (on protomer A), with the protonation states of the left two D571 residues (on protomers B and C) fixed at the protonated state and the other model with all three D571 residues on the homotrimer titratable. We defined protomer A as the one with the titratable D571 for the mono-D571-titratable model, and for the triple-D571-titratable model, the one with RBD opening to the up state first during the simulations would be defined as protomer A. Then, protomers A–C were defined as appearing counterclockwise when top-viewed along the homotrimer 3-fold axis. For clarity, we specify by subscripting the protomer of each subunit or residue. The Amber ff99SB force field was utilized to build the topological file for the spike protein via the tleap module in the Amber18 program.[69] The Generalized Born (GB) implicit solvent model[70] of igb = 8 with mbondi3 Born radii[71] was used to mimic the solvent environment. The all-RBD-down structure (PDB id: 7df3) was employed as the starting point for all simulations, and the RBD positioning during the simulations was quantified with reference to the experimental single-RBD-up structure (PDB id: 7dk3) to determine the occurrence of the RBD opening event.

Structure Equilibration

The system was first minimized with 10 kcal/(mol·Å2) positional restraints placed on the backbone atoms. Afterward, the system was heated gradually from 0.0 to 300.0 K over 0.5 ns with a 1 fs integration time step, again with 10.0 kcal/(mol·Å2) restraints placed on all of the heavy atoms. Then, the system was fully equilibrated for 1 ns MD simulation before the CpHMD production runs. No protonation state changes were attempted for the titratable residues during the equilibration procedure.

CpHMD Simulations

Starting from the equilibrated structure, each CpHMD simulation was run independently at a specific pH within the range of 4–7. The pH range has been validated to be sufficient in characterizing the whole ionization behavior of the investigated titratable residues. The CpHMD simulations were performed for both of the model homotrimeric spike systems (mono-D571-titratable and triple-D571-titratable). All production runs were carried out using the pmemd.cuda module of Amber18[69] with a time step of 2 fs. The GB implicit solvent model was used with a salt concentration of 0.1 M and an infinite cutoff for nonbonded interactions. The system temperature was maintained at 300 K using Langevin dynamics[72] with a collision frequency of 5.0 ps–1. The hydrogen-involved bonds were kept at their equilibrium lengths using the SHAKE algorithm.[73] Different random seeds were adopted to avoid synchronization artifacts.[74] The change of the protonation states for the titratable D571 was attempted every 100 fs. The protonation states of the other ionizable residues were assigned, as suggested from previous theoretical studies,[42] and were fixed throughout the simulations. The CpHMD simulation with mono-titratable D571 at each pH was run for at least 300 ns, and for simulation with triple-titratable D571, at least 400 ns production run was carried out to guarantee the converged result. The CpHMD trajectories were analyzed using the cpptraj module[75] of Amber18. The fractions of the protonation (or deprotonation) state, which can be evaluated as the time spent at the protonated (or deprotonated) state, were calculated using the cphstats program in AmberTools18.

To validate the reliability of the simulation protocol, we also performed CpHMD simulations on a model system Ace-Asp-Nme. The pKa of the titratable aspartic acid was estimated by fitting the fraction of deprotonation state fd at each pH to the Hill equation as followswhere n is the Hill coefficient, which should be close to 1.0 for noninteracting ionizable residues.[60]

Quantifying RBD Positioning

In the present work, we defined similar reaction coordinates using center-of-mass (COM) variables as that used in Fallon et al.’s study[42] to quantify the RBD motion during the CpHMD simulations. The definition of the reaction coordinates is shown in Table S1 and Figure S1 of the Supporting Information. An angle measuring the RBD opening up on the CTD1 pivot relative to the core central helices was defined. The three points used to define the RBD opening angle included the COM of the central β-strands in RBD, the COM of the β-strands in CTD1, and the COM of a short helical region in the upstream of the S2 domain. Moreover, a dihedral was also defined by adding the fourth point, which was the COM of the β-strand at the head of RBD, to describe the RBD rotation on its own long axis as with RBD opening. All of the four COM points were on the same protomer. The angle and dihedral were computed using the cpptraj module of Amber18.

Results and Discussion

To investigate the ionization of the internal ionizable residues on the spike RBD dynamics, we carried out two types of CpHMD simulations with mono- and triple-titratable D571 on the spike homotrimer, respectively. For each type, starting from the spike all-down structure, five independent CpHMD simulations at a specific pH (within a range of 4–7) were performed, which resulted in the sampling of structures in a total of >3.5 μs MD simulations. We first confirmed that the spike protein behaved reasonably in our simulation model. The backbone root-mean-square deviations (RMSDs) with reference to the initial down-state structure during MD simulations, as well as the root-mean-square fluctuations (RMSFs) of the backbone Cα atom, at different pH values are shown in Figure S2 of the Supporting Information. As can be seen from the figure, the low pH = 4 environment greatly stabilizes the spike all-down structure, with the RMSD of each protomer below 10.0 Å. The increased pH leads to larger RMSDs, indicating the conformational transitions from down to up configurations. As one RBD lifts away from the rest of the spike, the reduced packing of RBD leads to a more flexible spike structure, especially the open protomer with RMSD increasing to 12.0–13.0 Å. The RMSFs also show a similar pattern, with the RBD region of the opening protomer showing larger fluctuations at relatively higher pH conditions. These observations are consistent with the previous simulation results.[42,43] A molecular overlay, as shown in Figure S3 of the Supporting Information, demonstrates that the obtained RBD-up structure through our simulations agrees with the experimentally observed structure very well, with the RBD position and orientation perfectly matching the experiment, indicating that a successful down-to-up transition was captured from our CpHMD simulations. Hence, our simulation method is capable of describing the spike up/down structures, as well as the structural transitions, which can give a reasonable sampling of spike behaviors within the endosomal pH range.

Titration of the Internal D571 Residue

As observed from the experiment,[21] spike RBD opening is a pH-dependent manner, which indicates that the ionization of the internal ionizable residues in experiencing the environmental pH would have a remarkable effect on spike conformational transitions. The internal ionizable residues buried in the hydrophobic interior of the proteins often function as pH sensors to trigger the protein structural rearrangement,[48−52] and thus such internal ionizable residues on the spike were the primary focus of this study. The highly conserved D571 in the CTD1 domain is right under the cover of RBD (see Figure B), and it is completely surrounded by the CTD2, NTD, and S2 domains of the tightly packed homotrimeric spike. For the all-down state, the closed RBDs of the spike homotrimer tightly cover over D571 to shield it from the solvent forming a hydrophobic core. Therefore, the different protonation states of D571 were speculated to make a difference in the RBD dynamics and were validated in our present CpHMD simulations. Protein conformational changes coupled to the ionization of the internal ionizable residue can be reflected by its pKa value. The titration curves fitted based on the deprotonation fraction of D571 over the CpHMD simulations, in comparison with the titration result for a capped aspartic model system (Ace-Asp-NME), are shown in Figure A. As one can see from Figure , the titration curves for the mono- and triple-titratable D571 were notably shifted from that of the aspartic model system, indicating that the internal D571 experienced an electrostatic environment quite different from the solvent-accessible aspartic acid. The computed pKa of the model aspartic acid from our CpHMD simulation was 4.02, in excellent agreement with its normal pKa value of 4.0 in solution.[76] However, for the internal D571 on the spike homotrimer, both of the computed pKa values of 5.50 and 5.87 for mono- and triple-titratable D571, respectively, considerably deviated from its reference value. The shifted pKa values of D571 lie right within the endosomal pH range, indicating no single protonation state of D571 along the endosomal entry pathway.

Figure 2

(A) Titration curves for the mono-titratable (red), triple-titratable (blue) D571 on spike homotrimer, and the capped aspartic model system (black). (B) Conformation-specific titration curves for mono-titratable D571 buried (light blue), accessible to the solvent (violet), and for the complete ensemble (red).

We also computed the individual titration curves only for those structures unambiguously assigned as up or down states to obtain the conformation-specific pKa values. We select the spike up or down conformations based on the solvent-accessible surface area (SASA) of D571 (see the Supporting Information). The resulting titration curves are shown in Figure B. The computed pKa of the deeply buried D571 in the all-down conformations is further up-shifted to 6.10, while the pKa for the solvent-accessible D571 in the single-up conformations decreases to 4.87, which is closer to the normal pKa of a solvent-exposed aspartic acid. These calculations confirm that D571 has conformation-specific pKa values, as required to regulate pH-dependent allostery.

These computed pKa values indicate that D571 tended to be protonated in the hydrophobic core of the spike homotrimer under conditions with pH < 5.5 and thus help to stabilize the all-down state. The increase in pH (>5.5) eventually caused D571 to be deprotonated resulting in a negative charge in the hydrophobic core, which would be energetically unfavorable, thus promoting the destabilization of RBD and subsequent RBD opening. This supports the experimental observations that single-up structures could only occur at pH > 5.5.[21] Therefore, it could be expected that the ionization of D571 could provide a pH-regulated mechanism for RBD positioning.

RBD Opening

To engage with ACE2, spike RBD must undergo conformational transitions from down to up states to expose the binding interface first. Characterizing the full range of spike opening and related activation mechanisms are of vital importance for understanding pathogenesis and could provide insights into novel therapeutic options. To get a first glance at the behaviors of RBD, we projected the experimentally determined spike structures from PDB onto 2 order parameters, the opening angle, and the rotation of RBD. As shown in Figure C, these experimental structures are clearly separated into two groups, representing the RBD-down state with the angle mainly within 52–65° and rotation around 60–110°. The reduced packing of RBD, as it transitions from the down to up state, leads to a larger opening (70–111°) and rotation (35–110°) range, indicating looser constraint for RBD in the up state.

According to the previous studies,[42,43,47] capturing spike opening events was a challenging computational task, which usually required large-scale simulations and/or enhanced sampling techniques. Previous microsecond-time scale simulations using straight MD sampling approaches were unable to capture the RBD opening process.[53,77] Enhanced sampling techniques have significantly promoted the capturing of spike RBD opening event[42,43,47] but could not reveal the pH-dependent allosteric mechanism of the spike protein. In present CpHMD simulations, we find that the ionization of D571 in response to the environmental pH can significantly facilitate the spike RBD-up/down transitions, leading to the occurrence of RBD opening in hundreds of nanoseconds to microsecond time scale.

The transitions of RBD during the CpHMD simulations with mono-titratable D571 at different pH values are quantified by the RBD opening angle, as shown in Figure . At pH = 4, in which the internal D571 was predominantly protonated, the spike protein was very stable in the all-down state with the RBD angle of each protomer around 60°, which is in good agreement with the experimental observation (see Figure C). The deprotonation of D571 by pH increase activates the spike RBD (protomer A with titratable D571) transition from the down to up state. For simulations at pH ≥ 5, the RBD angles of the nontitratable protomers B and C were also very stable around 50–60°, whereas for the protomer A with titratable D571, the RBD angle significantly increases and finally becomes steady at around 90–100°, indicating a single-up structure. These simulation results are in excellent agreement with the experimental observations (see Figure C). As can be seen from Figure , all of the spike conformational transitions occur within 200 ns CpHMD simulations for pH values ≥5, and the full down-to-up transition process, corresponding to the region with a sharp angle increase, can be completed within ∼50 ns. Hence, the deprotonation of D571 in the hydrophobic interior of the all-down structure at serological and early-endosome pH values (6–7) could significantly prompt the RBD opening, facilitating the course of infection by recognizing and binding of ACE2. Whereas at the low late-endosome pH, the spike would be more probably stabilized in the all-down state with protonated D571, providing a potential means of immune evasion from RBD-up-recognizing antibodies. This view is in line with the experimental observation.[21] Our simulation results support the hypothesis that D571 functions as a pH sensor that guides the up/down conformation of RBD upon the pH change. As such, the pH-induced conformational switch has also been observed in many other viral fusion proteins,[36,38,39] the revealed mechanism in the present work would provide a general physicochemical insight for understanding the infection courses of other fatal viruses.

Figure 3

Evolution of the spike RBD opening angle of each protomer over the CpHMD simulations with mono-titratable D571 under different pH conditions.

Moreover, the increased pH could motivate the spike RBD to undergo a substantially larger opening scale, with the opening angle beyond 110° and within 140° (see pH = 7 in Figure ), which is larger than what has been observed in the experiment and indicates a fully open structure. The comparison of the simulated RBD-down, up, and open conformations is shown in Figure . This simulated large scale of RBD opening is consistent with the observations from previous simulation studies.[43,47] The dramatic opening indicates the flexibility of RBD and suggests that antibodies, as well as other therapeutics, can potentially bind to the regions of RBD that are deeply buried and seemingly inaccessible in existing experimental structures.[47]

Figure 4

Structures of a spike protomer sampled from the CpHMD simulations and colored in light blue, violet, and red, denoting the RBD-down, up, and open states, respectively.

For the CpHMD simulations with triple-titratable D571, a similar pattern is also observed for the RBD dynamics at different pH values. As the predicted pKa (5.87) for the triple-titratable D571 is relatively larger than that (5.5) for the mono-titratable one (see Figure A), the RBD opening (angle >70°) is observed at pH ≥ 5.5 for the simulations with triple-titratable D571, while at pH < 5.5, the spike protein is stabilized in the all-down state with the RBD angle of each protomer below 70° during the 400 ns MD simulations (see Figure S4 of the Supporting Information). At pH = 5.5, the RBD-down-to-up transitions (angle >70°) occur on both protomers A and B, indicating a double-RBD-up structure, but for pH 6 and 7, the RBD transitions are only observed on the protomer A, indicating single-RBD-up structures. These simulation results suggest that the spike single-up structure is always dominant no matter with mono- or triple-ionized D571, which is consistent with the experimental observations.[21,32]

Free-Energy Landscape of the pH-Dependent Conformational Changes

To provide more information about the RBD positions and their relative free energies, we calculated free-energy landscapes in the form of potentials of mean force (PMF)[78] for the RBD dynamics. The predefined RBD opening and rotation angles were used as the reaction coordinates to describe the RBD conformations. The PMFs based on the simulations with mono-titratable D571 at the low (pH = 4) and high (pH = 7) pH values are illustrated in Figure . At pH = 4, the close minima correspond to the all-down conformation of the spike homotrimer, with the RBD angle of each low-energy protomer within 50–65°. However, these minima are not located at the same position, indicating the asymmetric and nonidentical RBD structures of the spike homotrimer. The non-RBD transition at pH = 4 also indicates a significantly higher down-to-up energy barrier (>3.0 kcal/mol) for the spike with protonated D571. At pH = 7, the energy-minimum conformation of protomer A is shifted to a larger angle region (90–100°) representing the RBD single-up structure. Meanwhile, the opening RBD twists inward, as quantified by the decrease in the dihedral measuring RBD rotation. This is consistent with Zimmerman et al.’s simulations and further supports the conformational selection mechanism where the spike first opens and twists to an appropriate position before binding to its partners.[47] As shown in Figure , at pH = 7, the RBD-down-to-up transition requires an activation free energy of ∼3.0 kcal/mol, which is lower than the energy barrier required at pH = 4, and the free energy of the up state is ca. 0.6 kcal/mol higher than the down state. According to Fallon et al.’s simulations with fixed protonation states,[42] the predicted activation energy barrier for the RBD-down-to-up transition is ∼6.0 kcal/mol, which is significantly larger than the value obtained in this study. Therefore, the ionization of internal D571 could significantly decrease the transition barrier and thus make it much easier to capture the spike’s up/down conformational transitions. The PMFs obtained at other pH values also show similar patterns, as illustrated in Figure S5 of the Supporting Information.

Figure 5

Free-energy landscapes for RBD positioning at pH values of 4 and 7. The energy minima are labeled with the corresponding protomers.

Atomic Details of the Opening Mechanism

The ionization of internal D571 could trigger spike RBD conformational changes to tolerate the corresponding ionization, and RBD conformational changes, in turn, facilitate the ionization of internal D571. The conformational transitions coupled with the ionization of internal D571 provide an explanation for the pH-dependent allostery of the spike RBD. To identify critical components responsible for the different spike active states, the detailed atomic mechanisms of spike conformational transitions are also characterized based on MD simulations with mono-titratable D571 in this work. In the RBD-down state, the carboxylic groups of D427, D428, and E516 on RBD are right over the internal D571 with a distance of ∼17.0 Å, as shown in Figure A. Since these three ionizable residues (D427, D428, and E516) locate more closely to the protein–water interface and are more readily accessible to the solvent, their deprotonated states should be more preferred in comparison with the internal D571. Therefore, the repulsion between the carboxylic groups on RBD (D427, D428, and E516) and D571 due to the deprotonation of D571 would be non-negligible in the low-dielectric hydrophobic core, which would lead to the instability of RBD and thus opening transition. The RBD transition significantly increases the distance between the carboxylic groups on RBD (D427, D428, and E516) and D571 to ∼38.0 Å. Due to the RBD opening of one protomer with increasing the solvent accessibility of the former hydrophobic core, the environmental polarization, as well as the increased dielectric constant, weakens the repulsion between the carboxylic groups on the remaining two protomers, allowing the rest RBDs still stable in the down states even with deprotonated D571. This could potentially explain why the single-RBD-up conformation is dominant, with fewer double- and triple-RBD-up conformations, at the serological and early-endosome pH values, as observed from the experiment.

Figure 6

(A) The aspartic acids over D571 in the RBD-down conformation. (B) Interprotomer contacts between residues on the RBD of protomer A and S2 domain of the neighboring protomer B stabilizing the RBD-down state. The dashed lines represent the hydrogen-bonding interactions.

The RBD-down state features four interprotomer hydrogen bonds and an interprotomer salt bridge formed at the RBD-S2 (of the neighboring protomer) interface, as shown in Figure B. The S383A on the RBDA of the protomer A simultaneously forms three hydrogen bonds with R983B, D985B, and E988B on the S2B domain of the neighboring protomer B. The fourth hydrogen bond is formed between T385A on the RBDA and D985B on the S2B domain. Moreover, K386A on the RBDA can form a salt bridge with D985B on the S2B domain. These interactions play an important role in stabilizing the RBD-down state. The evolution of lengths of the four hydrogen bonds and the distance of the salt bridge over the simulations are shown in Figure S6 of the Supporting Information, indicating a concordance between these critical interactions (breaking coordinated interprotomer interactions) and the RBD transition (see Figure ). As shown in Figure S6, the interprotomer interactions between A and B are broken at pH ≥ 5, prompting the RBDA to twist upward, while the B–C and C–A interprotomer interactions are still very stable during all of the simulations, indicating a single-up transition. As a result, the deprotonation of D571 would interrupt these interprotomer interactions to prompt RBD opening.

The RBD opening is proceeded through a hinge region connecting the RBD and CTD1 domains. The RBD–CTD1 interface serves as a hinge around which the RBD can rotate up and outward. The hinge region undergoes local conformational changes as the RBD opens, and the new forming interdomain interactions (hydrogen bonds and salt bridges) at the RBD–CTD1 interface would maintain the RBD-up state. We find that the up RBD has a variety of styles and is more flexible, which is consistent with the experimental observation for the wider population of the RBD opening and rotation angles of the RBD-up structures (see Figure C). Sztain et al. also found quite a lot of independent RBD opening transition pathways, indicating the diverse poses for the RBD-up structures.[43] The flexible RBD is beneficial to virus evading immune recognition. The up RBD structures obtained at different pH environments are usually distinct, and accordingly, conformation-specific interactions may be required to stabilize the multistyle up RBDs (see Figure ). Overall, the N360A and A522A on the RBDA hinge region play an important part in stabilizing the various RBD-up poses. For the up structure obtained at pH = 5, the N360A and A522A residues on RBDA form two hydrogen bonds with P579A and N544A on CTD1A, respectively, while these two residues on RBDA could form hydrogen bonds with D574A and F543A on CTD1A to stabilize the up structure at pH = 6. Three hydrogen bonds (namely, S359A, N360A, and A522A on RBDA with Q580A, T581A, and A575A on CTD1A, respectively) formed at the RBDA–CTD1A interface could maintain the up structure at pH = 5.8, and two hydrogen bonds (namely, L335A and A522A on RBDA with T581A and F543A on CTD1A, respectively) and a salt bridge (K528A on RBDA with E583A on CTD1A) help to stabilize the up structure at pH = 7. The evolution of lengths of these hydrogen bonds and the distance of the salt bridge over the simulations are shown in Figure S7 of the Supporting Information. As shown in the figures, the spike RBD-down-to-up transitions, as well as the maintaining of the stability of the up state, are accompanied by the formation of these critical interactions. Our simulations provide atomic details on the role of these critical contacts in modulating RBD opening.

Figure 7

Conformation-specific interdomain interactions in the hinge region maintaining the multistyle RBD-up conformations obtained from pH values of 5 (A), 5.8 (B), 6 (C), and 7 (D). The dashed lines represent the hydrogen-bonding interactions.

Conclusions

As a prime target for neutralizing antibodies and drug design, SARS-CoV-2 spike functions along the endosomal entry pathway by utilizing a pH-dependent allosteric strategy, but less is known about the spike pH-sensitive elements gating the transitions between the observed states of the RBD. Here, we report extensive CpHMD simulations of the full-length spike homotrimer characterizing the transition from the down to up states of the RBD under different pH conditions. The present simulations successfully capture the RBD transition events and provide multiple independent RBD opening transition pathways. The simulated up structures align very well with the experimentally determined structure. Our simulations indicate a critical gating role of the highly conserved D571, which is protonated in the RBD-covered hydrophobic core of the spike homotrimer at low pH and can act as a pH sensor with a clear pKa shift. Deprotonation by a pH increase would result in a negative charge within the hydrophobic core, and thus the energetically unfavorable situation prompts RBD transition. The ionization of D571 coupling with SARS-CoV-2 spike up/down equilibrium provides a rational explanation for the experimental observation of the spike RBD’s pH-dependent switch.

We computed the free-energy landscapes of the pH-dependent conformational changes to quantify the distribution of the RBD-up/down states under different pH conditions. We find that the RBD-down state (with protonated D571) is always energetically more stable, with ∼0.6 kcal/mol free energy lower than the up state (with deprotonated D571), and the RBD-down-to-up transition with deprotonated D571 requires a lower activation energy (∼3.0 kcal/mol) than that with protonated D571. This energy barrier is also significantly lower than previously reported (∼6.0 kcal/mol) from MD simulations without considering the ionization of D571 under the pH conditions. Therefore, the ionization of internal D571 could significantly decrease the activation barrier and hence make it much easier for the spike to undergo down-to-up conformational transitions.

To probe more details of the RBD transitions, we also characterized the intramolecular interactions that modulate the RBD-up/down states. The RBD-down state is stabilized by the interprotomer interactions (hydrogen bonds and salt bridges) formed between the RBD and S2 domain of the neighboring protomer, and the RBD transition is accompanied by these broken interactions. The RBD opening is proceeded through a hinge region linking the RBD and CTD1 domains. The up RBD becomes more flexible due to the reduced structural packing of the spike homotrimer, and distinct up structures are sampled under different pH conditions. Therefore, conformation-specific interactions are required to stabilize the multistyle up RBDs. This work provides new insights into the mechanisms of viral infection and hence supports the discovery of novel therapeutics for COVID-19.