pH-Dependent Conformational Changes Due to Ionizable Residues in

Article Cite This: J. Phys. Chem. B 2019, 123, 5742−5754

pubs.acs.org/JPCB

pH-Dependent Conformational Changes Due to Ionizable Residues in a Hydrophobic Protein Interior: The Study of L25K and L125K Variants of SNase Ankita Sarkar,† Pancham Lal Gupta,‡ and Adrian E. Roitberg*,‡ †

Department of Physics, University of Florida, Gainesville, Florida 32611, United States Department of Chemistry, University of Florida, Gainesville, Florida 32603, United States

‡

Downloaded via NOTTINGHAM TRENT UNIV on August 13, 2019 at 15:12:10 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information *

ABSTRACT: Ionizable residues in the hydrophobic interior of certain proteins are known to play important roles in life processes like energy transduction and enzyme catalysis. These internal ionizable residues show experimental apparent pKa values having large shifts as compared to their values in solution. In the present work, we study the pH-dependent conformational changes undergone by two variants of staphylococcal nuclease (SNase), L25K and L125K, using pH replica exchange molecular dynamics (pH-REMD) in explicit solvent. Our results show that the observed pKa of Lys25 and Lys125 are significantly different than their pKa in solution. We observed that the internal lysine residues prefer to be water-exposed when protonated at low pH, but they remain buried within the hydrophobic pocket when deprotonated at high pH. Using thermodynamic laws, we estimate the microscopic conformation-specific pKa of the water-exposed and buried conformations of the internal lysine residues and explain their relation to the macroscopic observed pKa values. We present the differences in the microscopic mechanisms that lead to similar experimentally observed apparent pKa of Lys25 and Lys125, and explain the need of thermodynamic models of different complexities to account for our calculations. We see that L25K displays pH-dependent fluctuations throughout the entire β barrel and the α1 helix. In contrast, pH-independent fluctuations are observed in L125K, primarily limited to the α3 helix. The present computational study offers a detailed atomistic understanding of the determinants of the observed anomalous pKa of internal ionizable residues, bolstering the experimental findings.

■

INTRODUCTION Charged residues located within the hydrophobic core of a protein can destabilize its folded structure. Consequently, ionizable residues are found in abundance at the protein− water interface, where they interact with the aqueous media.1 While rare, some proteins have ionizable residues buried in their hydrophobic core. These internal ionizable residues play major roles in enzyme catalysis,2 electron transport in systems like ATPase3 and cytochrome c oxidase,4 as well as energy transduction in complex membrane proteins, such as bacteriorhodopsin.5 Interestingly, the experimentally observed apparent pKa of these buried ionizable amino acid residues are highly deviated with respect to their intrinsic pKa in bulk water.6 These anomalous pKa values of the internal ionizable residues have been found to be essential for biological functions associated with the proteins.7 Hence, the atomistic mechanisms driving the adaptations undergone by these proteins to tolerate the presence of buried ionizable residues are of special interest and require extensive study. Protein function is modulated by its structure and dynamics, which can be pH-dependent.8−13 The pH-dependence arises © 2019 American Chemical Society

due to the presence of amino acid residues with side chains capable of altering their protonation states with pH changes. Therefore, pH changes can lead to conformational reorganizations, modifying the pKa values of the ionizable residues. In proteins with internal ionizable residues, the pH-dependent conformational changes accompanied by alteration of protonation states of the internal ionizable residues are reflected in their observed apparent pKa values.11,13−15 Electrostatic interaction6,16,17 can affect the pKa of ionizable residues. An internal ionizable residue can undergo interactions with fixed dipoles in the protein backbone18 and with other ionizable residues, located both at the surface and the interior of the protein. Thermal fluctuations of charges and dipoles around their equilibrium positions19,20 and reorganization of permanent dipoles in the presence of charged residues21 can modify the interactions of key internal residues with its surroundings. Fluctuations of protein structures can lead to Received: April 23, 2019 Revised: June 17, 2019 Published: June 17, 2019 5742

DOI: 10.1021/acs.jpcb.9b03816 J. Phys. Chem. B 2019, 123, 5742−5754

Article

The Journal of Physical Chemistry B

undergo fluctuations that are relatively pH-independent and are limited to hundreds of microseconds, affecting residues only in the vicinity of Lys125. Previous computational studies with SNase variants have indicated the correlation between the magnitude of shifts in pKa of internal ionizable residues and different parameters of structural modifications, particularly the extent of water exposure of the key internal residues in their charged and neutral states.32 The study mentions the importance of migration of water molecules in the protein core due to the ionization of internal residues and emphasizes the need of constant pH molecular dynamics simulations (CpHMD) in explicit solvent to explain the pH-dependent structural responses of these proteins. CpHMD in explicit solvent using continuous protonation models have been carried out previously using four variants of SNase with lysine buried in different internal positions of the proteins, V66K, V99K, L125K, and I92K.33 The study emphasizes the importance of adequate sampling of “open-state” structures for the accurate prediction of pKa values of internal ionizable residues. However, it can be challenging to obtain fully converged ensembles using CpHMD as the simulations have a tendency to get trapped in local minima.15,34 Hence, CpHMD might not provide sufficient sampling of the conformational microstates traversed by the protein variants with changing pH and protonation states. Other hybrid methods35 have been used previously to predict the pKa values of internal ionizable residues in SNase mutants, taking into account the structural modifications of the ionizable residues that undergo large conformational changes. In this method, representative structures using protonated and deprotonated side chains of internal ionizable residues have been used to predict their pKa values. The rigid backbone approach within the multiconformation continuum electrostatics (MCCE) method has been used to report the pKa values of other residues. Zheng and Cui have lately carried out explicit molecular dynamics simulations with fixed charge and polarizable force fields to explore the microscopic mechanisms determining the titration response of internal ionizable residues of SNase.36 This study indicates that tweaking the dielectric constants of the protein interiors gives good agreement with experimentally proposed apparent pKa of internal ionizable residues. However, according to the experimental studies of L25K and L125K, the anomalous properties of the buried lysine residues are more to do with free energy changes associated with the pH-sensitive conformational changes undergone by the proteins and less to do with tweaking the polarizability and polarization of the protein interior surrounding the internal ionizable residue. It has been previously shown that structure based electrostatic calculations with low dielectric constants fail to reproduce the electrostatic interactions in proteins.37,38 Moreover, if there are peptides for instance, with pKa values close to 7.0, then using constant protonation states are destined to fail and the intermediate conformations will not be visible. Computational methods used to calculate pKa by assigning single protonation states to titratable residues do not include the possibility of electrostatic interactions with neighbors. If the pH lies in the vicinity of the pKa of titratable residue, or if the system lies in certain conformations in which the dominant protonation state is different, then these simulations are not necessarily accurate to describe the system of interest.14,15,39,40 In the light of the experimental insights of the pHdependent conformational changes accompanying the titration

water migration into their interiors and alter the dielectric environment around the key internal residue, without the need of extensive structural rearrangement.22 Proteins can have different conformational microstates with different occupational probabilities. If these conformational microstates are significantly different from each other in terms of their structures, the environment surrounding these ionizable groups are different because of the different orientations of charges and dipoles around them and also because of the presence or absence of water molecules in their vicinity. This pH-dependent statistical mechanical picture of protein residues has been described well in the thermodynamic linkage theory of Wyman23 and the work by Tanford.24 Frankel and Bevilacqua have lately studied the wavy rate−pH profiles for understanding how ribozymes respond to pH, when reverse protonated charged nucleobases interact.25 This study has uncovered salient features of the pH-dependent enzyme kinetics correlating the cooperativity of the general acid and general base in ribozymes and their “dark” pKa shifts. In the recent past, experimental studies have been carried out to understand the cause of the anomalous pKa of ionizable residues in the protein interiors. In these studies, variants of SNase have been used as model systems because of their relatively small size. NMR spectroscopy experiments with different variants of SNase carried out by Garcia-Moreno and colleagues,26−31 studied the conformational responses of these proteins by artificially replacing internal positions with lysine, glutamate and aspartate residues. The pH-dependent difference in free energies of structural changes between folded and unfolded wild type and the corresponding variants and their changes under similar pH conditions, have been used to report the experimentally observed apparent pKa of the internal residues. A majority of these variants have observed pKa values that deviate significantly from their intrinsic values in bulk water. Though some of the variants may tolerate the presence of charged residues due to surroundings with relatively high dielectric constants, the experiments conclude that structural rearrangements and conformational changes undergone by the protein play a significant role in the accurate pKa calculation of internal ionizable residues. The experimental study by Garcia-Moreno and collegues31 with SNase variants, L25K and L125K, has pointed out that the anomalous pKa values of the two internal lysine residues are a consequence of the diverse pH-dependent structural reorientations undergone by these proteins in a wide range of time scales. The work emphasizes the sensitivity of the observed pKa of the internal lysine residues to the nature and the time scale of conformational reorganizations that the mutant protein undergoes. Sophisticated experimental techniques like NMR relaxation dispersion spectroscopy and NMR exchange spectroscopy were used to characterize the nature, diversity and the wide range of time scales of the pHdependent conformational reorganizations undergone by L25K and L125K. According to these experiments, L25K undergoes pH-dependent fluctuations in time scales ranging from hundreds of microseconds to hundreds of milliseconds showing signatures of both local and subglobal conformational changes. Local conformational changes refer to small structural changes undergone by the protein affecting a majority of the residues in the neighborhood of the internal lysine. Subglobal conformational changes allow the exposure of the side chain of the lysine residue to bulk water, typically involving more residues. In contrast to L25K, L125K has been observed to 5743


Article


the structure was minimized using a restraint of 10 kcal mol−1 Å−2; however, only on the backbone atoms of the residues. The minimized structure was heated at constant volume from 100 to 300 K for 6 ns using the Langevin thermostat47 with a friction coefficient of 5 ps−1 and a restraint of 1 kcal mol−1 Å−2 on the backbone atoms, including stabilization at 300 K for the last 1 ns. The system was relaxed slowly in four steps. The temperature was fixed at 300 K using the Langevin thermostat with a friction coefficient of 1 ps−1, and the pressure was fixed at 1 atm using the Monte Carlo barostat. Initially, the structure was relaxed for 5 ns with a restraint of 1 kcal mol −1 Å−2 on the backbone atoms at constant temperature and pressure to ensure density stabilization. This was followed by a 5 ns relaxation at constant temperature and pressure without any restraint. The system was further relaxed for 5 ns at constant temperature and volume without any restraint. The relaxation dynamics, so far, were carried out at constant protonation state: all the lysine and histidine residues remained protonated, and all the glutamate and aspartate residues remained deprotonated throughout the relaxation dynamics. Finally, this system was used as a starting structure for the relaxation dynamics at each of the 22 equally spaced pH replicas at intervals of 0.5 units in the pH range 2.0 to 12.5. The final relaxation dynamics of the system were carried out at each of the 22 pH, in constant volume and temperature for 5 ns without any restraint, and the protonation change of the titratable residues was attempted every 100 steps with 200 steps of solvent relaxation after a successful protonation change. pH-REMD Simulations. pH-REMD simulations were carried out for 50 ns with 22 pH replicas from 2.0 to 12.5 at equally spaced intervals of 0.5 pH units for both L25K and L125K, titrating only the internal lysine residue in each case. We also did 80 ns simulations for both the variants, with 22 equally spaced pH replicas, titrating all the lysine, aspartate, histidine and glutamate residues. The last 70 ns of these simulations were used for analysis. At least five independent simulations of each kind were carried out. All the simulations were done in constant volume and temperature with a time step of 2 ns and without any restraints. The temperature was maintained at 300 K using the Langevin thermostat47 with a friction coefficient of 1 ps−1. One hundred steps of dynamics were allowed between replica exchange attempts; 100 steps were allowed between protonation change attempts, and 200 steps of solvent relaxation were performed after a successful protonation change. Random seeds were used to avoid synchronization artifacts in the simulations.48 The electrostatic interactions were treated with the particle mesh Ewald (PME) method using an 8 Å van der Waals cutoff. The GB model49 was used to attempt changes in the protonation state. The salt concentration of the solution was set to 0.1 M for all simulations, in consistence with the reference compound. All the simulations were run in GPUs using the pmemd.cuda module in AMBER 16. Analysis. Cphstats in AmberTools was used to calculate the resulting protonation state distribution from pH-REMD simulations. Titration data corresponding to a pH-ensemble was reconstructed from a replica time-series. The ensemble of protonation state change was extracted for every single pH for each of the 22 pH values simulated. The trajectories were analyzed using cpptraj50 and VMD.51 The observed pKa value

response in L25K and L125K, we used GPU-powered pH Replica Exchange Molecular Dynamics simulations (pHREMD) in explicit solvent in AMBER 16 to study our system of interest. All the calculations of the present work have been done in the Blue Waters. pH-REMD is an efficient technique coupling constant pH molecular dynamics simulations (CpHMD) and replica exchange molecular dynamics (REMD) simulations.15,41 pH-REMD enhances conformational sampling, spanning the pH phase space and accelerates the convergence of simulations with accurate pKa prediction. The details of CpHMD and pH-REMD are described in detail in references.15,39,42,43 In the present work, we study the pH-dependent conformational changes associated with the titration response of the internal lysine residues in two different internal positions of SNase, in the variants L25K and L125K. We use pH-REMD calculations to unveil the cause of the apparent pKa of these residues. The choice to specifically work with L25K and L125K is based on the detailed experimental insights into the systems by Garcia-Moreno and coauthors.31 The comparative investigation of L25K and L125K sheds light into the variety of the structural responses undergone by a protein as a result of the diversity of the dielectric environments encountered by the ionizable residue across different regions of the protein. Our work characterizes the pH-dependent fluctuations undergone by the protein variants, in consistence with the diverse experimental reorganization events proposed. More importantly, the thermodynamic models proposed explain the association of the anomalous pH-dependent titration responses of the protein variants, observed by experimental probes to the microscopic conformation-specific pKa of the internal residue attained at different pH. These thermodynamic models of various complexities have wide applications, beyond the atomistic understanding of the pH-dependent protein electrostatics associated with buried ionizable residues. They provide a generalized and robust methodology to predict the conformational-specific pKa values associated with pH sensor proteins that show conformational changes with slight changes in pH. This work provides theoretical tools for aiding the engineering of pH-sensitive conformational switches.

■

METHODOLOGY General Setup. The starting structure was prepared using the crystal structure of the Δ+PHS variant of SNase26 obtained from the Protein Data Bank (ID 3BDC). The Δ+PHS is a stable variant of SNase relative to the wild type (WT), having two substitutions G50F, V51N, and a loop deletion (residues 44−49). The ligands were not taken into account while preparing the starting structure. The L25K and L125K variants were created by replacing leucine (L) with lysine (K) in the 25th and 125th positions of the respective proteins. tleap in AmberTools44 was used to generate the starting topology and the input coordinate files. The systems were parametrized using the Amber ff99SB force field45 for all the protein residues. TIP3P was used to add water molecules explicitly into the simulation and the system was neutralized using counterions. All distances between the hydrogen atoms and heavy atoms were constrained using the SHAKE algorithm.46 Relaxation of the System. The system was minimized in two steps. Initially, the system was minimized with a restraint of 10 kcal mol −1 Å−2 on all the atoms of the protein residues using 1000 steps of steepest descent algorithm, followed by 1000 steps of conjugate gradient algorithm. In the second step, 5744


Article

The Journal of Physical Chemistry B of Lys25 and Lys125 was obtained by fitting the titration data to the Hill equation fd =

conformation-specific pKa values of internal lysine residues corresponding to buried and water-exposed conformations. This manuscript unveils the difference in the atomistic mechanisms adapted by Lys25 and Lys125 to account for their macroscopic anomalous behaviors. Furthermore, the pHsensitive fluctuations in L25K and L125K proteins have been characterized using the root-mean-square fluctuations of the backbone residues of the protein variants. L25K. Observed pKa of Lys25. The 25th and 125th locations in the wild type of SNase are occupied by leucine residues. These two positions in the wild type have entirely different neighboring residues and an inherent difference in water exposure. In the L25K variant, the leucine in the 25th location has been substituted with a lysine. Similarly, in the L125K variant, the leucine in the 125th location has been substituted with a lysine. Lys25 is located in the β2 strand of L25K, whereas Lys125 is located in the α3 helix of L125K (Figure 1). Owing to their locations, Lys25 and Lys125 encounter different surrounding electrostatic microenvironments and different accessibility to water molecules, when the respective proteins are in their native states. From our calculations, Lys25 has a pKa of 6.74 ± 0.31 and a Hill coefficient of 1.01 ± 0.08, when we allow only Lys25 to titrate. Lys25 has a pKa of 6.43 ± 0.20 and a Hill coefficient of 0.88 ± 0.10, when all the ionizable residues in the protein are allowed to titrate. Figure 2 shows a titration curve of Lys25 obtained from one independent pH-REMD simulation, when all the titratable residues are allowed to change their protonation states. It is important to note that the calculated pKa and the Hill coefficients have very close values in the two sets of calculation. The Hill coefficient which is a measure of cooperativity, is sufficiently close to unity for both the calculations. We conclude that Lys25 titrates independently, without significant contributions from its neighboring titratable residues. The observed pKa of Lys25 from our calculation is substantially depressed with respect to its intrinsic pKa of 10.4 in bulk water. Our results are in very good agreement with the experimentally observed apparent pKa value of 6.30.31 Water Molecules around the Lys25 Residue. We examined the trajectories of L25K at different pH values to check for the proximity of water molecules around the Lys25 side chain. At pH 3.0, the Lys25 side chain is protonated. In this charged state Lys25 adopts an orientation which makes it waterexposed during majority of the time (Figure 3A). In the waterexposed conformation, the lysine side chain is expected to have a microscopic conformation-specific pKa close to its intrinsic pKa of 10.4 in bulk water. At pH 12.0, Lys25 is deprotonated and neutral. At this pH, the Lys25 side chain remains snugly buried inside the hydrophobic pocket of the protein, reducing its wateraccessibility (Figure 3B). In its buried conformation within the protein moiety, the microscopic conformation-specific pKa of Lys25 is expected to be significantly lower than its intrinsic pKa of 10.4 in bulk water. To reveal a quantitative picture of our observation, we plotted the normalized probability distribution of the number of water molecules within 3.4 Å of the NZ terminal of the Lys25 side chain, at pH 2.0, 6.5, and 12.0, respectively. At pH 2.0, the probability density corresponds to a maximum of four water molecules around the Lys25 side chain (Figure 4). The distribution curve is skewed to the right, with a mean of about five water molecules around the side chain. This indicates the existence of a higher fraction of water-exposed snapshots of

1 1 + 10

n(pK a − pH)

(1)

where fd is the deprotonated fraction of the titratable residue and n is the Hill coefficient, a measure of cooperativity effects52 for residues interacting with their neighbors.

■

RESULTS AND DISCUSSIONS In this manuscript, we report the computed pKa of Lys25 and Lys125 side chains, corresponding to the midpoints of the

Figure 1. Location of Lys25 (blue sphere) and Lys125 (red sphere) in SNase. Lys25 is located in the β2 strand, and Lys125 is located in the α3 helix. SNase has been colored according to its secondary structure.

Figure 2. Titration curve of Lys25 obtained from one independent pH-REMD simulation, when all the titratable residues are allowed to change their protonation states. Data from the simulations (red dots) have been fitted to eq 1. In this particular calculation, Lys25 titrates with pKa of 6.42 and has a Hill coefficient of 0.87.

titration curves, obtained from pH-REMD calculations. We use the solvent accessible surface area (SASA) of the NZ terminal of the lysine side chain to measure the degree of burial of the titrating internal lysine with pH. In consistence with thermodynamic laws, we report estimates of the microscopic 5745


Article


Figure 3. (A) Snapshot of the exposed NZ terminal of Lys25 at pH 3.0. Lys25 prefers to remain water-exposed when protonated. (B) Snapshot of the buried NZ terminal of Lys25 at pH 12.0. Lys25 prefers to remain buried within the hydrophobic pocket when deprotonated.

Figure 6. Four-state thermodynamic cycle with the equilibrium constants and four chemical species of Lys25 in their protonated and deprotonated states with both buried and water-exposed conformations, respectively. Figure 4. Normalized probability density distribution of the number of water molecules within 3.4 Å of the NZ terminal of the Lys25 side chain at three pH values. At low pH (red), Lys is mostly solventexposed; at intermediate pH (green), it fluctuated between buried and water-exposed conformations, and at high pH (blue), it is mostly buried. The solid lines have been added to aid the eye.

Figure 7. Fraction of buried conformation as a function of pH from one independent pH-REMD simulation. The pKa of 6.18 corresponds to the midpoint of the curve. The n value of 1.02 signifies that Lys25 is primarily responsible for the change in conformation of the protein.

Lys25 at a low pH of 2.0. At pH 12.0, the probability density distribution has a relatively sharper peak corresponding to a maximum of one water molecule around Lys25. The mean of the distribution curve is shifted to the left as compared to that at pH 2.0. This behavior of the distribution curve is consistent with the decreased water-accessibility of the buried lysine side chain at a higher pH. Interestingly at intermediate pH of 6.5

Figure 5. Radial distributive function of the oxygen atom of water and the NZ terminal of Lys25 side chain at low (red), intermediate (green), and high (blue) pH.

5746


Article


occurrence of less than five water molecules within 3.4 Å of the NZ terminal of the Lys25. The data leads us to infer that at pH 6.5 the conformation ensemble of L25K is composed of a weighted fraction of the buried and water-exposed conformations of Lys25. It is important to point out that interpreting the observed anomalous pKa of the internal lysine of L25K only in terms of the change in free energies of deprotonated and protonated states is an insufficient description. We also need to account for the free energy differences associated with the pHdependent transition between the buried and water-exposed conformations. The data indicates that the observed pKa of the internal lysine is a result of the weighted probability distribution of the pH-sensitive conformations attained. To consolidate our inference, we plotted the radial distribution function of the oxygen atom of the water molecules within a radius of 10 Å of the NZ terminal of the Lys25 side chain. Consistent to our expectations, the radial distribution function (Figure 5) of the water molecules at pH 6.5 is intermediate between that at pH 2.0 and pH 12.5, which maximizes at a distance of 3 Å from the NZ terminal of Lys25. Simplest Required Thermodynamic Cycle. Looking at our data, we know that, in L25K, the titration of the Lys25 changes its protonation state and is accompanied by a change in the degree of burial of the Lys25 side chain. At low pH. when Lys25 is protonated, the water-exposed conformation dominates. At high pH, when Lys25 is deprotonated, the buried conformation dominates. To explain the energetics of this coupled, pH-dependent phenomena, we need at least four freeenergy terms. Following Tanford’s original idea,24 we present a four state thermodynamic cycle to explain the energetics of the L25K protein in association with the observed pKa of Lys25. Our group53−55 and others31,33 have used this idea in previous works to explain the coupled energetics in other systems. The relationship between the pH-dependent equilibrium constant for conformational change and the equilibrium constant for protonation and deprotonation of the internal lysine has been described using the simplest required four-state thermodynamic cycle shown in Figure 6. In the thermodynamic cycle, BuH+ and Bu represent the protonated and deprotonated states respectively of the buried lysine side chain, while ExH+ and Ex represent the protonated and deprotonated states respectively of the water-exposed lysine. The equilibrium constants for titration of the buried and the water-exposed conformations of the internal lysine are

Figure 8. Mole fraction of the four chemical species of Lys25 at different values of pH. The points with the error bars represent the mean and the standard deviation respectively, of the population fraction of the four chemical species of Lys25, obtained from five independent simulations. The solid lines have been obtained from the global fitting of eqs 4−7 where KD = 19.00, pK a(buried) = 5.96, and pK a(buried) = 7.66 are the three parameters.

K a(buried) =

[Bu][H+] [BuH+]

and K a(exposed) =

[Ex][H+] , [ExH+]

respec-

tively. These equilibrium constants are associated with the free energy difference between the protonated and deprotonated states at respective conformations. KH = KD =

[Bu] [Ex]

[BuH+] [ExH+]

and

correspond to the equilibrium constants determined

by the change in free energy associated with the structural transition between buried and the water-exposed conformations when lysine is protonated and deprotonated, respectively. The simple thermodynamic cycle holds true for Lys25 because our calculations show that its titration behavior remains unaffected by the titration behavior of its neighboring titratable residues. At low pH, Lys25 is protonated and mostly water-exposed and hence ExH+ is the major chemical species present in the system. For this to be true, KH < 1. At high pH however, Bu is present in majority, as lysine is deprotonated and in its buried conformation. For this to hold true, KD ≫ 1. According to the

Figure 9. (A) Root mean square fluctuation (Å) of the backbone atoms of the residues in L25K with respect to the average structure calculated at six pH values. (B) The heat map of the L25K protein indicating the variation of the flexibility of the protein with pH changes.

close to the observed pKa of Lys25, the probability density distribution has two peaks which seems to be an average of the ones at pH 2.0 and pH 12.0. Figure S1 shows the probability of 5747


Article


Figure 10. L25K colored according to the RMSF at pH values (A) 5.5, (B) 6.0, (C) 7.0, and (D) 8.0. Red symbolizes the highly fluctuating regions and blue symbolizes the inflexible regions in the protein.

Table 1. Observed pKa and the Hill Coefficient of Lys125: Dependence on the Protonation State of His121 protonation state of His121

observed pKa Of Lys125

Hill coefficient (n)

neutral charged titrating

7.00 ± 0.70 6.46 ± 0.80 7.20 ± 0.11

1.09 ± 0.07 1.06 ± 0.04 0.48 ± 0.04

thermodynamic cycle, the pH-dependent conformational change from the chemical species ExH+ at low pH to Bu at high pH can happen either via the path with BuH+ as the intermediate chemical species or via the path with Ex as the intermediate species. Since the total change in free energy (free energy for change in protonation state and conformation change) around the cycle is independent of the path, we have KHK a(buried) = KDK a(exposed)

(2)

Therefore, there are only three independent equilibrium constants, we need to account for to explain the present case of the anomalous observed pKa of Lys25. To quantify the extent of burial of Lys25 with pH change, we use a criterion based on the solvent accessible surface area (SASA) of the Lys25 side chain, similar to the treatment used in ref 55. The SASA of the NZ terminal of the lysine side chain has been calculated using a probe radius of 1.4 Å at each pH. The SASA histograms (Figure S2) show that the mean SASA of Lys25 reduces with the increase of pH. At pH 2.0 the conformational ensemble is dominated by water-exposed snapshots of Lys25, with widespread SASA values ranging

Figure 11. Normalized probability distribution of the number of water molecules within 3.4 Å of the NZ terminal of the Lys125 side chain. At low pH (red), the Lys residue is water-exposed, at intermediate pH (green), it fluctuates between buried and waterexposed conformations and at high pH (blue) it is mostly buried. The lines have been added to aid the eye.

from 0 to 70 Å2. At pH 12.5 the ensemble consists of a majority of buried conformations of Lys25, with a SASA values ranging from 0 to 30 Å2. A cutoff value of 15 Å2 has been used to separate the sampled conformations of L25K into buried and water-exposed. All the conformations with SASA lesser or 5748


Article


Figure 12. Radial distributive function of the oxygen atom of water and the NZ terminal of Lys125 side chain at low (red), intermediate (green), and high (blue) pH.

Figure 14. (A) Root mean square fluctuation (Å) of the backbone atoms of the residues in L125K with respect to the average structure calculated at six pH values. (B) The heat map of L125K indicating the variation of the protein’s flexibility with changes in pH.

Figure 13. Eight-state thermodynamic network. The network relates the 12 equilibrium constants with the eight chemical species of Lys125 and His121, in their protonated and deprotonated states in both buried and water exposed conformations.

This value is very close to its experimentally observed apparent pKa of 6.30. The Hill coefficient obtained from the buried fraction vs pH sigmoid is 0.89 ± 0.24, similar to that obtained from the titration curve. The experimentally observed apparent pKa of the internal lysine is a consequence of the weighed probability distribution of the pH-dependent conformations attained by Lys25. These anomalous pKa values observed in the internal ionizable residues, driven by pH-dependent conformational changes lead to conformation-specific pKa values via pH-dependent allosteric changes.24 Relationship between Conformation-Specific pKa and the Observed pKa of Lys25. According to the law of mass conservation, the sum of the mole fractions of all the chemical species in a system is one. Using the proposed four state thermodynamic cycle, we can obtain analytical expressions for the four chemical species of Lys25 in terms of the three equilibrium constantsKD, Ka(exposed), Ka(buried), and the hydrogen ion concentration in the system24,31,53,54

equal to 15 Å2 were considered to be buried. Likewise, all the conformations with SASA larger than 15 Å2 were considered to be water exposed. The fraction of buried conformation at each pH has been fitted to eq 3 (Figure 7). f = f min +

f max − f min 1 + 10 n(pK a − pH)

(3)

where fmin corresponds to the fraction of buried conformation at pH 2.0, fmax corresponds to the fraction of buried conformation at pH 12.5 and n is the Hill coefficient. The pKa value obtained from eq 3 corresponds to the midpoint of the SASA versus pH sigmoid. Looking at Figures 2 and 7, in the pH range 2.0−4.0 when Lys25 is protonated, the ensemble consists of about 25% of the buried conformation. In the pH range 9.0−12.5 when Lys25 is completely deprotonated, the ensemble consists of roughly 95% of buried conformation. The pH value corresponding to the midpoint of the buried fraction versus pH sigmoid of Lys25 is 6.33 ± 0.25, very close to the observed pKa of 6.43 ± 0.20, obtained from the titration data of our pH-REMD calculations.

[Ex] =

1 +

1 + KD + [H ] 5749

(

KD K a(buried)

+

1 K a(exposed)

)

(4)


Article


Figure 15. L125K colored according to RMSF at pH values (A) 5.5, (B) 6.0, (C) 7.0, and (D) 8.0. Red symbolizes the highly fluctuating regions and blue symbolizes the inflexible regions in the protein.

to our thermodynamic model, the observed pKa can be expressed as

[ExH+] +

[H ]

=

(

K

D + K a(exposed) 1 + KD + [H+]( K (buried) a

observed pK a = α([Ex] + [ExH+]) + β([Bu] + [BuH+])

)

1 ) K a(exposed)

(8)

(5)

[Bu] =

where α and β are limiting constants in the water-exposed and buried states, respectively. Using the analytic expression of the four chemical species, eqs 4−7, the observed pKa can be associated with the microscopic conformation-specific pKa values and the equilibrium constant KD53,54

KD 1 + KD + [H+]

(

KD K a(buried)

+

1 K a(exposed)

)

(6)

[BuH+] =

observed pK a = pK a(buried) + pK a(exposed)

[H+]KD

(

K

D + K a(buried) 1 + KD + [H+]( K (buried) a

1 ) K a(exposed)

− log10

)

KD + 1 KDK a(exposed) + K a(buried)

(7)

(9)

At low pH, Lys25 is protonated and water-exposed (implying ExH+ and BuH+ are the only species present and ExH+ is the dominant species), and hence, KH = 0.31 ± 0.10 was estimated from the fraction of buried conformation at pH 2.0 (0.23 ± 0.06). Similarly, KD = 61.00 ± 40.10 was estimated from the fraction of buried conformation at pH 12.5 (0.98 ± 0.02), where Bu and Ex are the only two species present, with Bu being the dominant one. It is important to point out that owing to the definition of KD, a small change in the fraction of the buried confirmation at pH 12.5 leads to a big leap in the estimated value of KD. The experiments involving the determination of the pKa of the Lys25 side chain measure the change in the free energy of the overall conformation of the protein and are dependent on the mole fraction of the chemical species of Lys25. According

Combining eqs 9 and 2, we can express the observed pKa in terms of the constantsKD, KH , andK a(exposed) observed pK a = pK a(exposed) − log10

KD + 1 KH + 1

(10)

Using the estimated values of KH = 0.31 ± 0.10, KD = 61.00 ± 40.10, and the observed pKa of 6.43 ± 0.20 for Lys25 from pHREMD calculations (which is very close to the experimentally predicted value of 6.30) in eq 10, we get estimates of the microscopic pKa values to be 5.79 ± 0.29 in the buried conformation and 8.02 ± 0.43 in the water-exposed conformation. Figure 8 shows the mole fraction of the four chemical species of Lys25 plotted as a function of pH. The points with the error bars represent the mean and the standard deviation 5750


Article


one for an ionizable residue titrating independently, which might indicate the interference of Lys125 with the titration behavior of a neighboring residue. We observed that among all the titrating residues within the vicinity of the Lys125 side chain, His121 is the only one showing a significantly low Hill coefficient of 0.48 ± 0.08. Interestingly, when the protonation state of His121 is kept fixed (both neutral and charged), the Hill coefficient of Lys125 is very close to one. The data from pH-REMD simulations show that the titration curve and hence the observed pKa of Lys125 is dependent on the protonation state of His121. The calculated pKa values of Lys125 and the corresponding protonation states of His121 have been reported in Table 1. Similar to Lys25, the observed pKa of Lys125 is significantly depressed with respect to its intrinsic pKa of 10.4 in bulk water. Water Molecules around the Lys125 Residue. Similar to Lys25, we plotted the normalized probability distribution of the number of water molecules within 3.4 Å of the NZ terminal of the Lys125 side chain at pH 2.0, 7.0, and 12.5. At pH 2.0, the probability density distribution of the number of water molecules around Lys125 (Figure 11) has a peak corresponding to four water molecules, indicating an increased waterexposure of Lys125 side chain when protonated. At this pH, Lys125 is expected to have a conformation-specific pKa close to its intrinsic pKa of 10.4. At pH 12.5, when Lys125 is deprotonated, the probability density attains a major peak with one water molecule around the Lys125 side chain. A small peak is seen corresponding to about three water molecules around the Lys125 side chain. At this high pH, Lys125 is buried and is expected to have a depressed conformationspecific pKa with respect to its intrinsic pKa of 10.4. At intermediate pH 7.0, Lys125 has a weighted distribution between buried and water-exposed conformation. The radial distributive function of the water molecules around the Lys125 side chain (Figure 12) also supports the fact that the water distribution around Lys125 at pH 7.0 is intermediate between that at pH 2.0 and pH 12.5 which maximizes at a distance of around 3 Å from the NZ terminal. Eight-State Thermodynamic Network. Our data suggests that the titration behavior of Lys125 is coupled to the titrating His121 residue. Also, their protonation states are affected by at least two conformations each, that is, buried and water exposed. The four-state thermodynamic cycle is not sufficient to define the pH-dependent coupled phenomena of protonation and conformation change driven by at least two key residues in a protein. To explain this relatively complex phenomena, the simplest model would be an eight-state thermodynamic network (Figure 13). The eight-state network which is an extension of the four-state thermodynamic cycle, can successfully explain the energetics of L125K.23 In the eight-state thermodynamic network, the front face represents the four chemical species of the Lys125 and His121 pair when the His121 is neutral, represented by K+BH0B, K+EH0E, K0BH0B, and K0EH0E. Similarly, the back face of the thermodynamic network represents the four chemical species of the Lys125 and His121 pair, when His121 is charged, represented by K+BH+B, K+EH+E, K0BH+B, and K0EH+E. At equilibrium, the His121 side chain and Lys125 side chain can either remain simultaneously buried or simultaneously water-exposed, eliminating the possibility of the existence of eight more combinations of the Lys−His pair. The 12 equilibrium constants corresponding to the edges of the cubic network are associated either with the free energy

respectively, of the population fraction of the four chemical species of Lys25 at each pH values obtained from five independent simulations. It is interesting to note that due to the nature of the limiting behavior of the four analytical functions (eq 4−7), the difference in the two microscopic pKa values has to necessarily have a tight bound, despite the large uncertainty in the value ofKD. The data from the simulations are consistent with this expected behavior of the analytical functions. The solid lines have been obtained by the global fitting of the four analytical functions, eqs 4−7, which give us the three parameters, KD = 19.00, pK a(buried) = 5.96, and pK a(buried) = 7.66. In Figure 8, for pH values less than 6.0, the protonated lysine is the dominant chemical species with the majority of them being in the water-exposed conformation. At pH 8.0 and beyond, Lys25 is completely deprotonated with an insignificant existence of the water-exposed species. The results obtained are consistent with the experimental observations31 according to which Lys25 is neutral above pH 7.2 and charged below pH 5.8. Mobility and pH-Dependent Conformational Changes in L25K. To evaluate the flexibility of L25K, we calculated the root mean square fluctuation (RMSF) of the backbone atoms of all the residues, with respect to the average structure at six values of pH (Figure 9A). In L25K the highly fluctuating regions are spread throughout the entire protein and are not limited to the vicinity of the internal lysine. Residues 4−13 (located in the β1 sheet), 16−30 (located in the β2 and β3 sheets), and 100−110 (located in the turn before the α3 helix) are highly flexible as displayed by the corresponding RMSF values. Residues 50−65 and 80−90 show intermediate mobility in comparison. The regions of high and intermediate flexibility, distributed throughout the entire β barrel correspond to subglobal conformational changes in L25K, as noted by the experimentalists. The simulations show that the entire β barrel moves in a coordinated manner in response to the titration of Lys25, corresponding to the experimentally detected longer time scale conformational changes.31 Figure 9B depicts a heat map of the pH-dependent flexibility of L25K. Residues 20−29, 50−70, and 100−110, all located in or near the β barrel and the α helix, show pH-dependent fluctuations of varying magnitudes. In the vicinity of Lys25, residues 15−30 fluctuate with a higher intensity at a low pH which corresponds to the increase in fluctuations around the β1 and other sites of the β sheet at acidic pH. These conformational changes around β1 sheet at low pH lead to the exposure of Lys25 to water molecules, leading to the protonation of the internal lysine. Some pH-dependent fluctuations are also observed in the region around the α3 helix. The signatures of the flexibility of the L25K protein are in agreement with experiments.31 Figure 10 shows the regions of the L25K protein that undergo pH-sensitive conformational changes and their variation with increasing pH. It can be seen that the flexibility of the protein decreases with increase in pH, evident by the change in the intensity of the red color in the β strand. L125K. Observed pKa of Lys125. Similar to L25K, we carried out two set of titration calculations for L125K. We observed that Lys125 has a Hill coefficient of 1.09 ± 0.07, when we allow only Lys125 to titrate. However, Lys125 has a significantly depressed Hill coefficient of 0.48 ± 0.04, when all the ionizable residues in the protein are allowed to titrate. This value is highly depressed with respect to the expected value of 5751


Article


5.55 ± 0.99, and pKE+ a (K) = 8.18 ± 0.18, from the calculations when His121 is neutral and charged, respectively. Mobility and the pH-Dependent Conformational Changes in L125K. In L125K, the high fluctuating regions are located around residues 100−120, indicating the high mobility of the α3 helix (Figure 14A). These flexible regions in L125K are primarily localized in the vicinity of Lys125, although some less intense fluctuations are seen throughout the protein (Figure 14B). The structural fluctuations observed in the L125K protein are relatively pH-independent as seen in Figure 15. These events are in coherence with the experimentally observed pH-independent reorganizations in the microsecond−millisecond time scale near Lys125.

difference of protonation change or free energy difference of E0 conformation change. KB0 a (K) and Ka (K) are associated with the free energy difference between the protonated and deprotonated states of Lys125, when His121 is deprotonated buried and deprotonated water-exposed, respectively. Similarly, E+ KB+ a (K) and Ka (K) are associated with the free energy difference between the protonated and deprotonated states of Lys125, when His121 is protonated buried and protonated water-exposed, respectively. These equilibrium constants give the conformation-specific pKa values of Lys125 at two fixed protonation states of His121. Similarly, we can define KB0 a (H), B+ E0 KE+ a (H), Ka (H), K a (H), where His121 is the titrating residue with two respective fixed protonation states of Lys125. The equilibrium constants K DD , K DH , K HD , and K HH correspond to the free energy difference associated with the simultaneous conformational change from buried to waterexposed for the Lys125 and His121 pair, at different protonation states. The observed apparent pKa of L125K as measured by the experimentalists corresponds to the Ka (eff), in our eight-state thermodynamic model. It is important to note that this is not the same as the observed pKa of Lys125 from pH-REMD calculations. Following the laws of thermodynamics, the total free energy change associated with the transition from K+EH+E at low pH to K0BH0B at high pH in the thermodynamic network is path-independent. There are only six routes that can take us from K+EH+E to K0BH0B. Therefore, for a given Ka(eff), only six equilibrium constants are independent and sufficient to define the whole thermodynamic network. It is important to emphasize that any randomly chosen six equilibrium constants may not be independent. Also, the choice of these six independent equilibrium constants is not unique. For example, E0 the equilibrium constants KE0 a (K), KDD, KHD, Ka (H), KDH, B+ and Ka (K), represent a set of six independent equilibrium constants, sufficient to define the whole thermodynamic network. Estimating Useful Constants from the Eight-State Thermodynamic Network. Keeping in mind that at least two key residues trigger the pH-dependent conformational change in L125K, leading to the anomalous observed experimental pKa, we use a modified criterion to separate the conformational ensemble of L125K and hence to estimate the value of Ka(eff). We use a combined criterion based on the SASA of the Lys125 side chain and the His121 side chain. The SASA of the NZ terminal of the lysine side chain and the histidine side chain has been calculated using a probe radius of 1.4 Å. A cut off SASA of 6 Å2 corresponding to lysine, and a cut off SASA of 60 Å2, corresponding to histidine have been used to separate the conformational snapshots into fully buried and fully water-exposed conformations. All the structures in the ensemble that had a lysine SASA value lesser or equal to 6 Å2 and simultaneously a histidine SASA value less or equal to 60 Å2 were considered to be a part of the fully buried conformational ensemble. Likewise, all the structures that had a lysine SASA larger than 6 Å2 and a histidine SASA value larger than 60 Å2 were considered to be a part of the waterexposed conformational ensemble. When the fraction buried of the conformational ensemble of L125K is plotted as a function of pH, the calculations give us an estimate of the pKa(eff) = 5.78 ± 0.29. This value is in good agreement to the experimentally reported value of 6.20 for Lys125. Finally, we estimate the microscopic conformational specific values to be, B0 B+ pKE0 a (K) = 8.78 ± 0.11, pKa (K) = 5.92 ± 0.59, pKa (K) =

■

CONCLUSION

■

ASSOCIATED CONTENT

The present work explores the pH-dependent conformational changes of ionizable lysine residues Lys25 and Lys125, buried in the hydrophobic interior of L25K and L125K variants of SNase, respectively, in association to their experimentally observed anomalous pKa values. To study the problem of interest, we calculated the pKa values of both the residues using the GPU powered pH-REMD methodology. The calculated observed pKa values displayed by these internal lysine residues (from the midpoint of the titration curve) were established to be a consequence of the pH-dependent structural reorganizations undergone by the titrating residues to accommodate the varying microenvironments in their vicinity. In L25K, Lys25 is the only key residue governing the pH-dependent conformational change. Interestingly, L125K offers a more challenging scenario, where two key residues, namely, the internal Lys125 and His121 seem to be governing the conformational changes undergone by the protein. A range of analyses were carried out to characterize the nature and pH-dependence of the conformational changes undergone. More importantly, thermodynamic models were used to establish that the observed pKa values of the internal residues are a reflection of the coupled nature of conformational equilibria and ionization equilibria of the key residues. An estimate of these microscopic conformation-specific pKa values have been proposed. The results obtained in our work can be applied to understand the atomistic mechanisms of complex proteins, possessing internal ionizable residues undergoing pH-dependent conformational changes, via allosteric changes, governing their functions. Further, the apparent pKa values of internal ionizable residues can be used to engineer pH sensing regions using an artificially introduced single residue into the hydrophobic environments of the protein that are capable of inducing conformational changes, leading to microscopic pKa values.

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.9b03816. Figures showing the probability of occurrence of water molecules near the NZ terminal of the Lys25 side chain at different pH values, normalized histogram of the solvent exposed surface area of Lys25 at low pH and high pH, titration curve of Lys125 and His121, and normalized histogram of the solvent exposed surface area of Lys125 at low and high pH and fraction of buried conformation at each pH for L125K (PDF) 5752


Article


■

(15) Swails, J. M.; Roitberg, A. E. Enhancing Conformation and Protonation State Sampling of Hen Egg White Lysozyme Using PH Replica Exchange Molecular Dynamics. J. Chem. Theory Comput. 2012, 8 (11), 4393−4404. (16) McIntosh, L. P.; Hand, G.; Johnson, P. E.; Joshi, M. D.; Korner, M.; Plesniak, L. A.; Ziser, L.; Wakarchuk, W. W.; Withers, S. G. The PK(a) of the General Acid/Base Carboxyl Group of a Glycosidase Cycles during Catalysis: A 13C-NMR Study of Bacillus Circulans Xylanase. Biochemistry 1996, 35 (31), 9958−9966. (17) Chivers, P. T.; Prehoda, K. E.; Volkman, B. F.; Kim, B. M.; Markley, J. L.; Raines, R. T. Microscopic PK(a) Values of Escherichia Coli Thioredoxin. Biochemistry 1997, 36 (48), 14985−14991. (18) Warshel, A. Energetics of Enzyme Catalysis. Proc. Natl. Acad. Sci. U. S. A. 1978, 75 (11), 5250−5254. (19) Li, Z.; Raychaudhuri, S.; Wand, A. J. Insights into the Local Residual Entropy of Proteins Provided by NMR Relaxation. Protein Sci. 1996, 5 (12), 2647−2650. (20) Jarymowycz, V.; Stone, M. Fast Time Scale Dynamics of Protein Backbones: NMR Relaxation Methods, Applications, and Functional Consequences. Chem. Rev. 2006, 106, 1624−1671. (21) Sham, Y. Y.; Chu, Z. T.; Warshel, A. Consistent Calculations of p K a ’s of Ionizable Residues in Proteins: Semi-Microscopic and Microscopic Approaches. J. Phys. Chem. B 1997, 101 (22), 4458− 4472. (22) Dwyer, J. J.; Gittis, A. G.; Karp, D. A.; Lattman, E. E.; Spencer, D. S.; Stites, W. E.; García-Moreno E., B. High Apparent Dielectric Constants in the Interior of a Protein Reflect Water Penetration. Biophys. J. 2000, 79 (3), 1610−1620. (23) Wyman, J. Heme Proteins. Adv. Protein Chem. 1948, 4, 407− 531. (24) Tanford, C. Ionization-Linked Changes in Protein Conformation. I. Theory. J. Am. Chem. Soc. 1961, 83 (7), 1628−1634. (25) Frankel, E. A.; Bevilacqua, P. C. Complexity in PH-Dependent Ribozyme Kinetics: Dark PKaShifts and Wavy Rate-PH Profiles. Biochemistry 2018, 57 (5), 483−488. (26) Castañeda, C. A.; Fitch, C. A.; Majumdar, A.; Khangulov, V.; Schlessman, J. L.; García-Moreno, B. E. Molecular Determinants of the PKavalues of Asp and Glu Residues in Staphylococcal Nuclease. Proteins: Struct., Funct., Genet. 2009, 77 (3), 570−588. (27) Chimenti, M. S.; Khangulov, V. S.; Robinson, A. C.; Heroux, A.; Majumdar, A.; Schlessman, J. L.; García-Moreno, B. Structural Reorganization Triggered by Charging of Lys Residues in the Hydrophobic Interior of a Protein. Structure 2012, 20 (6), 1071− 1085. (28) Harms, M. J.; Castañeda, C. A.; Schlessman, J. L.; Sue, G. R.; Isom, D. G.; Cannon, B. R.; García-Moreno E., B. The PKaValues of Acidic and Basic Residues Buried at the Same Internal Location in a Protein Are Governed by Different Factors. J. Mol. Biol. 2009, 389 (1), 34−47. (29) Isom, D. G.; Castaneda, C. A.; Cannon, B. R.; Garcia-Moreno E., B. Large Shifts in PKa Values of Lysine Residues Buried inside a Protein. Proc. Natl. Acad. Sci. U. S. A. 2011, 108 (13), 5260−5265. (30) Isom, D. G.; Castaneda, C. A.; Cannon, B. R.; Velu, P. D.; Garcia-Moreno E., B. Charges in the Hydrophobic Interior of Proteins. Proc. Natl. Acad. Sci. U. S. A. 2010, 107 (37), 16096−16100. (31) Richman, D. E.; Majumdar, A.; García-Moreno E., B. Conformational Reorganization Coupled to the Ionization of Internal Lys Residues in Proteins. Biochemistry 2015, 54 (38), 5888−5897. (32) Damjanović, A.; Brooks, B. R.; García-Moreno E., B. Conformational Relaxation and Water Penetration Coupled to Ionization of Internal Groups in Proteins. J. Phys. Chem. A 2011, 115 (16), 4042−4053. (33) Goh, G. B.; Laricheva, E. N.; Brooks, C. L. Uncovering PHDependent Transient States of Proteins with Buried Ionizable Residues. J. Am. Chem. Soc. 2014, 136 (24), 8496−8499. (34) Williams, S. L.; De Oliveira, C. A. F.; McCammon, J. A. Coupling Constant PH Molecular Dynamics with Accelerated Molecular Dynamics. J. Chem. Theory Comput. 2010, 6 (2), 560−568.

AUTHOR INFORMATION

Corresponding Author

*E-mail: roitberg@ufl.edu. ORCID

Adrian E. Roitberg: 0000-0003-3963-8784 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This research is part of the Blue Waters sustained-petascale computing project, which is supported by the National Science Foundation (awards OCI-0725070 and ACI-1238993) and the state of Illinois. Blue Waters is a joint effort of the University of Illinois at Urbana−Champaign and its National Center for Supercomputing Applications. A.S. was supported in part by the College of Liberal Arts and Sciences Dissertation Fellowship funded by the Charles Vincent and Heidi Cole McLaughlin Endowment, University of Florida.

■

REFERENCES

(1) Perutz, M. F.; Kendrew, J. C.; Watson, H. C. Structure and Function of Haemoglobin: II. Some Relations between Polypeptide Chain Configuration and Amino Acid Sequence. J. Mol. Biol. 1965, 13 (3), 669−678. (2) Ho, M. C.; Ménétret, J. F.; Tsuruta, H.; Allen, K. N. The Origin of the Electrostatic Perturbation in Acetoacetate Decarboxylase. Nature 2009, 459 (7245), 393−397. (3) Abrahams, J. P.; Leslie, A. G. W.; Lutter, R.; Walker, J. E. Structure at 2.8 Â Resolution of F1-ATPase from Bovine Heart Mitochondria. Nature 1994, 370 (6491), 621−628. (4) Iwata, S.; Ostermeier, C.; Ludwig, B.; Michel, H. Structure at 2.8 Å Resolution of Cytochrome c Oxidase from Paracoccus Denitrificans. Nature 1995, 376 (6542), 660−669. (5) Luecke, H.; Schobert, B.; Richter, H. T.; Cartailler, J. P.; Lanyi, J. K. Structure of Bacteriorhodopsin at 1.55 Å Resolution. J. Mol. Biol. 1999, 291 (4), 899−911. (6) Harris, T. K.; Turner, G. J. Structural Basis of Perturbed PK a Values of Catalytic Groups in Enzyme Active Sites. IUBMB Life 2002, 53, 85−98. (7) Hosler, J. P.; Ferguson-Miller, S.; Mills, D. A. Energy Transduction: Proton Transfer Through the Respiratory Complexes. Annu. Rev. Biochem. 2006, 75 (1), 165−187. (8) Bierzynski, a.; Kim, P. S.; Baldwin, R. L. A Salt Bridge Stabilizes the Helix Formed by Isolated C-Peptide of RNase A. Proc. Natl. Acad. Sci. U. S. A. 1982, 79 (April), 2470−2474. (9) Shoemaker, K. R.; Kim, P. S.; Brems, D. N.; Marqusee, S.; York, E. J.; Chaiken, I. M.; Stewart, J. M.; Baldwin, R. L. Nature of the Charged-Group Effect on the Stability of the C-Peptide Helix. Proc. Natl. Acad. Sci. U. S. A. 1985, 82 (8), 2349−2353. (10) Antosiewicz, J.; Briggs, J. M.; McCammon, J. A. Orientational Steering in Enzyme-Substrate Association: Ionic Strength Dependence of Hydrodynamic Torque Effects. Eur. Biophys. J. 1996, 24 (3), 137−141. (11) Hunenberger, P. H.; Helms, V.; Narayana, N.; Taylor, S. S.; Mccammon, J. A. Determinants of Ligand Binding to CAMPDependent Protein Kinase. Biochemistry 1999, 38 (8), 2358−2366. (12) Demchuk, E.; Genick, U. K.; Woo, T. T.; Getzoff, E. D.; Bashford, D. Protonation States and PH Titration in the Photocycle of Photoactive Yellow Protein. Biochemistry 2000, 39 (5), 1100−1113. (13) Dillet, V.; Dyson, H. J.; Bashford, D. Calculations of Electrostatic Interactions and PK(a)s in the Active Site of Escherichia Coli Thioredoxin. Biochemistry 1998, 37 (28), 10298−10306. (14) Sabri Dashti, D.; Meng, Y.; Roitberg, A. E. PH-Replica Exchange Molecular Dynamics in Proteins Using a Discrete Protonation Method. J. Phys. Chem. B 2012, 116 (30), 8805−8811. 5753


Article

The Journal of Physical Chemistry B (35) Witham, S.; Talley, K.; Wang, L.; Zhang, Z.; Sarkar, S.; Gao, D.; Yang, W.; Alexov, E. Developing Hybrid Approaches to Predict PK a Values of Ionizable Groups. Proteins: Struct., Funct., Genet. 2011, 79 (12), 3389−3399. (36) Zheng, Y.; Cui, Q. Microscopic Mechanisms That Govern the Titration Response and PKavalues of Buried Residues in Staphylococcal Nuclease Mutants. Proteins: Struct., Funct., Genet. 2017, 85 (2), 268−281. (37) Schutz, C. N.; Warshel, A. What Are the Dielectric “Constants” of Proteins and How to Validate Electrostatic Models? Proteins: Struct., Funct., Genet. 2001, 44 (4), 400−417. (38) Antosiewicz, J.; McCammon, J. A.; Gilson, M. K. The Determinants of p K a s in Proteins †. Biochemistry 1996, 35 (24), 7819−7833. (39) Swails, J. M.; York, D. M.; Roitberg, A. E. Constant PH Replica Exchange Molecular Dynamics in Explicit Solvent Using Discrete Protonation States: Implementation, Testing, and Validation. J. Chem. Theory Comput. 2014, 10 (3), 1341−1352. (40) Meng, Y.; Roitberg, A. E. Constant PH Replica Exchange Molecular Dynamics in Biomolecules Using a Discrete Protonation Model. J. Chem. Theory Comput. 2010, 6 (4), 1401−1412. (41) Itoh, S. G.; Damjanović, A.; Brooks, B. R. PH Replica-Exchange Method Based on Discrete Protonation States. Proteins: Struct., Funct., Genet. 2011, 79 (12), 3420−3436. (42) Mongan, J.; Case, D. A.; McCammon, J. A. Constant PH Molecular Dynamics in Generalized Born Implicit Solvent. J. Comput. Chem. 2004, 25 (16), 2038−2048. (43) Mongan, J.; Case, D. A. Biomolecular Simulations at Constant PH. Curr. Opin. Struct. Biol. 2005, 15 (2), 157−163. (44) Case, D. A.; Cheatham, T. E.; Darden, T.; Gohlke, H.; Luo, R.; Merz, K. M.; Onufriev, A.; Simmerling, C.; Wang, B.; Woods, R. J. The Amber Biomolecular Simulation Programs. J. Comput. Chem. 2005, 26 (16), 1668−1688. (45) Hornak, V.; Abel, R.; Okur, A.; Strockbine, B.; Roitberg, A.; Simmerling, C. Comparison of Multiple Amber Force Fields and Development of Improved Protein Backbone Parameters. Proteins: Struct., Funct., Genet. 2006, 65 (3), 712−725. (46) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of n-Alkanes. J. Comput. Phys. 1977, 23 (3), 327−341. (47) Pastor, R. W.; Brooks, B. R.; Szabo, A. An Analysis of the Accuracy of Langevin and Molecular Dynamics Algorithms. Mol. Phys. 1988, 65 (6), 1409−1419. (48) Sindhikara, D. J.; Kim, S.; Voter, A. F.; Roitberg, A. E. Bad Seeds Sprout Perilous Dynamics: Stochastic Thermostat Induced Trajectory Synchronization in Biomolecules. J. Chem. Theory Comput. 2009, 5 (6), 1624−1631. (49) Onufriev, A.; Bashford, D.; Case, D. A. Exploring Protein Native States and Large-Scale Conformational Changes with a Modified Generalized Born Model. Proteins: Struct., Funct., Genet. 2004, 55 (2), 383−394. (50) Roe, D. R.; Cheatham, T. E., III PTRAJ and CPPTRAJ: Software for Processing and Analysis of Molecular Synamics Trajectory Data. J. Chem. Theory Comput. 2013, 9 (7), 3084−3095. (51) Humphrey, W.; Dalke, A.; Schulten, K. {}: {Visual} Molecular Dynamics. J. Mol. Graphics 1996, 14 (1), 33−38. (52) Williams, D. H.; Stephens, E.; O’Brien, D. P.; Zhou, M. Understanding Noncovalent Interactions: Ligand Binding Energy and Catalytic Efficiency from Ligand-Induced Reductions in Motion within Receptors and Enzymes. Angew. Chem., Int. Ed. 2004, 43 (48), 6596−6616. (53) Di Russo, N. V.; Estrin, D. A.; Martí, M. A.; Roitberg, A. E. PHDependent Conformational Changes in Proteins and Their Effect on Experimental PKas: The Case of Nitrophorin 4. PLoS Comput. Biol. 2012, 8 (11), e1002761. (54) Di Russo, N. V.; Martí, M. A.; Roitberg, A. E. Underlying Thermodynamics of PH-Dependent Allostery. J. Phys. Chem. B 2014, 118 (45), 12818−12826.

(55) Liu, J.; Swails, J.; Zhang, J. Z. H.; He, X.; Roitberg, A. E. A Coupled Ionization-Conformational Equilibrium Is Required to Understand the Properties of Ionizable Residues in the Hydrophobic Interior of Staphylococcal Nuclease. J. Am. Chem. Soc. 2018, 140 (5), 1639−1648.

5754


pH-Dependent Conformational Changes Due to Ionizable Residues in

Recommend Documents