This Article Abstract Full Text (PDF) Supplemental File All Versions of this Article: biophysj.105.061341v1 89/1/141    most recent Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Khandogin, J. Articles by Brooks, C. L. Search for Related Content PubMed PubMed Citation Articles by Khandogin, J. Articles by Brooks, C. L., III

Constant pH Molecular Dynamics with Proton Tautomerism

Department of Molecular Biology, The Scripps Research Institute, La Jolla, California

Correspondence: Address reprint requests to Jana Khandogin, Tel.: 858-784-8070; Fax: 858-784-8688; E-mail: janakhan{at}scripps.edu.

	ABSTRACT

TOP ABSTRACT INTRODUCTION THEORY METHODS RESULTS AND DISCUSSION CONCLUDING REMARKS SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
The current article describes a new two-dimensional -dynamics method to include proton tautomerism in continuous constant pH molecular dynamics (CPHMD) simulations. The two-dimensional -dynamics framework is used to devise a tautomeric state titration model for the CPHMD simulations involving carboxyl and histidine residues. Combined with the GBSW implicit solvent model, the new method is tested on titration simulations of blocked histidine and aspartic acid as well as two benchmark proteins, turkey ovomucoid third domain (OMTKY3) and ribonuclease A (RNase A). A detailed analysis of the errors inherent to the CPHMD methodology as well as those due to the underlying solvation model is given. The average absolute error for the computed pK_a values in OMTKY3 is 1.0 pK unit. In RNase A the average absolute errors for the carboxyl and histidine residues are 1.6 and 0.6 pK units, respectively. In contrast to the previous work, the new model predicts the correct sign for all the pK_a shifts, but one, in the benchmark proteins. The predictions of the tautomeric states of His¹² and His⁴⁸ and the conformational states of His⁴⁸ and His¹¹⁹ are in agreement with experiment. Based on the simulations of OMTKY3 and RNase A, the current work has demonstrated the capability of the CPHMD technique in revealing pH-coupled conformational dynamics of protein side chains.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION THEORY METHODS RESULTS AND DISCUSSION CONCLUDING REMARKS SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
The stability and function of proteins are dependent on the environmental pH (1). Well-known examples of pH-dependent biological processes include fibril formation from one of 16 normally soluble and functional human amyloidogenic proteins such as amyloid ß-peptides, transthyretin, and prion proteins (2,3), membrane insertion of diphtheria toxin (4) and influenza virus hemagglutinin (5), and proton gradient-driven ATP synthesis (6), as well as the catalytic pathway of dihydrofolate reductase (7). In a conventional molecular dynamics (MD) simulation the protonation state of the protein is preset and kept fixed according to the known pK_a values of the ionizable side chains. Two problems arise from the assumption of static protonation. First, because the pK_a value of an ionizable side chain is typically unknown in its protein environment one has to resort to the value of an isolated amino acid (model pK_a), which sometimes deviates significantly from the value in the protein. Second, the protonation/deprotonation of an ionizable side chain is a process that is coupled to the dynamics of the local environment. By artificially fixing the protonation state in a MD simulation, the dynamics of the ionization equilibrium is neglected, thereby precluding detailed understanding of the specific interactions that give rise to the pH-coupled biological phenomena.

The history of pH-coupled molecular dynamics is relatively short. A decade ago, Mertz and Pettitt (8) demonstrated a titration simulation of acetic acid using an open system Hamiltonian. Sham et al. (9) showed that the pK_a values of lysozyme can be accurately modeled using the protein dipoles Langevin dipoles model, which treats the protein relaxation in the microscopic framework of the linear response approximation. In more recent years, several simulation techniques that allow pH-coupled molecular dynamics have emerged. Most of these techniques are based on the combination of molecular dynamics (MD) and Monte Carlo (MC) sampling of protonation states, with the major differences being the solvent representation and the way the free energy change accompanying the switch of the protonation state is computed. In the methods of Baptista et al. (10,11) and Bürgi et al. (12), explicit solvent is used to generate the MD trajectory. While the method of Bürgi et al. (12) utilizes thermodynamic integration to determine the protonation free energy change of a single residue, the method of Baptista et al. (11) makes use of static pK_a calculations based on the implicit solvent Poisson-Boltzmann (PB) model. More recently, these hybrid MD/MC schemes were extended to implicit solvent simulations: Dlugosz and Antosiewicz (13) make use of the analytic continuum electrostatic and the PB models for generating the MD trajectory and evaluating the protonation free energy change respectively; Mongan et al. (14) use the Generalized-Born (GB) model throughout the simulation.

There are several issues that have to be dealt with in the hybrid MD/MC constant pH simulation schemes. First, the periodic abrupt switch in the protonation state introduces a discontinuity in energy and force, which may result in conformational and energetic instabilities. Second, since only one protonation site at an instantaneous conformation of the protein is randomly chosen during a MC step, sampling convergence with regard to protonation and conformational states may be slow for proteins containing many ionizable sites. Third, some issues remain relating to the determination of the free energy change accompanying the protonation state switch. In the methods of Baptista et al. (11) and Dlugosz and Antosiewicz (13), which use the instantaneous protein conformation for the evaluation of the free energy difference, a somewhat ad hoc protein dielectric constant has to be chosen to mediate the electrostatic interactions in the otherwise static protein configuration. On the other hand, in the method of Bürgi et al. (12), which employs a short run (tens of picoseconds) of free energy simulation such as thermodynamic integration, the increased computational cost associated with the large of number of MC steps has to be considered.

In contrast to the discrete methods, the methods of Börjesson and Hünenberger (15), and of Lee et al. (16) are based on simultaneous propagation of continuous titration and conformational degrees of freedom. Both approaches use linear combinations of the deprotonated and protonated states to model nonbonded interactions. Although the method of Börjesson and Hünenberger relies completely upon explicit solvent simulation and uses the -variable to represent the extent of protonation for each ionizable site (15), the method of Lee et al., which has its roots in the -dynamics method (17), utilizes the implicit GB solvent model (18) and attaches physical meaning only to the end points. To minimize the possible distortion due to the unphysical mixed states that are nonetheless necessary for the transition between the end-point states, an artificial titration barrier is utilized in the continuous pH molecular dynamics method (CPHMD) of Lee et al. (16), which serves to prolong the residency time of the end-point states. The combined advantage of a stable dynamics trajectory and low computational cost makes the CPHMD approach of Lee et al. particularly suitable and attractive as a means of exploring pH-dependent conformational processes in proteins.

One prerequisite for the application of the constant pH molecular dynamics method to pH-coupled protein dynamics is reasonable accuracy in the prediction of protonation states or pK_a values for ionizable side chains. Among the aforementioned methods, so far, only three have been tested on proteins. While two of them, the hybrid MD/MC methods of Bürgi et al. (12) and Mongan et al. (14) were tested only on hen egg white lysozyme, the CPHMD method of Lee et al. (16) was tested on four different proteins, including turkey ovomucoid third domain (OMTKY3), bovine pancreatic ribonuclease A (RNase A), bovine trypsin inhibitor, and hen egg white lysozyme. In this pilot study, the total average absolute error over all proteins was 1.6 pK units (16).

A quick survey of the computed pK_a values in the work of Lee et al. (16) reveals that the average absolute error is largest for the histidine groups (3.2 pK units), followed by the carboxyl (1.6 pK units) and amino groups (1.2 pK units). One obvious reason for the increased errors in predicting the pK_a of histidine is the simplistic treatment of proton isomerism, in which case the imidazole ring is split into two titration halves involving the ND1 and NE2 sites with an added penalty potential to prevent double deprotonation (16). This model will be referred to as the split model in the rest of the article. In a constant pH MD simulation, the split model gives rise to two problems. First, since the two titration halves have common atoms, proper atomic charge states cannot be assigned to the neutral histidine tautomers, HSE (protonated on NE2) and HSD (protonated on ND1). Second, the penalty function to eliminate double deprotonated histidine, K₁₂, introduces an undesired biasing force.

In the present work, a new tautomeric state titration model, also referred to as the double-site model, is developed. This model allows simultaneous titration at two competing protonation sites, such as the NE2 and ND1 atoms in histidine. The new model utilizes an additional dynamic variable that is treated on an equal footing to the titration degrees of freedom to control the tautomer interconversion process. Besides histidine, which exhibits the true proton tautomerism in solution, carboxyl groups such as Asp, Glu, and the C-terminus can be considered as having two tautomeric protonation sites, since the exchange of two titrating oxygens due to rotation around the single bond is slow on the routine MD simulation timescale. In essence, our tautomeric state titration model is a two-dimensional -dynamics approach that can, in general, be extended to additional dimensions. For example, the titration of amino groups can be viewed as a tautomeric equilibrium among three states and can be dealt with by constructing a three-dimensional model. However, since most amino groups do not exhibit large pK_a shifts, as seen for example in hen egg white lysozyme, the necessity for the development of a more sophisticated model is questionable. In fact, as mentioned previously, the errors in the computed pK_a values for amino groups are the smallest in the existing CPHMD method of Lee et al. (16) Finally, it is worthwhile noting that although the two-dimensional -dynamics approach is applied to the coupled protonation and tautomeric interconversion equilibria in the present work, the approach is general and can be applied to other problems.

The CPHMD method is further extended in the present study to incorporate the newly developed GBSW implicit solvent model (19) because of the added capability of this method for simulations in a membrane environment (20) and for the greater computational efficiency (21) relative to the GBMV method (18) that was implemented in the original CPHMD code. Since the accuracy of the CPHMD simulations is intimately linked to conformational sampling and the underlying implicit solvent model, the applications of CPHMD will benefit from the improved sampling schemes as well as the solvent model.

The purpose of the present work is to extend and improve the existing CPHMD method developed by Lee et al. (16). In their work (16) it was noted that CPHMD simulations using the GBMV solvent model gave an overestimation of 0.3 pK units for the blocked aspartic acid. This has prompted us to examine the details of the titration simulations of model compounds aimed at understanding the convergence behavior and errors inherent to the CPHMD methodology. Other errors in the prediction of protein ionization equilibrium arising from the force field and implicit solvent model will be explored by revisiting two proteins: OMTKY3 and RNase A. Since the predictions of pK_a shifts in the carboxyl and histidine residues in these proteins proved particularly problematic in the previous work (16), they provide excellent test beds for the tautomeric state model as well as the GBSW solvent model.

The rest of the article is organized as follows. In the Theory section, following a brief review of the CPHMD methodology, the formalism of the two-dimensional -dynamics approach is outlined and discussed for two special cases where tautomerism manifests in either the protonated or the deprotonated state. In Methods, the derivations of model potential functions as well as the macroscopic pK_a values for the titration simulations involving histidine and carboxyl side chains are given. Force-field issues related to the construction of the protonation state models are discussed. The last subsection in Methods details the simulation protocol. In Results and Discussion, the convergence behavior and errors in the titration simulations of blocked aspartic acid and histidine are discussed. The pK_a values of OMTKY3 and RNase A obtained from the CPHMD simulations using the new model are analyzed and compared to those from the previous work of Lee et al. (16), from the static PB calculations and experimental data. The issues related to the force field and implicit solvation, as well as the dynamic coupling between protonation and the local conformation. The last section of the article draws conclusions and makes suggestions for future directions of research.

	THEORY

TOP ABSTRACT INTRODUCTION THEORY METHODS RESULTS AND DISCUSSION CONCLUDING REMARKS SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
Continuous constant pH molecular dynamics
To facilitate derivations and discussions of the new model, we briefly review the methodology of the continuous implicit solvent constant pH molecular dynamics (CPHMD) (16). In the CPHMD method, the protonation/deprotonation (titration) process at a protein titrating residue i is described by a set of continuous coordinates,

(1)

the end points of which define the deprotonated (_i = 1) and protonated states (_i = 0). Thus, the CPHMD simulation is a special case of -dynamics (17), where the -variables define the titration coordinates.

In classical simulations, the deprotonation free energy of a protein side chain is obtained through modeling the difference between the free energy of deprotonating a side chain in the protein environment and that in solvent. Since the latter energy, namely, the deprotonation energy of a blocked amino acid (model compound) in solvent, is experimentally accessible, the deprotonation energy in the protein environment can be obtained as

(2)

where G_exp(mod) is given as

(3)

and the pK_a value reflects the pH condition at which the populations of protonated and deprotonated states are equal.

Based on the above considerations, the total potential energy of the protein in the CPHMD simulation can be written as

(4)

where {_i;i = 1, n} is a set of titration coordinates for the n titrating residues as defined in Eq. 1. In Eq. 4, the two terms are the protonation-independent internal (bonded) energy and nonpolar solvation energy, and the rest are the protonation-dependent nonbonded and biasing potential energies.

The van der Waals energy for the interaction between the titrating proton i and another atom j is given by

(5)

where U^vdW(i, j) represents the protonation independent (full-strength) van der Waals interaction energy. The protonation dependence of the Coulomb and GB electrostatic energies is incorporated through the atomic partial charges on the titrating residue, which are linearly interpolated between the values in the deprotonated and protonated states,

(6)

where q_i, represents the partial charge of atom on residue i; and and are the charges in the deprotonated and protonated states, respectively.

The first biasing potential –U^mod in Eq. 4 originates from the term –G_class(mod) in Eq. 2 and is a potential of mean force (PMF) along the -coordinate. Due to the nature of the nonbonded terms (Eq. 5 and Eq. 6), it is a quadratic function of the titration coordinate,

(7)

and can be obtained via thermodynamic integration. The second biasing potential in Eq. 4 originates from the term G_exp(mod) in Eq. 2 and models the pH-dependence of the deprotonation free energy. Thus,

(8)

where pK_a(i) is the pK_a value of the titrating group i.

The last term in Eq. 4 is called a barrier potential,

(9)

where ß_i is the barrier at the center of the titration coordinate. The barrier potential does not alter the energy difference between the end-point states ( = 0 or = 1) but serves as an umbrella to prolong sampling time in the region of the end-point coordinates. This is necessary because the electrostatic and van der Waals interactions of the titrating residue associated with a mixed state (0 < < 1) are unphysical.

Two-dimensional -dynamics method
Tautomeric state model
As in most pK_a calculations, the existing CPHMD method assumes a particular protonating atom for residues that contain multiple titration sites, such as histidine and carboxyl groups. In what follows, a new model is introduced that removes this bias and treats the titration of the proton tautomers on an equal footing.

In the new model we assign an additional continuous coordinate x, which offers control of interconversion between the tautomeric states. This coordinate will be treated in the same manner as the titration coordinate . Thus, the titration of histidine involves interconversion among three states, which can be defined in terms of the titration and tautomeric coordinates as (0) for the doubly protonated state HSP; (1,1) for the NE2-protonated state HSE; and (1,0) for the ND1-protonated state HSD, where the first number in the parenthesis refers to the -value and the second the x value (Fig. 1 A). Notice that HSP can take on any x-value since its energy is degenerate with respect to x. Analogously, the three states in the titration of aspartic acid are defined as (1) for the deprotonated ASP; (0,1) for the OD1-protonated state ASPP1; and (0,0) for the OD2-protonated state ASPP2 (Fig. 1 B). Notice that ASP can take on any value of x.

View larger version (16K): [in this window] [in a new window]   FIGURE 1  Acid-base and proton tautomeric equilibria of histidine (A) and aspartic acid (B).

Nonbonded interactions
Following the addition of the tautomeric state coordinate x, the atomic charge on a titrating group can be written as a linear combination of four charge states,

(10)

where the values for the - and x-coordinates are given in parentheses. Since histidine has only one protonated state, q_i,(01) = q_i,(00). Analogously, for aspartic acid, q_i,(11) = q_i,(10). Notice that Eq. 6 can be recovered from Eq. 10 if the charge states of the tautomers are identical.

In the case of the deprotonated tautomerism (Fig. 1 A), the van der Waals energy for the interaction of titrating atom i and another atom j is given by

(11)

where the superscripts m and n indicate the nature of the titrating atoms i and j, respectively. For example, in the titration of histidine, a value of 1 is associated with the titrating proton HE2 and 0 is associated with HD1. For a nontautomeric titrating proton, a value of –1 is assigned. Thus, the function f^m(x_i) can be written as

(12)

where the function fⁿ(x_j) can be defined analogously. Notice that Eq. 5 can be recovered from Eq. 11 using the last definition of f^m(x_i), or fⁿ(x_i), in the absence of tautomerism.

In the case of the protonated tautomerism (Fig. 1 B), the van der Waals energy is given by

(13)

where the function f^m(x_i) (or fⁿ(x_i)) m (or n) takes on the value of 1 for the titrating proton HD1, 0 for HD2, and –1 for a nontautomeric titrating proton. Eq. 5 can again be recovered from Eq. 13 in the absence of tautomerism.

Biasing potentials
In the presence of proton tautomerism, the pH-dependent potential function becomes

(14)

where and refer to the microscopic pK_a values for the two titration sites. To suppress the population of mixed states in terms of x, a tautomer interconversion barrier potential in analogy to the titration barrier potential (Eq. 9) is applied to each coordinate x_i.

The functional form of the model potential in terms of the titration and tautomeric interconversion coordinates is determined by the functional forms of the Coulomb, GB electrostatic and van der Waals interaction energies. In the absence of tautomerism this leads to a quadratic function as mentioned earlier (Eq. 7). In the presence of proton tautomerism the coupling between and x results in a complex functional form for the model potential. However, since linear interpolation is used in the expressions for both the atomic partial charges and van der Waals interaction energy, the model potential function is expected to remain quadratic in either the - or x-variable.

Before attempting the derivation of the model potential function in the two-dimensional coordinate space (, x), it is instructive to consider the scenarios where either the titration or tautomeric interconversion coordinate is fixed at one of its end points. In the case of the deprotonated tautomerism, each process (or equilibrium) displayed in Fig. 1 A, is associated with a quadratic function. Since the energy of the protonated form, HSP, is independent of x, the potential function at HSP must be a constant. The following equations present the boundary conditions for the model potential function,

(15)

where A₁, B₁, A₀, B₀, A₁₀, and B₁₀ are the parameters in the quadratic functions that describe the one-dimensional processes and are related to each other through the relationship

(16)

In other words, the difference in the free energies of deprotonation for HSPHSE and HSPHSD is the same as the free energy of tautomer interconversion HSEHSD. We will show later that the model potential parameters can be derived analytically by making use of the boundary conditions given in Eq. 15. It should be noted that the model potential function contains an arbitrary constant since the difference deprotonation energy (see Eq. 2) is the interesting quantity in this work. For the case of the protonated tautomerism, there exists an analogous set of equations that describes the one-dimensional conditions.

Since the model potential function is quadratic in either the - or x-variable, a general expression in the form of a bivariate polynomial containing eight parameters can be found:

(17)

In the case where tautomerism exists for one protonation form (either U^mod(0, x_i) = 0 or U^mod(1, x_i) = 0) there are, at most, six nonvanishing parameters in Eq. 17. One way to identify the nonzero terms is to follow the procedure below,

• Step 1, at fixed values of x, obtain parameters for the quadratic function of : f(;x) = A(x)[ – B(x)]²;
  
• Step 2, fit the parameters A(x) and B(x) to the second-order polynomial: b₁x² + b₂x + b₃;
  
• Step 3, at fixed values of , obtain parameters for the quadratic function of x: f(x;) = A()[x – B()]²; and
  
• Step 4, fit the parameters A() and B() to the second-order polynomial: c₁² + c₂ + c₃,
  

	METHODS

TOP ABSTRACT INTRODUCTION THEORY METHODS RESULTS AND DISCUSSION CONCLUDING REMARKS SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
Macroscopic pK_a values from the tautomeric state model
In the CPHMD simulation, the fractional population of unprotonated states (unprotonated fraction) for a titrating residue is given as

(18)

where ^unprot and ^prot are the probabilities for the unprotonated and protonated states, which can be related to the number of times unprotonated (N^unprot) and protonated (N^prot) states are observed in the simulation, respectively.

Since bounded continuous coordinates are utilized for the propagation of protonation and tautomeric states, some cutoff values have to be used to define the pure protonation and tautomeric states. In this work we extend the convention employed previously by Lee et al. (16) to include the tautomeric states. Thus, the number of unprotonated and protonated states are defined as

(19)

Notice that the mixed tautomeric states (0.1 < x < 0.9) are discarded from the above definition, although the value of x is irrelevant for the doubly protonated state of histidine (Fig. 1 A) or the deprotonated carboxylate residue (Fig. 1 B). However, since the protonation states in the simulation are defined through some cutoff values, their energies do depend on x. Therefore, the mixed tautomeric states are not included in the statistics.

To obtain the pK_a value for a titrating residue, the unprotonated fractions at different pH were used to fit to the generalized form of the Henderson-Hasselbach equation,

(20)

where n is the Hill coefficient. The deviation of the value of n from unity reflects the degree of cooperativity (coupling) between groups that interact and ionize over the same pH range. In this work the Hill coefficient is included into the fitting procedure to account for the couplings between ionizable sites. It is noted, however, that the computed pK_a is only marginally affected by the value of n.

Given the equilibrium constants k₁ and k₂ for the microscopic deprotonation processes HSP HSE, and HSP HSD, respectively (Fig. 1 A), the equilibrium constant k for the macroscopic deprotonation of histidine is given as

(21)

Thus, by applying the definition of pK_a (pK_a = –log₁₀k, the macroscopic pK_a for the histidine titration can be related to the microscopic counterparts, pK₁ and pK₂ by

(22)

Provided the reference microscopic value of 6.6 (pK₁) and 7.0 (pK₂) for the ND1- and NE2-titration in the blocked histidine, respectively, Eq. 22 leads to a macroscopic pK_a of 6.45.

Given the equilibrium constants k₁ and k₂ for the microscopic deprotonation processes of aspartic acid (Fig. 1 B), the equilibrium constant k for the macroscopic deprotonation can be written as

(23)

which leads to the relation

(24)

Assuming a reference microscopic value 4.0 for both pK₁ and pK₂, Eq. 24 leads to a macroscopic pK_a of 4.3. Since the experimental macroscopic pK_a for the blocked aspartic acid is 4.0, a postcorrection is made (see later discussions).

Protonation state models
In the CPHMD method, the transition between protonation states of a titrating residue is modeled by linear attenuation of nonbonded energies while keeping the topology, bonded, and van der Waals parameters of the residue intact. Consequently, in a deprotonated state, the titrating hydrogen atom becomes a dummy atom with zero charge and has no van der Waals interaction with other atoms. This strategy is therefore similar to the single-topology scheme in free energy simulations.

The titration of histidine in the tautomeric state model is based on the topology and associated bonded as well as van der Waals parameters of the doubly protonated histidine (residue type HSP in CHARMM). The atomic charges of the three protonation states (Fig. 1 A) are identical to those in CHARMM type HSP, HSE, and HSD, respectively. In the intermediate or mixed states, the nonbonded interactions are computed by linear interpolation via the - and x-variables as described in the Theory section.

For the titrations of the carboxyl residue, a new CHARMM residue type is constructed based on carboxylic acid such that a second hydrogen atom is placed on the unprotonated oxygen. Let us consider aspartic acid. The titration of aspartic acid in the tautomeric state model is based on the topology, bonded and van der Waals parameters of the new hybrid residue, which will be referenced as ASDP. The force-field parameters of ASDP differ from those of aspartic acid (CHARMM type ASPP, protonated on OD2) in that the CG, OD1, and OD2 atoms are assigned with the corresponding bonded parameters of the deprotonated aspartate (CHARMM type ASP) whereas the titrating (dummy) hydrogens HD1 and HD2 are assigned with the bonded parameters of the aspartic acid (CHARMM type ASPP). From the hybrid residue ASDP, three protonation states are constructed: the deprotonated state ASP, with both dummy hydrogens uncharged and having no van der Waals interactions; two protonated states (ASP1 and ASP2), with either HD1 or HD2 charged and having van der Waals interactions. Notice that there is no interaction between the two dummy hydrogens due to the unphysical nature of the doubly protonated state, although both hydrogens are allowed to interact with other atoms in a mixed state.

One drawback of the dummy atom representation of the protonation states for the carboxyl residues is that the uncharged dummy hydrogen can lose the ability to gain charge once it is rotated to the anti position (16,14). A remedy used previously by Lee et al. (16) was to significantly lower the barrier to the rotation around the carboxylate bond such that the rotation back to the syn position is facilitated. A caveat is, however, that this can lead to nonnative hydrogen bonds as a result of stabilization of the unphysical out-of-plane conformations (16). The syn conformation is thermodynamically more favorable relative to the anti one by 2.8 kcal/mol, according to a solid-state NMR experiment (22). Quantum calculations and molecular dynamics simulations have shown a free energy difference of 1.6–1.9 kcal/mol (22,14). On the basis of these considerations, our solution is to kinetically disfavor the anti position by raising the C-O bond rotation barrier from 4.1 to 6.0 kcal/mol and placing the dummy carboxyl hydrogens at the syn positions at the beginning of the simulation.

Simulation protocol
All of the simulations described in this work were performed using the PHMD module with the all-atom CHARMM22 force field (24) for proteins including the dihedral cross term corrections (CMAP) (25) and the GBSW solvent model (20,19) in the CHARMM molecular dynamics program (26). The original PHMD module (16), in which the CPHMD method was implemented, was extended to allow tautomeric state titrations and to include pH-dependent solvation energy and force with the GBSW solvent model. The SHAKE algorithm was applied to hydrogen bonds and angles, such that a timestep of 2 fs could be used. A 22 Å distance truncation was applied to both the nonbonded and GB energy/force evaluations. All simulations utilized a Nosé-Hoover thermostat to ensure a canonical distribution of atomic velocities at a constant temperature of 298 K. For the GB calculations, a smoothing length of 0.6 Å at the dielectric boundary with 24 radial integration points up to 20 Å and 38 angular integration points were used. The nonpolar solvation energy was computed using the surface tension coefficient of 0.03 kcal mol^–1 Å^–2.

For the propagation of titration and tautomer interconversion coordinates, the velocities of and x were coupled to a three-mass (30, 50, and 70 amu) Nosé-Hoover chain (27,28) to ensure a canonical distribution at the target temperature of 298 K. Unless otherwise specified, the default initial velocity seed was used for the simulation. Model compounds used in this work were the blocked amino acids obtained by acetylating the N-terminus and amidating the C-terminus. To obtain the one-dimensional model potential functions, thermodynamic integrations based on 1-ns MD simulations (100 ps equilibration and 1 ns production) were performed at several fixed points along the titration or tautomer interconversion coordinates, as described previously (16). In single-site titration simulations, a titration barrier of 1.25 kcal mol^–1 was used (16). For double-site titrations, a titration barrier of 1.75 and a tautomer interconversion barrier of 2.75 kcal mol^–1 for the carboxyl groups and 2.5 kcal mol^–1 for histidines were applied, which in general keeps the fractional population of mixed titration and tautomeric states below 25% and 15%, respectively.

The microscopic reference pK_a values of 4.0, 4.4, and 3.8 were used for the Asp, Glu, and C-terminus carboxyl residues, respectively. This introduced a positive error of 0.3 units when both the carboxylate oxygens participated in the titration process (see Methods). Thus, we included a correction that accounts for this and other errors that occur in the titration simulations of model compounds in the reported pK_a values for the test proteins, OMTKY3 and RNase A (see the legends of Tables 3 and 5). For histidine, the microscopic reference values of 6.6 and 7.0 were used for the ND1- and NE2-titration, respectively. A correction that accounted for the error in the titration of blocked histidine was included in the reported pK_a values for the histidines in RNase A (see Table 5 legend). The reference pK_a value of 7.5 was used for the N-terminus amino group.

Development of atomic radii for use with the GB model is a continuous effort. Previous radii parameterizations for the PB methods were targeted toward the experimental solvation energies of small molecules as well as the calculated interaction energy of small molecules containing charged side chains with explicit waters. However, this strategy is not ideal, for at least two reasons. First, the solvation treatment of small molecules is not entirely transferable to a protein environment. Second, reproducing the explicit solvent behavior using the TIP3P (32), or any other explicit water model carries along inherent problems of these models (33,34). Thus, current effort in the group is directed toward reoptimization of the existing PB radii by performing control GB simulations of model peptides and mini-proteins aimed at finding the balance among stability, foldability, and the underlying physical principles (31). Since the strength of the electrostatic interactions in a microscopic approach such as the current CPHMD method is sensitive to changes in solvation or desolvation energies, the accuracy of the constant pH simulations will greatly benefit from the improvement of the GB model.

In the current method, a fixed set of solvation radii was used throughout the titration simulations. This is a drawback, because the radii of ND1 and NE2 on histidine in the GBSW model vary according to the charge state. Although for the neutral state (HSD and HSE), the radius is 1.8 Å, it is 0.5 Å larger in the charged state (HSP). A solution would be to linearly interpolate the radius between these two states, which would, however, result in a deviation from the quadratic behavior in the model potential function. Thus, in the present work, we assumed a radius of 2.0 Å for ND1 and NE2 in the pH range, expanding one unit below and one unit above the reference pK_a value, and the default radii for the neutral or charged states under other pH conditions. Since the pK_a values of carboxyl groups in OMTKY3 and RNase A are a few units below that of the histidine, the error due to the fixed histidine radii approach is minor.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION THEORY METHODS RESULTS AND DISCUSSION CONCLUDING REMARKS SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
Model compounds
Fig. 2 shows the two-dimensional potential of mean force (PMF) maps for the blocked histidine and aspartic acid along the titration and tautomer interconversion coordinates. Notice that the PMF of the doubly protonated histidine HSP (Fig. 2, A = 0) and the deprotonated aspartate (Fig. 2, B, = 1) is constant. From the PMF map for histidine (Fig. 2, A), it is seen that the ND1-protonated form HSD (lower-right corner where = 0, x = 1) is higher in energy than the NE2-protonated form HSE (upper-right corner where = 0, x = 0). The barrier for HSP HSD is higher and occurs later relative to that for HSP HSE. From the PMF map for the aspartic acid (Fig. 2, B), it is seen that the OD1-protonated form ASP1 and OD2-protonated form ASP2 have almost identical energy. The barriers for ASP1 ASP and ASP2 ASP are very small (<5 kcal/mol, not visible in Fig. 2) and occur very early in the deprotonation process.

View larger version (27K): [in this window] [in a new window]   FIGURE 2  Two-dimensional potential of mean force (PMF) map along the titration coordinate and tautomeric interconversion coordinate x for the blocked His (left) and Asp (right) residues. The color bars shown on the right indicate that the PMF increases from dark blue to dark red.

As mentioned in Methods, the difference between the solvation radii for ND1 and NE2 in the neutral and charged states of histidine is 0.5 Å. Since our approach uses a fixed set of radii for the -dependent GB energy and force, the protonation state whose assigned solvation radii are closer to the "true" values (as given in the PB radii set; see Refs. 28,29) may be artificially favored. For example, in the extreme case, if the true radius for the charged state HSP is used in the single-site titration simulations, the titrating nitrogen can become undersolvated resulting in a strong electrostatic attraction between the titrating proton and the partially negatively charged backbone carbonyl oxygen, which prevents the deprotonation to occur. Curiously, this has only occurred in titration simulations at the ND1 site. To better understand this bias, we performed a standard GBSW simulation of the blocked doubly protonated histidine. The top plot of Fig. 3 shows that up until 1 ns of MD simulation, the distance from the backbone oxygen to ND1 (solid) is 3.4 Å, which is >1 Å smaller than the distance to NE2 (shaded). At slightly after 1 ns, a conformational switch occurs, leading to the merge of the ND1-O and NE2-O distances. The middle plot of Fig. 3 shows a switch of the ₁ dihedral angle from 50° to 59° at 1 ns. The bottom plot of Fig. 3 shows that the ₂ dihedral angle remains approximately constant during the entire 2-ns simulation. Thus, the difference in the electrostatic environment between the ND1 and NE2 provides an explanation for the failure of titration at ND1. As the nitrogen radius for solvation is decreased to 2.0 Å, both the ND1-O and NE2-O distances increase but the former remains 1 Å smaller even after a conformational switch occurs at 1 ns (data not shown).

View larger version (43K): [in this window] [in a new window]   FIGURE 3  2-ns standard GBSW simulation of the blocked doubly protonated histidine. The top plot shows the time series of the distance from the backbone carbonyl oxygen to ND1 (solid) and that to NE2 (shaded). The middle and bottom plots show the time series of the 1 and 2 angles, respectively.

Blocked aspartic acid
Table 2 summarizes the running pK_a values for blocked aspartic acid in the double-site and single-site titration simulations. Given the target pK_a of 4.3 (see Methods), the double-site titrations lead to an underestimation of 0.5–0.6 pK units, with a standard deviation of 0.33 pK units in the first 2-ns simulation time. The magnitude of error and standard deviation decreases as the simulation continues. The last 2-ns gives an average error of –0.2 pK units with a standard deviation of 0.25.

If the above explanation were true, it would predict a positive error for single-site titrations using the bonded parameters of the neutral aspartic acid that has one single and one double C-O bond, as in the work of Lee et al. (16). This is indeed the case. The last two rows of Table 2 show a persistent overestimation by 0.2 pK units in the single-site titrations at OD1 or OD2 using the force-field parameters that are designed for aspartic acid. The same single-site model when used with the GBMV solvation method also gave a overestimation of 0.3 pK units in the predicted pK_a (16). Notice that the single-site titrations using the force-field parameters for the neutral species exhibit standard deviations twice those relative to using the parameters for the ionic species. Based on the above considerations and given the limitation of the current CPHMD methodology that allows only single topology, titration simulations of carboxyl group with the bonded parameters designed for the ionic species seems to be a better choice. It should be noted that the problem related to the force-field parameters does not exist in the titration simulations of histidine, since the latter has identical bonded parameters for both charge states in the CHARMM22 force field.

OMTKY3
Table 3 summarizes the computed pK_a values for turkey ovomucoid third domain (OMTKY3) from the double-site and single-site simulations and compares them with the results from experiment and other simulations. The double-site titration simulations with the crystal structure (column Cryst.) yield an average absolute error of 1.0 pK unit based on the first 2-ns simulations and 1.2 pK units based on the second 2-ns simulations. The double-site simulations with the NMR structure (column NMR) yield pK_a shifts (relative to the model compound reference values) of the same sign and the average absolute errors of similar magnitude. Thus we do not observe a strong initial structure dependence to the overall pK_a comparisons. The single-site OD1 and OD2 titrations with the crystal structure (column OD1 (OD2)) result in pK_a shifts of the same sign as from the double-site simulations, but much larger average absolute errors (1.9 and 1.5, respectively). As expected, the magnitude of these errors is similar to that from the single-site simulations with the GBMV solvent model (column GBMV). Compared to experimental data, the current simulations predict correct signs of pK_a shifts for all but one residue (Asp²⁷, see later discussions), whereas the GBMV simulations predict only one correct sign of the pK_a shift (Glu¹⁹). Finally, compared to the CPHMD simulations, the two sets of PB calculations employing a protein dielectric constant of 20 (column Static) based on the static crystal structure yield smaller average absolute errors (0.6 pK units).

To characterize the local electrostatic environment of titrating residues in OMTKY3, Table 4 summarizes the distances among titrating residues from the double-site simulations at three pH values and compares them with the distances measured in the crystal structure. These three pH values represent the conditions in which the specified residue is fully protonated (acidic), partially protonated, or fully deprotonated (basic). Below we will discuss the correlation between the pK_a value and the dominant electrostatic interactions involving the titrating residue as revealed by the simulations and the crystal structure data. In a related study, Li et al. (35) showed that pK_a values in OMTK3 could be reproduced using QM methods with minimal models including key interacting residues. It should be noted that since the crystal structure was obtained under basic (pH = 10) conditions, where all carboxyl groups were deprotonated, the measured distances between oppositely charged residues represent lower bounds to those measured at lower pH.

The pK_a (2.8) for Asp⁷ obtained from the first 2-ns simulations using the double-site model is between the values computed using the OD1 (pK_a = 3.9) and OD2 (pK_a = 1.9) single-site model, and agrees well with the experimental value of 2.7 (Table 3). This is consistent with the reasonable agreement of the simulated major electrostatic interactions (Asp⁷ with Glu¹⁰ and Lys³⁴) under basic conditions with those from experiment (crystal structure).

The electrostatic environment of Glu10 is more complex, and includes the titratable groups Asp⁷, Lys¹³, Arg²¹, and Lys³⁴ (Table 4). In the crystal structure, Lys¹³ is the closest residue with the Lys¹³-Glu¹⁰ distance of 4.8 Å, followed by Lys³⁴, Asp⁷, and Tyr¹¹, that are >6 Å away. Arg²¹ is 18 Å from Glu¹⁰ in the crystal structure. In contrast, MD simulations under neutral to basic conditions show that Arg²¹ moves to be <7 Å away from Glu¹⁰ (Table 4).

The top plot of Fig. 4 shows how the motion of Glu10 is correlated with that of Asp⁷, Lys³⁴, and Arg²¹ at pH 3, where Glu¹⁰ is partially deprotonated. At the beginning of the simulation, Lys³⁴ is in salt-bridge contact with Glu¹⁰ (dashed curve) and is 6.5 Å away from Asp⁷. Due to the attractive electrostatic force and the diminished mobility of Lys³⁴, Asp⁷ is driven closer to form a salt bridge with Lys³⁴ at 320 ps (shaded dashed curve), which also results in a sharp drop in the Asp⁷-Glu¹⁰ distance from 9 Å down to below 6 Å (solid curve). However, this local environment is not stable, due to the electrostatic repulsion between Asp⁷ and Glu¹⁰. At the same time, the electrostatic attraction between Arg²¹ and Glu¹⁰ facilitates the departure of Glu¹⁰ from Lys³⁴ and the formation of a loose salt-bridge contact between Glu¹⁰ and Arg²¹ at 450 ps (dotted curve), which then breaks off at 520 ps due to the motion of Glu¹⁰. The dynamics in the following period of time (until the end of the first half of the simulation), is dominated by the motion of Glu¹⁰, driven by the attraction from Lys³⁴ and Arg²¹. This continues at the beginning of the second half of the simulation and starting at 1.1 ns, Glu¹⁰ arrives at a position that allows it to form salt-bridge contacts with both Lys³⁴ and Arg²¹, although not as tight as the Lys³⁴-Asp⁷ salt bridge. The dynamics of the last part of the simulation (between 1.6 and 2 ns) is again characterized by the breaking and making of the electrostatic contacts of Glu¹⁰-Lys³⁴ and Glu¹⁰-Arg²¹.

View larger version (19K): [in this window] [in a new window]   FIGURE 4  Coupling between Asp7 and Glu10 in a 2-ns simulation of OMTKY3 at pH 3. The top plot shows the time series of the running average distances between Asp7 and Glu10 (solid), Asp7 and Lys34 (shaded dashed), Glu10 and Lys34 (dashed), and Glu10 and Arg21 (dotted) computed from 100-ps time windows. The definitions of the residue-residue distances are given under Table 4. The bottom plot shows the running unprotonated fraction for Asp7 (solid) and Glu10 (dashed) computed from 100-ps time windows.

In contrast to the simulation, which reveals two on and off salt bridges (Glu¹⁰···Lys³⁴ and Glu¹⁰···Arg²¹), the crystal structure shows only one salt bridge (Glu¹⁰···Lys¹³). This is the reason for the negative difference (–1.1 pK units) in the computed pK_a for Glu¹⁰ relative to experiment. In the last 2-ns simulations, the simultaneous formation of the two salt bridges gives rise to a larger difference (–1.9 pK units) relative to the first 2-ns simulations (–1.1 pK units). The salt-bridge effects seem to be more pronounced in single-site titrations. As a result, the differences are larger: –1.8 for the OD1 and –1.5 pK units for the OD2 titrations. Interestingly, the PB calculation also gives an underestimation of 0.8 pK units.

In the crystal structure under basic conditions, Glu¹⁹ interacts with three positively charged groups (Lys¹³, Arg²¹, and Lys³⁴) at a distance over 8 Å (Table 4). However, in the simulation at pH >2, Glu¹⁹ forms a persistent salt bridge with Arg²¹, thereby preventing the protonation of Glu¹⁹. This is the reason for the negative difference of Glu¹⁹'s pK_a (–0.9 pK units) using the double-site simulations (Table 3) relative to experiment. The magnitude of the underestimation using the OD1 and OD2 titrations is 2 pK units larger.

The prediction of the pK_a for Asp²⁷ is most problematic, since it is partially sequestered from solvent and yet it is the most acidic residue in OMTKY3 according to experiment (Table 4). In the simulation at pH 4, Asp²⁷ is a hydrogen-bond acceptor to its backbone amide nitrogen and to that of Lys²⁹ with an occupancy of 20%. In addition, Asp²⁷ interacts with the hydroxyl oxygen of Tyr³¹ as both a hydrogen-bond acceptor and donor with an occupancy of 30%. In the first 2-ns, Asp²⁷ is >40% deprotonated due to the salt-bridge interaction with Lys²⁹ (Table 4). Although the mobility of Asp²⁷ is low, the fluctuation of the side-chain position of Lys²⁹ is large as seen from the >1 Å root-mean square (RMS) deviation of the Asp²⁷-Lys²⁹ distance. In the last 2-ns, Asp²⁷ loses the salt-bridge contact with Lys²⁹ and dominantly occupies the protonated state (Table 4). The convergence of the protonation state sampling for Asp²⁷ is incomplete, which can be seen from a large percentage of mixed states (36% in the first and 42% in the last 2-ns of simulation time).

The double-site titration model overestimates the pK_a for Asp²⁷ by 2.4 pK units, the largest error among all the residues in OMTKY3 (Table 3). Interestingly, Asp²⁷ shows the largest computational error in the static PB calculations as well, although to a lesser degree (overestimated by 1.7 and 1.3 pK units). As compared to the single-site model, the double-site model does not offer any advantage for the titration of Asp²⁷. One obvious reason for the large deviation from experiment is the aforementioned issue of incomplete sampling of protonation states. Conformational relaxation of buried residues in response to ionization requires much longer simulation time as demonstrated in a recent work by Simonson et al. (36). Another possible reason for the difference may be the overestimation of the desolvation effect due to undersolvation for the protonated aspartic acid relative to the charged aspartate, since the current GB model employs one atomic solvation radius for the both the neutral and charged forms.

Glu⁴³ does not have strong electrostatic interactions with other charged groups and has therefore a small experimental pK_a upshift of 0.4 units (Table 3). Although the OD1 and OD2 titrations give a downshift of 0.4 and a upshift of 0.1 pK units, respectively, the double-site model predicts an upshift of 0.9 units, consistent with the experimental values. The deviations from experiment are within statistical errors, as shown in the titration simulations for model compounds. Interestingly, the static PB calculation is unable to capture the sign of the experimental pK_a shift of Glu⁴³.

The local electrostatic environment around the -carboxyl group CT-Cys includes two positively charged residues His⁵² and Lys⁵⁵ (Table 4). The pK_a downshift of CT-Cys is dominated by the interaction with His⁵², since the simulation data reveals the formation of a salt bridge between CT-Cys and His⁵² that does not exist in the crystal structure (Table 4). This explains why the double-site simulations result in a pK_a downshift that is 1.2-pK-units greater relative to experiment and the single-site titrations yield even larger errors.

RNase A
Bovine pancreatic ribonuclease A (RNase A) is an enzyme that catalyzes the transphosphorylation and hydrolysis reactions of RNAs (37). Table 5 summarizes the computed pK_a values from the double-site simulations with the GBSW solvation model (current model), in comparison with the split-model simulations (imidazole ring is split into two titration halves involving ND1 and NE2) with the GBMV solvation model (previous work) (16), and the PB static calculations (38), as well as experimental data. The current model gives an average absolute error of 0.6 pK units for the histidine groups, representing a significant improvement over the previous work using the split-model (2.5 pK units). In fact, the error from the current model is even smaller than that from the static PB calculations (0.8 pK units). It is noteworthy that results from the split-model simulations with the GBSW solvation model yield even larger differences (data not shown). The current model also gives smaller average absolute error (1.6 pK units) for the carboxyl groups relative to the previous work (2.0 pK units). For the N-terminus group, the current model yields a difference of 0.6 pK units, as compared to 1.0 pK unit from the previous work. Again, the current model predicts the correct sign for all pK_a shifts in RNase A, in contrast to the previous work that predicts the correct sign for only seven pK_a values (with a total of 17).

The protonation equilibria of active-site residues His¹² and His¹¹⁹ are responsible for the pH-dependent catalytic mechanism of RNase A. A traditional view is that His¹² and His¹¹⁹ act as general acid and base in RNA hydrolysis, respectively, whereas their roles are reversed in the transphosphorylation step (37). The crystal structure of the phosphate-free RNase A (PDB:7RSA, see Ref. 39; pH = 5.3), shows that His¹² is completely buried whereas His¹¹⁹ is partially buried in the protein interior. Thus, desolvation effects are expected to result in their pK_a downshifts. Distance measurements in the crystal structure reveal an electrostatic interaction between His¹² and His¹¹⁹ and a hydrogen bond between the NE2 of His¹¹⁹ and OD1 of Asp¹²¹ (Table 6).