A Coupled Ionization-Conformational Equilibrium Is Required To

Jan 8, 2018 - Ionizable residues in the interior of proteins play essential roles, especially in biological energy transduction, but are relatively ra...
50 downloads 7 Views 3MB Size
Article Cite This: J. Am. Chem. Soc. XXXX, XXX, XXX−XXX

pubs.acs.org/JACS

A Coupled Ionization-Conformational Equilibrium Is Required To Understand the Properties of Ionizable Residues in the Hydrophobic Interior of Staphylococcal Nuclease Jinfeng Liu,†,‡,§ Jason Swails,⊥ John Z. H. Zhang,†,∥ Xiao He,*,†,∥ and Adrian E. Roitberg*,‡ †

School of Chemistry and Molecular Engineering, Shanghai Engineering Research Center of Molecular Therapeutics and New Drug Development, East China Normal University, Shanghai, 200062, China ‡ Department of Chemistry, University of Florida, Gainesville, Florida 32611, United States § Department of Basic Medicine and Clinical Pharmacy, China Pharmaceutical University, Nanjing, 210009, China ∥ NYU-ECNU Center for Computational Chemistry at NYU Shanghai, Shanghai, 200062, China ⊥ Department of Chemistry and Chemical Biology and BioMaPS Institute, Rutgers University, Piscataway, New Jersey 08854, United States S Supporting Information *

ABSTRACT: Ionizable residues in the interior of proteins play essential roles, especially in biological energy transduction, but are relatively rare and seem incompatible with the complex and polar environment. We perform a comprehensive study of the internal ionizable residues on 21 variants of staphylococcal nuclease with internal Lys, Glu, or Asp residues. Using pH replica exchange molecular dynamics simulations, we find that, in most cases, the pKa values of these internal ionizable residues are shifted significantly from their values in solution. Our calculated results are in excellent agreement with the experimental observations of the Garcia-Moreno group. We show that the interpretation of the experimental pKa values requires the study of not only protonation changes but also conformational changes. The coupling between the protonation and conformational equilibria suggests a mechanism for efficient pH-sensing and regulation in proteins. This study provides new physical insights into how internal ionizable residues behave in the hydrophobic interior of proteins.



INTRODUCTION

have demonstrated that the charge-state alteration of internal residues can trigger protein structural reorganization of varying amplitude,10,12,13 which would be a dominant factor that determines the unusual pKa values of the buried ionizable group. Many proteins harness conformational changes, which are coupled to the protonation states of internal ionizable residues, for functional purposes.14,15 This kind of biochemical process is often associated with environmental pH changes, and pH-dependent allosteric regulation has important effects on protein function.16−20 Conformational changes coupled to the ionization of the internal ionizable residues can also be reflected in their pKa values.11,21−23 For proteins that depend on internal ionizable residues for biological functions, it is important to know the charge state of the residue at the physiological pH and understand the factors that determine them. Therefore, examining how the conformational responses to pH or the ionization of the internal ionizable residues take place is of great importance.

Internal ionizable amino acid residues are relatively rare, but they are essential for many fundamental biological processes in all living organisms. For example, they are necessary for catalysis,1 proton transport,2,3 ion homeostasis,4 light-activated processes,5 and most forms of biological energy transduction.2,6 During mechanistic cycling, these internal ionizable residues can experience different microenvironments, and their protonation states and pKa values would adjust accordingly.7 To understand the functional roles played by internal ionizable residues, it is necessary to know their pKa values and to understand the molecular determinants of these values. The internal ionizable residues often titrate with anomalous pKa values because charged species are usually not compatible with the hydrophobic interior of proteins.8,9 Here for the ionizable residues, we refer primarily to Lys, Glu, Asp, and Arg residues. The internal ionizable residues, which are buried inside a less polar or polarizable environment, prefer to stay at the neutral state, and their pKa values tend to be shifted from the normal values in water, where the charged state of the ionizable residues is predominant.7,10,11 Previous experiments © XXXX American Chemical Society

Received: August 11, 2017 Published: January 8, 2018 A

DOI: 10.1021/jacs.7b08569 J. Am. Chem. Soc. XXXX, XXX, XXX−XXX

Article

Journal of the American Chemical Society

indicated the importance of conformational relaxation and desolvation energy for accurate prediction of the pKa values of these internal titratable residues.50 In this work, we carried out the pH replica exchange MD (pH-REMD)51−53 simulations to investigate the ionizable residues in the interior of SNase. pH-REMD is an efficient method which couples the constant pH MD (CpHMD)47 and replica exchange MD (REMD).54 The primary advantage of CpHMD is that it allows proper sampling of the coupled protein conformations and protonation states of the titratable residues along an MD trajectory.55−59 Protonation changes are attempted regularly and accepted or rejected after evaluating the energetic cost of the proposed change based on the titratable residue’s identity, environment, and solvent pH. As a result, the simulation provides an ensemble of properly weighted protonation states and conformations at the pH of choice.60 This is critically important when the pKa values of the titratable residues of a protein are near the solvent pH, since no single protonation state properly describes the ensemble at that pH, or when the pKa is either unknown or very shifted from the standard values.61 REMD is a family of extended ensemble techniques that have been shown to dramatically improve conformational sampling.62−66 In REMD simulations, a series of independent replicas periodically attempt to exchange information, such as temperature54,62 and pH,51 to sample from an expanded ensemble covering multiple states. In pHREMD, each replica is simulated at a different pH using CpHMD, and the solvent pH is swapped between neighboring replicas. This method has been shown to enhance conformational sampling and accelerate convergence with better pKa prediction precision.52,53 Twenty-one variants of SNase with internal Lys, Glu, and Asp residues were selected in this study. The pKa values of these internal ionizable residues were calculated by performing pH-REMD simulations. Subsequently, the simulated pKa values were compared to experimental results. The general trends in the response of proteins to internal ionizable residues were unraveled at the atomistic level. We also investigated the behaviors of the internal ionizable residues in both charged and neutral forms, and the pH-dependent protein conformational changes. The relation between the protonation of internal ionizable residue and conformation equilibria was discussed in the context of a thermodynamic model. Our study provides valuable insights into pH-driven conformational reorganization coupled to the ionization of internal ionizable residues in the hydrophobic interior of the protein.

To examine determinants of the unusual properties of internal ionizable residues in proteins, Garcia-Moreno and coworkers have systematically studied the internal ionizable residues and measured their pKa values with a family of variants of staphylococcal nuclease (SNase), in which 25 internal positions were replaced with Lys, Glu, or Asp, one at a time.7,10,24−26 The conformational changes due to the ionization of internal groups have been examined with calorimetry, NMR, Trp-fluorescence, and CD spectroscopy experiments.10,12,13,23,27,28 They found that the pKa values of the majority of the internal ionizable residues in SNase are very different from their normal pKa values in water, and some are shifted by as many as 5 pKa units. These shifts are in the direction that favors the neutral form of the ionizable residues, which shows the interior of SNase is not a good solvent for charges.11 SNase is a protein with high intrinsic stability.24 In Garcia-Moreno and co-workers’s experiments, they used a hyperstable variant of SNase Δ+PHS, and in only a few of the mutated cases the internal ionizable residues unfold the protein globally.7,10 For most of cases, the internal ionizable residues have detectable effects on the conformation of the protein,7,10,23,24 and the stability of the variant is highly dependent on pH, indicating that the pKa values of the internal ionizable residues are shifted from their normal values in water.7 Despite the wealth of discoveries obtained from the experiments, the detailed mechanisms of protein conformational responses to the presence of internal ionizable residues remain poorly understood. The theoretical methods applied to examine the structural basis of conformational changes coupled to the ionization of internal residues at the microscopic level would be helpful to rationalize these experimental findings and understand the unusual properties of the internal ionizable residues. The pKa-cooperative in 2009, organized to advance development of accurate and useful computational methods for structure-based calculation of pKa values, challenged computational laboratories to make blind predictions on pKa values experimentally determined by Garcia-Moreno’s group, and many theoretical researchers explored these pKa predictions using various approaches without knowing the answers.29−45 Garcia-Moreno’s group previously performed multiple molecular dynamics (MD) simulations to study the patterns of conformational reorganization triggered by the ionization of internal residues.11 Recently, Zheng and Cui carried out microsecond MD simulations on several mutants of SNase with explicit solvation and provided insights into the spatial and temporal scales of protein and solvent motions that dictate the diverse titration behaviors of buried protein residues.46 In their simulations, those internal ionizable residues were at constant protonation states (protonated or deprotonated), and there were no protonation changes during their simulations. Hence, those simulations could not fully capture the effect of pH on the protonation state, and how the protonation state influences the protein conformation and dynamics. Williams and coworkers once utilized the constant pH MD (CpHMD)47 method for the blind prediction of pKa values of titratable residues in wild type and mutated structures of the SNase protein, and showed the necessity to accurately and sufficiently sample conformational space in order to obtain pKa values consistent with experimental results.34 Shen and co-workers also performed continuous constant pH molecular dynamics (CPHMD)48,49 simulations to study the titratable residues in various hydrophobic regions of SNase and variants, and



RESULTS AND DISCUSSION Titrations of the Internal Ionizable Residues. The structure of SNase and the locations of investigated internal ionizable residues are shown in Figure 1. The 21 variants and pKa values are shown in Table 1. In general, our calculated pKa values from pH-REMD simulations are in good agreement with the experimental values. Most notably, the pKa of I92K is decreased by more than 5 units and gives the largest shift in all 7 cases, which is consistent with the experimental observation. As shown in Figure 1, Lys92 is buried more deeply inside the protein than other Lys residues, and in such a hydrophobic environment the equilibrium between the charged and the neutral form of Lys is shifted in favor of the neutral form, which accounts for the largest decrease of the pKa value. Lys118 is located in a flexible loop, where it can be exposed to a water solution more easily (see Figure 1). Therefore, the calculated B

DOI: 10.1021/jacs.7b08569 J. Am. Chem. Soc. XXXX, XXX, XXX−XXX

Article

Journal of the American Chemical Society

Glu at other positions is more preferred. This can be explained by the local environment of those Glu residues being buried inside the protein, or being less exposed to bulk water, as shown in Figure 1. The pKa shifts of these Glu residues show a similar pattern to those for Lys residues, showing a correlation between burial and pKa shift. The titration curves of those Glu residues are given in Figure 2b. The correlation coefficient (R) between the calculated pKa values of those internal Glu residues and experimental values is 0.93, and the slope is also close to 1.0 (Figure 3b). The experimental pKa values for the internal ionizable Asp residues are not available. Because Asp has the same ionizable group as Glu, the internal Asp residues would show similar behavior as the internal Glu residues to a certain extent.23 Karp et al. once experimentally measured the pKa values of Asp and Glu buried at position 66 in SNase and obtained similar pKa values (8.95 vs 9.07).23 The titration curves of the Asp residues are shown in Figure 2c. According to our predictions, the pKa values of these calculated Asp residues are similar to those of the Glu residues at the same positions in SNase, but not the same. Because the interactions of the side chain terminus with the backbone and surrounding environments would be different, their ability to be buried or exposed to the protein surface would be different, which would give rise to the difference between the pKa values of Asp and Glu residues at the same position. Conformational Changes due to Ionization of Internal Residues. The ionization of internal residues usually triggers protein conformational changes to tolerate the ionization. This means that one should use extreme care when interpreting the pKa shift data just in terms of solvent exposure, since that itself is probably a function of pH. To examine the possible conformational changes coupled to the ionization of internal residues, we first calculated the solvent accessible surface area (SASA) of the ionizable group of each investigated internal residue to determine the pH dependence of the burial. SASA was determined with a probe radius of 1.4 Å, and it was defined in the unit of Å2. The deep buried internal ionizable residues (Lys and Glu) at site 92 of SNase have the largest pKa shifts. In contrast, the site 118 is very close to the surface of the protein, where the internal ionizable residues usually have normal pKa values. The difference in SASAs between these two sites is the most conspicuous, and this comparison would better help to understand the large pKa shift at site 92. So here we mainly show the SASAs of the ionizable groups at sites 92 and 118 in the lowest and highest pH values, which correspond to two completely different protonation states. The SASAs for other ionizable residues studied in this work are shown in Figure S1 of the Supporting Information. The SASA histograms for Lys92, Lys118, Glu92, and Glu118 at the lowest and highest pH values used during the pH-remd runs are depicted in Figure 4. The ionizable groups of all internal residues show increases in SASA, when the ionizable groups are charged (at low pH for Lys, high pH for Glu and Asp), indicating that the ionizable groups are substantially exposed to water. In contrast, when these ionizable groups are in their neutral forms (at high pH for Lys, low pH for Glu and Asp), the SASAs are much lower, which denotes that they are mainly buried inside the protein and inaccessible to water. For the convenience of the following analysis, we define the mode SASA as the value that appears most often in the distribution of the SASA values derived from the MD simulation in a specific pH.

Figure 1. Locations of the internal ionizable residues in the structure of SNase.

Table 1. Calculated and Experimental pKa values of the 21 Mutants of SNase Variants I92K L25K L36K V74K N100 K L38K N118 K N118E L37E L38E F34E V74E L36E I92E N118D L37D L38D F34D V74D L36D I92D

Calculated pKa

Experimental pKa7,10

± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

5.3 6.3 7.2 7.4 8.6 10.4 10.4 4.5 5.2 6.8 7.3 7.8 8.7 9.0 a a a a a a a

4.4 4.9 6.0 7.0 7.3 9.1 9.8 4.4 6.2 6.9 8.1 8.9 8.5 8.6 4.1 6.4 8.1 9.3 7.6 9.3 9.6

0.18 0.17 0.26 0.17 0.13 0.18 0.10 0.08 0.11 0.19 0.09 0.25 0.13 0.13 0.20 0.28 0.24 0.21 0.14 0.11 0.25

a

The experimental pKa values for Asp-containing variants are not available.

pKa of Lys118 is closest to the normal value in water. Lys38 is also built in a loop; however, it is less accessible to water than Lys118. During the REMD simulation, the local flexibility of the loop region also allows Lys38 to be exposed to water. Our calculated result is consistent with Garcia-Moreno and coworkers’ previous studies.28,67 The simulated titration curves for these Lys residues are given in Figure 2a. The correlation coefficient (R) between the calculated pKa values for internal Lys residues and experimental values is 0.98, and the slope is very close to 1.0 (Figure 3a), which demonstrates that our pHREMD simulation is able to accurately reproduce experimental pKa shifts. For the Glu residues, most of the calculated pKa values are much higher than the normal pKa of 4.5 for Glu in water, except for N118E (see Table 1). The upward shifts in pKa for the other internal Glu residues indicate that the neutral form of C

DOI: 10.1021/jacs.7b08569 J. Am. Chem. Soc. XXXX, XXX, XXX−XXX

Article

Journal of the American Chemical Society

Figure 2. Titration curves for the (a) Lys-containing variants, (b) Glu-containing variants, and (c) Asp-containing variants.

Figure 3. Correlation of the experimental and calculated pKa values for the (a) Lys-containing variants and (b) Glu-containing variants. The normal pKa values for Lys and Glu residues in solution are 10.4 and 4.5 respectively.

each case, that the amount of solvent exposure is highly pH dependent, and that the pH corresponding to the midpoint of these distributions agrees with the pKa values obtained if we look only at protonation versus pH. Conformational changes coupled to the ionization of an internal ionizable residue should be reflected in pKa changes during the dynamics. To see this effect, we can compute the changes in pKa during sampling, at a pH close to the macroscopic pKa. For the ionizable residue deeply buried inside with little chance to access the solution or conversely to be fully exposed to the solution, it usually prefers to stay at a fixed adaptive protonation state during the dynamics, and its pKa does not change much accordingly. However, for V74K and L38E, their observed pKa values are near the midpoint of the pKa ranges of the studied Lys and Glu residues located at different sites. And this means that Lys74/Glu38 could be half buried inside the protein and half exposed to the solution during the dynamics, reflecting the large amplitude of pKa and conformational changes. So we select these two variants to

For Lys92, which has the largest pKa shift among the seven investigated Lys residues, the majority of snapshots in MD simulations have SASA values (the mode SASA) around 26.0 Å2 at the lowest pH (3.0), where Lys92 is mostly charged. However, for Lys118, which has a normal pKa, most of snapshots in the simulations have SASA values around 68.5 Å2 when the Lys118 is charged at pH = 3.0. In the highest pH (12.0), where the internal Lys residues stay in their neutral form, the SASAs of the ionizable groups for the majority of snapshots in MD simulations are very low, meaning they are mostly buried. The smaller increase in SASA for Lys92 from neutral to charged state indicates that it is buried more deeply than Lys118, which accounts for the larger pKa shift for Lys92. The pair of Glu92 and Glu118 shows a similar pattern to that of Lys92 and Lys118. At the highest pH (11.0), where the Glu residues are charged, the peak appears when SASA is around 17.5 Å2 for Glu92, and it is around 107.5 Å2 for Glu118. We can follow the change in solvent exposure versus pH by integrating over the histograms in Figure 4. Figure 5 shows, in D

DOI: 10.1021/jacs.7b08569 J. Am. Chem. Soc. XXXX, XXX, XXX−XXX

Article

Journal of the American Chemical Society

Figure 4. SASA of the ionizable group in (a) Lys92, (b) Lys118, (c) Glu92, and (d) Glu118.

in which the internal ionizable groups were unambiguously buried or exposed to water. For V74K and L38E, we discriminated the internal ionizable groups using the SASA criterion: when the SASA of the ionizable group is larger than the specific mode SASA (for Lys, it corresponds to the mode SASA in the lowest pH, and for Glu, it corresponds to the mode SASA in the highest pH), we regard it as exposed to water; when the SASA of the ionizable group is smaller than the other specific mode SASA (for Lys, it corresponds to the mode SASA in the highest pH, and for Glu, it corresponds to the mode SASA in the lowest pH), we regard it as buried inside the protein. The resulting titration curves given in Figure 8 were also fitted to the Hill equation. For Lys74, the pKa for the conformations with the ionizable group buried inside the protein is 5.5, while the pKa for the conformations with the ionizable group exposed to water is 9.3. For Glu38, the pKa values are 8.9 and 5.1, respectively, when the internal ionizable group is buried inside the protein or exposed to water. As expected, based on the local chemical environment of these internal ionizable residues, when the internal ionizable residues are buried deeply inside the protein, their pKa values are considerably shifted in the direction that promotes the neutral state. When the internal ionizable residues are accessible to water, their pKa values are close to the normal values in the water solution. Our result demonstrates that the internal ionizable residues with notable shifted pKa values usually have conformation-specific pKa values, which are important for protein functions via pH-dependent allostery.68,69 Previous studies have demonstrated that the observed pKa is not simply an average of the conformation-specific pKa values because the conformational equilibrium is itself pH-dependent

calculate the chunk pKa and investigate the conformational changes in the following analysis. We computed pKa values over sections of 100 000 MD steps (chunk pKa) and show them versus time in Figure 6 for V74K and L38E. The observed pKa for Lys74 is 7.0, but the chunk pKa fluctuation ranges from 5.5 to 8.3. On the other hand, the chunk pKa fluctuation ranges from 6.0 to 8.5 for Glu38, whose observed pKa is 6.9. As shown in Figure 6, the chunk pKa fluctuates significantly during the REMD simulations, indicating that the protein undergoes large local conformational changes or even global changes due to the ionization of the internal residues. Snapshots, in the replica under pH = 7 (the closest pH to their pKa values), were extracted from the pH-REMD simulations when the internal ionizable residues were neutral and ionized, respectively (see Figure 7). The protein conformational changes, especially the local ones, caused by the ionization of the internal residues are remarkable. For V74K, the loop connected to Lys74 is in a halfopen state when Lys74 is in its neutral form, compared to the more largely open state when Lys74 is ionized. L38E also shows a similar pattern. The neutral Glu38 is completely covered by a loop which is on top of it, and this loop remarkably opens up when Glu38 becomes ionized. And this loop conformational change also induces the reorientation of the connected helical structure, as shown in Figure 7b. The ionized Lys74 and Glu38 induced opening of neighbor crevices that connected the ionizable groups to bulk water or the protein−water interface. This is probably a common regulatory mechanism for the internal ionizable residues existing in proteins. In order to obtain the conformation-specific pKa values, we titrated the internal ionizable residues based on the structures E

DOI: 10.1021/jacs.7b08569 J. Am. Chem. Soc. XXXX, XXX, XXX−XXX

Article

Journal of the American Chemical Society

Figure 5. Number of frames which have mode SASA in varying pH for (a) Lys92, (b) Lys118, (c) Glu92, and (d) Glu118. The mode SASA is the approximate value that appears most often in the distribution of the SASA values derived from the MD simulation in each pH.

Figure 6. Chunk pKa for (a) Lys74 and (b) Glu38.

and must be taken into account.16,60 This can be accomplished through a four-state thermodynamic model, as summarized in Figure 9. This model has four parameters: the conformationspecific pKa values when the investigated ionizable residue is exposed to water or buried inside the protein, and the conformational equilibrium constants when the investigated ionizable residue is in the neutral (K0) or ionized state (K*). Only three of these four parameters are independent, and their values are constrained. Using the conformation-specific pKa values obtained from pH-REMD simulations, we calculated the conformational equilibrium constants using the following equation,60

pK a(obs) = pK a(buried) + pK a(exposed) K* + 1 − log K *K a(buried) + K a(exposed)

(1)

The conformational equilibrium constants (K*), when the internal ionizable residues are ionized, are 4.4 × 10−3 and 98.1, respectively, for V74K and L38E. This model is applicable to any pH-dependent conformational changes and regulation by a single amino acid.



CONCLUSIONS In this work, 21 variants of SNase with Lys, Glu, or Asp residues located in the interior of the protein were studied F

DOI: 10.1021/jacs.7b08569 J. Am. Chem. Soc. XXXX, XXX, XXX−XXX

Article

Journal of the American Chemical Society

Figure 7. Conformational changes due to the ionization of the internal residues for (a) V74K and (b) L38E. The structure with the ionized internal residue is shown in cyan, and the structure with the neutral internal residue is shown in pink. The internal ionizable residues are shown in sticks, with the C atoms of the ionized residue in cyan color and the C atoms of the neutral residue in pink color. The O, N, and H atoms are labeled in red, blue, and white colors, respectively.

Figure 8. Titration curves of the ensembles with the internal ionizable groups buried (red), exposed to water (blue), and for the complete ensemble (black), respectively, for (a) Lys74 and (b) Glu38.

direction where the neutral state was favored. The calculated pKa shifts are in excellent agreement with the experimental observation. Protein conformational changes triggered by the ionization of the internal ionizable residues were also examined. SASA values of the internal ionizable groups showed that all internal ionizable residues were buried when they were at the neutral states. The ionization of the internal residue leads to the increase in solvent accessibility through opening small crevices that connected the ionizable groups to bulk water or the protein−water interface. pH-REMD simulations help to explore the protein conformational changes, which are linked to the protonation of internal ionizable residues. The ionizable residues, which have large shifted pKa values, usually have two significantly different conformation-specific pKa values. On the basis of our results, we presented a general thermodynamic model showing protein conformational changes regulated by a single amino acid. The four-state model based on a single ionizable residue and two protein conformations allows us to analyze the basic requirement for efficient conformational switching, driven by a single amino acid that is solvent exposed in one of the conformations. The results from this study help to better understand the properties of internal ionizable residues in proteins. This allows the study of biological processes where internal ionizable groups are fundamental structural motifs necessary for protein

Figure 9. Thermodynamic cycle. pKa (buried) is obtained using the conformations with the internal ionizable residue buried inside the protein. pKa (exposed) is obtained using the conformations with the internal ionizable residue exposed to water or protein−water interface. Res0 and Res* denote the neutral and ionized form of the internal ionizable residue, respectively, and K0 and K* are the corresponding equilibrium constants.

through pH-REMD simulations, aiming to better understand the unusual properties of the internal ionizable residues. Our result demonstrates that, in most of the cases we investigated, the pKa values of these internal ionizable residues were shifted from their normal values in water, and the shifts were in the G

DOI: 10.1021/jacs.7b08569 J. Am. Chem. Soc. XXXX, XXX, XXX−XXX

Article

Journal of the American Chemical Society

deviating from 1 in the case of interacting titratable residues because of cooperativity.79 The fits were performed with Grace using the nonlinear fitting method.

function. It suggests that some degree of conformational change may be expected to facilitate the charge transfer reactions, which has an important influence on the kinetic and thermodynamic properties of proton transport mechanisms.





S Supporting Information *

METHODS

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/jacs.7b08569. SASAs of the ionizable groups in the internal residues (PDF)

System Setup. The crystal structure of SNase was downloaded from the Protein Data Bank (PDB id: 3BDC). The ligands, water molecules, and ions were removed from the PDB file, and only the protein was kept and used as the template. For each of the 21 starting structures used in our simulations, one specific internal residue in the template was mutated to Lys, Glu, or Asp, one at a time. In GarciaMoreno and co-workers’ experiments, they engineered variants of SNase with Lys, Glu, and Asp, respectively, at 25 different internal positions and then measured their corresponding pKa values. In our study, 7 Lys-containing, 7 Glu-containing, and 7 Asp-containing systems, which have notable differences in their measured pKa values, were constructed. Topology files were built using the tleap module in AmberTools 14.70 The Amber ff99SB force field71 was used for the protein. The SHAKE algorithm72 was employed to keep bonds involving hydrogen atoms at their equilibrium lengths. We utilized the Generalized Born (GB) implicit solvent model73 (with igb = 8) to mimic the solvent environment, and the mbondi3 Born radii74 were used in GB simulation. All simulations adopted different random seeds to avoid synchronization artifacts.75 Structure Equilibration. Before the pH-REMD simulations, an equilibration procedure was employed for each starting system. The protein was first minimized with a positional restraint of 10.0 kcal mol−1 Å−2 placed on the backbone atoms. Then the system was heated gradually from 10.0 to 300.0 K over a 1.0 ns MD simulation, where a 1.0 kcal mol−1 Å−2 restraint was placed on the backbone atoms. After heating, each system was equilibrated for 1.0 ns with the backbone restraint reduced to 0.1 kcal mol−1 Å−2. During the equilibration procedure, no protonation state changes were attempted. pH-REMD Simulations. All pH-REMD simulations were run with 10 equally spaced pH replicas. For the Lys-containing systems, the pH ranged from 3 to 12 with an interval of 1, while, for the Glu- and Aspcontaining systems, the pH ranged from 2 to 11 with an interval of 1. The pH ranges were sufficiently large to describe the behavior of the investigated titratable residues. We used a 0.1 M salt concentration and an infinite cutoff for nonbonded interactions. Langevin dynamics76 was utilized in order to maintain the temperature at 300 K with a collision frequency of 5.0 ps−1. We attempted a change in protonation state every 10 steps. In each system, only the investigated internal ionizable residue was allowed to change its protonation state. The atomic partial charges for the neutral form of these investigated internal ionizable residues were taken from the corresponding residue definitions in the Amber ff99SB force field. For the other titratable residues, the lysine, arginine, glutamic, and aspartic amino acid residues were modeled in their ionized form; the cysteine, tyrosine, and histidine (with the hydrogen placed at the ε position) residues were modeled in their neutral form. Each replica was initially simulated at its corresponding pH for 1 ns, using the pmemd.cuda.MPI module, before any exchange between replicas was allowed. During the 50 ns production run, we attempted to exchange between replicas every 1000 steps. We performed five fully independent pH-REMD simulations to estimate error bars. The trajectories were analyzed using the cpptraj module77 of amber, and VMD.78 The program cphstats in AmberTools 14 was utilized to calculate the resulting protonation state distributions from the pH-REMD simulations. All simulations were performed with the AMBER 14 program.70 The pKa of each investigated residue was obtained by fitting the titration data at each pH to the Hill equation shown below,

fd =



AUTHOR INFORMATION

Corresponding Authors

*[email protected] (X.H.) *roitberg@ufl.edu (A.E.R.) ORCID

Jinfeng Liu: 0000-0002-0614-5097 John Z. H. Zhang: 0000-0003-4612-1863 Xiao He: 0000-0002-4199-8175 Adrian E. Roitberg: 0000-0003-3963-8784 Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This work was supported by the Ministry of Science and Technology of China (Grant No. 2016YFA0501700), National Natural Science Foundation of China (Grant Nos. 21703289, 21761132022, 21673074, 21433004), Young Top-Notch Talent Support Program of Shanghai, Shanghai Putuo District (Grant 2014-A-02), the NYU-ECNU Center for Computational Chemistry at NYU Shanghai, and NYU Global Seed Grant for Collaborative Research. This research is part of the Blue Waters sustained-petascale computing project, which is supported by the National Science Foundation (Awards OCI0725070 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. This work is also part of the PRAC allocation support by the National Science Foundation (Award Number ACI 1515572). We also thank the High-Performance Computing Center of China Pharmaceutical University and Supercomputer Center of East China Normal University for providing us with computational time.



REFERENCES

(1) Warshel, A.; Sharma, P. K.; Kato, M.; Parson, W. W. Biochim. Biophys. Acta, Proteins Proteomics 2006, 1764, 1647−1676. (2) Luecke, H.; Richter, H. T.; Lanyi, J. K. Science 1998, 280, 1934− 1937. (3) Yoshikawa, S.; Shinzawa-Itoh, K.; Nakashima, R.; Yaono, R.; Yamashita, E.; Inoue, N.; Yao, M.; Fei, M. J.; Libeu, C. P.; Mizushima, T.; Yamaguchi, H.; Tomizaki, T.; Tsukihara, T. Science 1998, 280, 1723−1729. (4) Doyle, D. A.; Cabral, J. M.; Pfuetzner, R. A.; Kuo, A. L.; Gulbis, J. M.; Cohen, S. L.; Chait, B. T.; MacKinnon, R. Science 1998, 280, 69− 77. (5) Luecke, H.; Lanyi, J. K. Adv. Protein Chem. 2003, 63, 111−130. (6) Le, N. P.; Omote, H.; Wada, Y.; Al-Shawi, M. K.; Nakamoto, R. K.; Futai, M. Biochemistry 2000, 39, 2778−2783. (7) Isom, D. G.; Castaneda, C. A.; Cannon, B. R.; Garcia-Moreno, B. E. Proc. Natl. Acad. Sci. U. S. A. 2011, 108, 5260−5265. (8) Robinson, A. C.; Majumdar, A.; Schlessman, J. L.; GarciaMoreno, B. Biochemistry 2017, 56, 212−218.

1 1 + 10n(pKa − pH)

ASSOCIATED CONTENT

(2)

where fd is the fraction of the total simulation that the investigated residue spent in the deprotonated state and n is the Hill coefficient, which should approach 1 for noninteracting ionizable residues, H

DOI: 10.1021/jacs.7b08569 J. Am. Chem. Soc. XXXX, XXX, XXX−XXX

Article

Journal of the American Chemical Society (9) Harris, T. K.; Turner, G. J. IUBMB Life 2002, 53, 85−98. (10) Isom, D. G.; Castaneda, C. A.; Velu, P. D.; Garcia-Moreno, B. Proc. Natl. Acad. Sci. U. S. A. 2010, 107, 16096−16100. (11) Damjanovic, A.; Brooks, B. R.; Garcia-Moreno, E. B. J. Phys. Chem. A 2011, 115, 4042−4053. (12) Chimenti, M. S.; Khangulov, V. S.; Robinson, A. C.; Heroux, A.; Majumdar, A.; Schlessman, J. L.; Garcia-Moreno, B. Structure 2012, 20, 1071−1085. (13) Richman, D. E.; Majumdar, A.; Garcia-Moreno, E. B. Biochemistry 2015, 54, 5888−5897. (14) Vorburger, T.; Ebneter, J. Z.; Wiedenmann, A.; Morger, D.; Weber, G.; Diederichs, K.; Dimroth, P.; von Ballmoos, C. FEBS J. 2008, 275, 2137−2150. (15) Fillingame, R. H.; Angevine, C. M.; Dmitriev, O. Y. FEBS Lett. 2003, 555, 29−34. (16) Di Russo, N. V.; Estrin, D. A.; Marti, M. A.; Roitberg, A. E. PLoS Comput. Biol. 2012, 8, e1002761. (17) Laricheva, E. N.; Goh, G. B.; Dickson, A.; Brooks, C. L. J. Am. Chem. Soc. 2015, 137, 2892−2900. (18) Huang, Y. D.; Chen, W.; Dotson, D. L.; Beckstein, O.; Shen, J. Nat. Commun. 2016, 7, 12940. (19) Lang, E. J. M.; Heyes, L. C.; Jameson, G. B.; Parker, E. J. J. Am. Chem. Soc. 2016, 138, 2036−2045. (20) Gayen, A.; Leninger, M.; Traaseth, N. J. Nat. Chem. Biol. 2016, 12, 141−145. (21) Denisov, V. P.; Schlessman, J. L.; Garcia-Moreno, B.; Halle, B. Biophys. J. 2004, 87, 3982−3994. (22) Karp, D. A.; Gittis, A. G.; Stahley, M. R.; Fitch, C. A.; Stites, W. E.; Garcia-Moreno, B. Biophys. J. 2007, 92, 2041−2053. (23) Karp, D. A.; Stahley, M. R.; Garcia-Moreno, E. B. Biochemistry 2010, 49, 4138−4146. (24) Isom, D. G.; Cannon, B. R.; Castaneda, C. A.; Robinson, A.; Bertrand, G. M. E. Proc. Natl. Acad. Sci. U. S. A. 2008, 105, 17784− 17788. (25) Robinson, A. C.; Castaneda, C. A.; Schlessman, L. J.; GarciaMoreno, E. B. Proc. Natl. Acad. Sci. U. S. A. 2014, 111, 11685−11690. (26) Fitch, C. A.; Karp, D. A.; Lee, K. K.; Stites, W. E.; Lattman, E. E.; Garcia-Moreno, B. Biophys. J. 2002, 82, 3289−3304. (27) Chimenti, M. S.; Castaneda, C. A.; Majumdar, A.; GarciaMoreno, B. J. Mol. Biol. 2011, 405, 361−377. (28) Harms, M. J.; Castaneda, C. A.; Schlessman, J. L.; Sue, G. R.; Isom, D. G.; Cannon, B. R.; Garcia-Moreno, B. J. Mol. Biol. 2009, 389, 34−47. (29) Alexov, E.; Mehler, E. L.; Baker, N.; Baptista, A. M.; Huang, Y.; Milletti, F.; Nielsen, J. E.; Farrell, D.; Carstensen, T.; Olsson, M. H. M.; Shen, J. K.; Warwicker, J.; Williams, S.; Word, J. M. Proteins: Struct., Funct., Genet. 2011, 79, 3260−3275. (30) Nielsen, J. E.; Gunner, M. R.; Garcia-Moreno, E. B. Proteins: Struct., Funct., Genet. 2011, 79, 3249−3259. (31) Witham, S.; Talley, K.; Wang, L.; Zhang, Z.; Sarkar, S.; Gao, D. Q.; Yang, W.; Alexov, E. Proteins: Struct., Funct., Genet. 2011, 79, 3389−3399. (32) Wang, L.; Li, L.; Alexov, E. Proteins: Struct., Funct., Genet. 2015, 83, 2186−2197. (33) Carstensen, T.; Farrell, D.; Huang, Y.; Baker, N. A.; Nielsen, J. E. Proteins: Struct., Funct., Genet. 2011, 79, 3287−3298. (34) Williams, S. L.; Blachly, P. G.; McCammon, J. A. Proteins: Struct., Funct., Genet. 2011, 79, 3381−3388. (35) Arthur, E. J.; Yesselman, J. D.; Brooks, C. L. Proteins: Struct., Funct., Genet. 2011, 79, 3276−3286. (36) Gunner, M. R.; Zhu, X. Y.; Klein, M. C. Proteins: Struct., Funct., Genet. 2011, 79, 3306−3319. (37) Czodrowski, P. Proteins: Struct., Funct., Genet. 2011, 79, 3299− 3305. (38) Warwicker, J. Proteins: Struct., Funct., Genet. 2011, 79, 3374− 3380. (39) Kitahara, R.; Hata, K.; Maeno, A.; Akasaka, K.; Chimenti, M. S.; Garcia-Moreno, B.; Schroer, M. A.; Jeworrek, C.; Tolan, M.; Winter,

R.; Roche, J.; Roumestand, C.; de Guillen, K. M.; Royer, C. A. Proteins: Struct., Funct., Genet. 2011, 79, 1293−1305. (40) Word, J. M.; Nicholls, A. Proteins: Struct., Funct., Genet. 2011, 79, 3400−3409. (41) Meyer, T.; Kieseritzky, G.; Knapp, E. W. Proteins: Struct., Funct., Genet. 2011, 79, 3320−3332. (42) Shan, J. F.; Mehler, E. L. Proteins: Struct., Funct., Genet. 2011, 79, 3346−3355. (43) Burger, S. K.; Ayers, P. W. Proteins: Struct., Funct., Genet. 2011, 79, 2044−2052. (44) Webb, H.; Tynan-Connolly, B. M.; Lee, G. M.; Farrell, D.; O’Meara, F.; Sondergaard, C. R.; Teilum, K.; Hewage, C.; McIntosh, L. P.; Nielsen, J. E. Proteins: Struct., Funct., Genet. 2011, 79, 685−702. (45) Olsson, M. H. M. Proteins: Struct., Funct., Genet. 2011, 79, 3333−3345. (46) Zheng, Y. Q.; Cui, Q. Proteins: Struct., Funct., Genet. 2017, 85, 268−281. (47) Mongan, J.; Case, D. A.; McCammon, J. A. J. Comput. Chem. 2004, 25, 2038−2048. (48) Wallace, J. A.; Shen, J. K. J. Chem. Theory Comput. 2011, 7, 2617−2629. (49) Goh, G. B.; Hulbert, B. S.; Zhou, H. Q.; Brooks, C. L. Proteins: Struct., Funct., Genet. 2014, 82, 1319−1331. (50) Wallace, J. A.; Wang, Y. H.; Shi, C. Y.; Pastoor, K. J.; Nguyen, B. L.; Xia, K.; Shen, J. K. Proteins: Struct., Funct., Genet. 2011, 79, 3364− 3373. (51) Itoh, S. G.; Damjanovic, A.; Brooks, B. R. Proteins: Struct., Funct., Genet. 2011, 79, 3420−3436. (52) Dashti, D. S.; Meng, Y. L.; Roitberg, A. E. J. Phys. Chem. B 2012, 116, 8805−8811. (53) Swails, J. M.; Roitberg, A. E. J. Chem. Theory Comput. 2012, 8, 4393−4404. (54) Sugita, Y.; Okamoto, Y. Chem. Phys. Lett. 1999, 314, 141−151. (55) Baptista, A. M.; Teixeira, V. H.; Soares, C. M. J. Chem. Phys. 2002, 117, 4184−4200. (56) Baptista, A. M.; Martel, P. J.; Petersen, S. B. Proteins: Struct., Funct., Genet. 1997, 27, 523−544. (57) Burgi, R.; Kollman, P. A.; van Gunsteren, W. F. Proteins: Struct., Funct., Genet. 2002, 47, 469−480. (58) Borjesson, U.; Hunenberger, P. H. J. Phys. Chem. B 2004, 108, 13551−13559. (59) Khandogin, J.; Brooks, C. L. Biophys. J. 2005, 89, 141−157. (60) Di Russo, N. V.; Marti, M. A.; Roitberg, A. E. J. Phys. Chem. B 2014, 118, 12818−12826. (61) Mcgee, T. D.; Edwards, J.; Roitberg, A. E. J. Phys. Chem. B 2014, 118, 12577−12585. (62) Pitera, J. W.; Swope, W. Proc. Natl. Acad. Sci. U. S. A. 2003, 100, 7587−7592. (63) Chodera, J. D.; Shirts, M. R. J. Chem. Phys. 2011, 135, 194110. (64) Nadler, W.; Meinke, J. H.; Hansmann, U. H. E. Phys. Rev. E 2008, 78, 061905. (65) Meng, Y. L.; Roitberg, A. E. J. Chem. Theory Comput. 2010, 6, 1401−1412. (66) Meng, Y. L.; Dashti, D. S.; Roitberg, A. E. J. Chem. Theory Comput. 2011, 7, 2721−2727. (67) Harms, M. J.; Schlessman, J. L.; Chimenti, M. S.; Sue, G. R.; Damjanovic, A.; Garcia-Moreno, B. Protein Sci. 2008, 17, 833−845. (68) Tanford, C. J. Am. Chem. Soc. 1961, 83, 1628−1634. (69) Tanford, C.; Taggart, V. G. J. Am. Chem. Soc. 1961, 83, 1634− 1638. (70) Case, D. A.; Cheatham, T. E.; Darden, T.; Gohlke, H.; Luo, R.; Merz, K. M.; Onufriev, A.; Simmerling, C.; Wang, B.; Woods, R. J. J. Comput. Chem. 2005, 26, 1668−1688. (71) Hornak, V.; Abel, R.; Okur, A.; Strockbine, B.; Roitberg, A.; Simmerling, C. Proteins: Struct., Funct., Genet. 2006, 65, 712−725. (72) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327−341. (73) Onufriev, A.; Bashford, D.; Case, D. A. Proteins: Struct., Funct., Genet. 2004, 55, 383−394. I

DOI: 10.1021/jacs.7b08569 J. Am. Chem. Soc. XXXX, XXX, XXX−XXX

Article

Journal of the American Chemical Society (74) Nguyen, H.; Roe, D. R.; Simmerling, C. J. Chem. Theory Comput. 2013, 9, 2020−2034. (75) Sindhikara, D. J.; Kim, S.; Voter, A. F.; Roitberg, A. E. J. Chem. Theory Comput. 2009, 5, 1624−1631. (76) Pastor, R. W.; Brooks, B. R.; Szabo, A. Mol. Phys. 1988, 65, 1409−1419. (77) Roe, D. R.; Cheatham, T. E. J. Chem. Theory Comput. 2013, 9, 3084−3095. (78) Humphrey, W.; Dalke, A.; Schulten, K. J. Mol. Graphics 1996, 14, 33−38. (79) Weiss, J. N. FASEB J. 1997, 11, 835−841.

J

DOI: 10.1021/jacs.7b08569 J. Am. Chem. Soc. XXXX, XXX, XXX−XXX