Uncovering pH-dependent transient states of proteins with buried ionizable residues.

The role of pH in regulating biological activity is ubiquitous, and understanding pH-mediated activity has traditionally relied on analyzing static biomolecular structures of highly populated ground states solved near physiological pH. However, recent advances have shown the increasing importance of transiently populated states, the characterization of which is extremely challenging but made plausible with the development of techniques such as relaxation dispersion NMR spectroscopy. To unlock the pH dependence of these transient states with atomistic-level details, we applied the recently developed explicit solvent constant pH molecular dynamics (CPHMD(MSλD)) framework to a series of staphylococcal nuclease (SNase) mutants with buried ionizable residues and probed their dynamics in different pH environments. Among our key findings is the existence of open states in all SNase mutants containing "buried" residues with highly shifted pKa's, where local solvation around the protonation site was observed. The calculated pKa demonstrated good agreement with experimental pKa's, with a low average unsigned error of 1.3 pKa units and correlation coefficient R(2) = 0.78. Sampling both open and closed states in their respective pH range, where they are expected to be dominant, was necessary to reproduce experimental pKa's, and in the most extreme examples of pKa shifts measured, it can be interpreted that the open-state structures are transient at physiological pH, contributing a small population of 1-2%. This suggests that buried ionizable residues can trigger conformational fluctuations that may be observed as transient-state structures at physiological pH. Furthermore, the coupled relationship of both open and closed states and their role in recapitulating macroscopic experimental observables suggest that structural analysis of buried residues may benefit from looking at structural pairs, as opposed to the conventional approach of looking at a single static ground-state conformation.

One of the critical regulators of biological activity is pH, and misregulated pH is a hallmark of cancer-related physiology.[5] Ionizable groups that change their protonation states as a function of pH and their local electrostatic environment are often the key residues that govern pH-dependent activity. Due to their charged nature, the majority of them are expressed near the protein surface; however, in some systems, such as membrane proteins, they can be found in more hydrophobic environments. Until recently, the dominant approach to probe the mechanism of such pH-dependent activity has been focused on the analysis of static structures typically solved at physiological pH. However, recent studies have demonstrated the increasingly important role of biomolecular excited-state structures, typically referring to populations that are 0.1% to a few percent of the total population under physiological conditions[6] (hereon referred to as “transient states” to avoid confusion with electronically excited states), in protein folding[7] and ligand binding.[8] While the list of examples of pH-dependent transient states is not that long, there is precedence for their importance, such as in membrane fusion involving the influenza hemagglutinin HA2 subunit.[10]

The characterization of transient states is beyond the detection limits of conventional experimental techniques, but recent developments in relaxation dispersion NMR spectroscopy[11] and room-temperature X-ray crystallography[12] have made this endeavor feasible. Probing pH-dependent transient states will undoubtedly be similarly challenging. In this context, we report on the application of novel computational techniques to model pH-dependent dynamics of proteins and to probe their transiently populated minor states, which can serve as complementary tools in the field since simulations are not subject to the same detection limits as experiments. One major approach to model accurate pH-dependent dynamics is the constant-pH molecular dynamics (CPHMD) framework, where MD simulation is coupled to the protonation state of the titrating residue.[13] In the continuous variant of CPHMD used in this study, the titration coordinate represents an instantaneous microstate, and it is propagated continuously between the protonated and unprotonated states using the λ dynamics approach.[16] We focus our investigation on modeling the pH-dependent dynamics of a series of engineered mutants of staphylococcal nuclease (SNase) with Lys residues buried in the hydrophobic core[18] using the recently developed explicit solvent CPHMDMSλD framework for proteins.[17] While the effectiveness of CPHMD has been demonstrated on numerous systems,[19] almost all work reported to date was based on an implicit solvent model.[13] Applications of explicit solvent CPHMD on biomolecules thus far have only a few reported successes.[17,24] More importantly, none of the existing works have attempted a comprehensive investigation of buried ionizable residues, where, we hypothesize, the contrast in electrostatic environment between a hydrophobic pocket and a solvent-exposed environment will provide a key driving force for pH-mediated conformational fluctuations and, possibly, transient states formation.

As shown in Figure 1a, the mutants for this study comprise a set of diverse mutation sites, with residues varying in the magnitudes of pKa shifts. In addition, experiments have demonstrated that, for the mutants that we selected, the titration of the internal Lys is decoupled from the ionization of Asp and Glu residues.[18] Thus, to facilitate convergence within a reasonable time, we titrated only the buried Lys and performed an initial check of the accuracy of our simulations by calculating its pKa value, using the crystallographic structures as the starting models (see Supporting Information (SI) for details). As summarized in Figure 1b and Table S1, most of the buried residues have a predicted pKa < 1, which is shifted by nearly 10 pKa units from the reference value of 10.4. pKa calculations on the V66D and V66E mutants produced a predicted pKa > 14, which is dramatically shifted from their reference values of 4.0 and 4.4, respectively. In both cases, the pKa shifted toward the direction stabilizing the neutral state of each titrating species (upward for Asp/Glu; downward for Lys). The dramatic pKa shifts observed in our simulations mirror the experimental spectrophotometric measurements of Brønsted acids and bases, where a similar pKa shift of 10–20 units was recorded when moving from an aqueous to an organic environment.[28] For example, acetic acid, which has the same functional group as Asp and Glu, has its pKa value shifted from 4.8 to ∼23 in water (ε = 80) vs acetonitrile (ε = 37).[29] Even though the dielectric constant of a typical hydrophobic pocket in SNase, which ranges from 4 to 20,[30] is lower than that of acetonitrile, such dramatic pKa shifts of more than 10 units have never been observed in the protein. In fact, some of the largest pKa shifts for Asp, Glu, and Lys reported in the literature are perturbed by a mere 5 units from their reference values.[18,33] Thus, there is an inconsistency in the magnitude of experimental pKa shifts observed in SNase compared to historical experimental pKa shifts observed in organic solvents, even though both environments have similarly low dielectric values, which indicates that additional factors have to account for the apparent inconsistency.

Figure 1

(a) Distribution of internal Lys residues of SNase mutants, color-coded depending on the pKa shift: yellow, not shifted; orange, shifted by 1–2 units; red, shifted by >2 units. Comparison between experimental and calculated pKa’s from explicit solvent pH-REX CPHMDMSλD simulations using (b) only closed crystallographic like structures and (c) both open and closed states.

Strong coupling between pKa and conformational sampling has been previously reported,[34] and, while the importance of conformational sampling has been observed in some SNase mutants, no connection to experimental pKa values has been reported.[38] Based on these observations, we performed an additional series of extended-run explicit solvent pH-REX CPHMDMSλD simulations on the structures that were pre-equilibrated for 50 ns at high and low pH (see SI for details). The findings, summarized in Figure 1c and Table S1, demonstrate significantly improved results, with the averaged unsigned error (AUE) reduced from 3.9 to 1.2 pKa units. We observed that the longer equilibration allowed the sampling of “open” solvated structures that were critical for reproducing the experimental pKa values reported for all Lys mutants with highly shifted pKa values. This is in contrast to the “closed” crystallographic-like structure where there is little to no water present within 3 Å of the protonation site (Figure 2a). Using the V66K mutant as an example, we demonstrate that SNase conformations sampled in our simulations interconvert between closed and open states (Figure 2b,c) when the external pH is close to the calculated pKa value. While the coupling between pKa and conformational sampling has been reported for a few isolated examples,[36,37] and the importance of sampling such alternative conformational states has been previously postulated by experiments,[39] there has been no comprehensive proof for their role. Our work presents the first interpretation that pH-dependent transient states exist and may be of general importance for proteins with buried ionizable groups. For all mutants with highly shifted pKa values, we observed that protonation of the internal Lys was concomitant with an increase in local solvation around the protonation site, as illustrated in Figure 2a. The backbone RMSD of the entire structure between both closed and open states in these mutants is small, ranging from 1.2 to 1.9 Å (also see Figure S1), which suggests that the conformational relaxation to accommodate a buried charged residue does not require significant structural rearrangement. Our observations are consistent with experimental measurements indicating that buried ionizable residues in SNase are readily accommodated without any special structural adaptation or distortion to the overall protein conformation.[18,33]

Figure 2

(a) Structures of the four most highly shifted Lys mutants at high-pH “closed” conformation (blue) and low-pH “open” conformation (red), with the backbone RMSD between the two structures in parentheses. Water molecules within 3 Å of the protonation site at low pH (green) are included; no water molecules were observed at high pH. Using V66K, we show the conversion between the two states at external pH close to the calculated pKa in the time evolution of (b) number of waters within 3 Å of the protonating site and (c) backbone RSMD relative to the crystallographic structure.

Having established that buried ionizable residues can trigger pH-dependent structural fluctuations, we extend our analysis to determine the relevance of these alternative states at pH 7. Using a two-state model that assumes a conversion between one dominant open state and one dominant closed state, we derive an equation (see SI for derivation) describing the ratio of open to closed states, ROC:

This function can be decomposed into pH-dependent (KpH) and pH-independent (K0) terms. The KpH term depends on the microscopic pKa of each state (pKclosed, pKopen) and external pH. The K0 term can be physically related to the free energy difference of the open and closed states in their protonated form (SI). However, as each form favors opposing protonation states, this free energy is usually not measured by experiments. Thus, it may be advantageous to express K0 as a function of the system’s readily measured macroscopic or apparent pKa (pKapp) and the microscopic pKa of each state (pKclosed, pKopen). Based on the ROC equation, one may also derive the fraction of the open state (Fopen), a more intuitive metric for open states at a specific pH:

Due to the rapid dynamics and/or low population of the minor state, it is often beyond the detection limits of experiments to establish the pKa values of each microscopic state.[39] However, because the open state is solvent-exposed, pKopen may be approximated as the reference pKa of the free amino acid. Using the range of pKa shifts of 10–20 units recorded from spectrophotometric measurements of organic acids and bases in a low dielectric solvent,[28] we have conservatively assigned pKclosed to be shifted by 10. As shown in Figure 3, one can use the ROC equation to derive the pH-dependent fraction of the open state for a series of hypothetical pKapp for buried Lys or Glu.

Figure 3

pH-dependent distribution of the fraction of open states (Fopen) for buried (a) Lys and (b) Glu residues. Several color-coded hypothetical apparent pKa (pKapp) values are illustrated, shifted by 1−5 pKa units. In the limits of extreme (i.e., 10) or null pKa shift, the expected population of 100% closed or open states is recovered. For this plot, pKopen = 10.4 and 4.4, and pKclosed = 0.4 and 14.4, for Lys and Glu, respectively.

Our analysis indicates that pH-dependent transient open states may contribute as much as 2% of the total population at pH 7 when the apparent pKa of Lys is shifted by as much as 5 units, which appears to be the upper limit of pKa shift as recorded in current literature.[18,33] For Asp and Glu, a pKapp shifted by 5 pKa units represents a 1% contribution of the transient state at pH 7. Thus, for the residues with highly shifted pKa values, the low-population transient states are likely to contribute significantly to the apparent pKa and, thus, need to be elucidated to correctly compute the apparent pKa. Although the existence of transient states involving buried ionizable groups does not necessarily imply a functional relevance, there is increasing precedence that the inclusion of transient states is needed to fully account for biological properties.[6−8] In the context of our work, we suggest that the effect of such pH-dependent transient states will be pronounced when an ionizable group transitions between hydrophilic and hydrophobic environments, such as in membrane fusion processes, where activated/transient states have been postulated to play a crucial role.[10] In addition, traditional studies of catalytic mechanisms have always assumed that crystallographic structures correlate to measured pKa, but, as we have shown, that may not be true for buried residues with highly shifted pKa values. Moreover, the coupled relationship between open and closed states and their role in recapitulating macroscopic experimental observables suggests that structural analysis of buried residues should be performed from the perspective of looking at structural pairs, as opposed to the conventional approach of looking at a single static ground-state conformation. For such analyses, the equations we have provided will prove useful for a quick, “back of the envelope” estimation of the population of proposed transient states. For example, one could use the experimentally measured pKapp to select the appropriate curve in Figure 3, and use it to estimate the fraction of the open state at a given pH, as a means to evaluate the plausibility of experimental characterization within the detection limits of the methods employed. Alternatively, it could also be used to estimate the pH range where the population of proposed transient states will enter the detection range of experiments.

Lastly, we investigate the various computational models, specifically previous CPHMD implementations, to ascertain the robustness of predictions obtained from simulations. One insightful observation is that our calculated pKa for the N100K mutant of SNase is 6.6, which is close to the calculated pKa of 7.0 reported by Shen and co-workers,[37] despite the differences in the solvation model and cutoff schemes utilized. This suggests that the CPHMD framework is not overly sensitive to specifics of the simulation setup. This is because CPHMD calculates the free energy of protonation relative to a reference compound, and, as long as the simulations of the protein and the reference are performed under identical conditions, the differences originating from the simulation setup approximately cancel out. To test this hypothesis, we compare our pKa predictions from explicit solvent pH-REX CPHMDMSλD simulations with those obtained from the GBSW[14] implicit solvent pH-REX CPHMD framework, using the refined protocol presented in this paper. Our results (Figure 4a,b and Table S2) show excellent correlation of R2 = 0.89 between the pKa’s predicted from both explicit and implicit models. However, unlike the explicit solvent CPHMDMSλD simulations, which resulted in a pKa shift of more than 10 units when only the closed state was used, the implicit solvent CPHMD simulations produced pKa shifts that were smaller in magnitude. The van der Waals surface representation used to define the solute–solvent dielectric boundary in GBSW is known to form small crevices of high dielectric region between atoms, which results in the underestimation of the Born radii and overestimation of the solvation free energy.[40] Therefore, despite the actual hydrophobicity of the environment near the protonating site, the GBSW model may be too “wet”, causing a smaller pKa shift. Alternatively, the same effect can be achieved as a result of the faster conformational dynamics in the GBSW model, which may sample both open and closed states more frequently.

Figure 4

(a) Accuracy of calculated pKa from implicit solvent pH-REX CPHMD simulations is similar to that of explicit solvent results, and (b) both sets of calculated pKa values are highly correlated with each other. (c) Calculated pKa values from explicit solvent pH-REX CPHMDMSλD simulations using the CHARMM36 force field result in improved correlation with experimental pKa values. All simulations were initiated using both closed and open structures.

To distinguish between these two possibilities, additional simulations, with the structures of the V66K, V99K, L125K, and I92K mutants rigidified by applying harmonic restraints to all heavy atoms, were performed, and the resulting pKa values were similar to those calculated without restraints (Table S2). These results suggest that smaller pKa shift in the implicit solvent primarily stems from the “wetness” of the GBSW model, rather than its faster conformational dynamics. Finally, we investigated the effect of using the recently released CHARMM36 all-atom force field for proteins on our predictions, as it has been previously shown to yield superior reproduction of experimental dynamic data.[41] We recalculated the biasing potentials for the common titrating residues of proteins using the CHARMM36 force field (Table S4), and revised pKa values, as shown Figure 4c, indicate an improvement in the predictive performance over the older CHARMM22/CMAP force field, with R2 increasing from 0.52 to 0.78, and the slope of the regression moving from 0.52 to 1.09, while maintaining the same level of accuracy with an average unsigned error of 1.3 pKa units.

In conclusion, we have applied the recently developed explicit solvent constant pH molecular dynamics framework to simulate the pH-dependent dynamics of a comprehensive set of SNase mutants with buried ionizable residues that have varying degrees of pKa shifts. Among our key findings are that a buried charged residue cannot be accommodated inside a purely hydrophobic pocket, and that an open-state structure for these “buried” residues, characterized by local solvation around the protonating site, was observed in all SNase mutants with highly shifted pKa’s. At physiological pH, buried ionizable groups with large pKa shifts have transiently populated open states, where they contribute a small but non-zero population of 1–2%. Nevertheless, sampling these open states is a necessary condition for accurately reproducing experimental pKa measurements, with which calculated pKa from our explicit solvent CPHMDMSλD simulations demonstrated good agreement, with a low average unsigned error of 1.3 pKa units and correlation coefficient R2 = 0.78. This work provides the first validation that buried ionizable residues can readily trigger pH-mediated conformational fluctuations that may be observed as transient-state structures at physiological pH. Lastly, the discovery of a coupled relationship of both open and closed states and their role in recapitulating macroscopic experimental observables suggests that structural analysis of buried residues may benefit from looking at structural pairs, as opposed to the conventional approach of looking at a single static ground-state conformation.