constant-pH Simulations Research Articles

Constant pH Molecular Dynamics in Explicit Solvent with Enveloping Distribution Sampling and Hamiltonian Exchange

We present a new computational approach for constant pH simulations in explicit solvent based on the combination of the enveloping distribution sampling (EDS) and Hamiltonian replica exchange (HREX) methods. Unlike constant pH methods based on variable and continuous charge models, our method is based on discrete protonation states. EDS generates a hybrid Hamiltonian of different protonation states. A smoothness parameter s is used to control the heights of energy barriers of the hybrid-state energy landscape. A small s value facilitates state transitions by lowering energy barriers. Replica exchange between EDS potentials with different s values allows us to readily obtain a thermodynamically accurate ensemble of multiple protonation states with frequent state transitions. The analysis is performed with an ensemble obtained from an EDS Hamiltonian without smoothing, s = ∞, which strictly follows the minimum energy surface of the end states. The accuracy and efficiency of this method is tested on aspartic...

Effects of system net charge and electrostatic truncation on all-atom constant pH molecular dynamics.

Constant pH molecular dynamics offers a means to rigorously study the effects of solution pH on dynamical processes. Here, we address two critical questions arising from the most recent developments of the all-atom continuous constant pH molecular dynamics (CpHMD) method: (1) What is the effect of spatial electrostatic truncation on the sampling of protonation states? (2) Is the enforcement of electrical neutrality necessary for constant pH simulations? We first examined how the generalized reaction field and force-shifting schemes modify the electrostatic forces on the titration coordinates. Free energy simulations of model compounds were then carried out to delineate the errors in the deprotonation free energy and salt-bridge stability due to electrostatic truncation and system net charge. Finally, CpHMD titration of a mini-protein HP36 was used to understand the manifestation of the two types of errors in the calculated pK(a) values. The major finding is that enforcing charge neutrality under all pH conditions and at all time via cotitrating ions significantly improves the accuracy of protonation-state sampling. We suggest that such finding is also relevant for simulations with particle mesh Ewald, considering the known artifacts due to charge-compensating background plasma.

Constant pH Molecular Dynamics in Explicit Solvent with Enveloping Distribution Sampling and Hamiltonian Exchange.

We present a new computational approach for constant pH simulations in explicit solvent based on the combination of the enveloping distribution sampling (EDS) and Hamiltonian replica exchange (HREX) methods. Unlike constant pH methods based on variable and continuous charge models, our method is based on discrete protonation states. EDS generates a hybrid Hamiltonian of different protonation states. A smoothness parameter s is used to control the heights of energy barriers of the hybrid-state energy landscape. A small s value facilitates state transitions by lowering energy barriers. Replica exchange between EDS potentials with different s values allows us to readily obtain a thermodynamically accurate ensemble of multiple protonation states with frequent state transitions. The analysis is performed with an ensemble obtained from an EDS Hamiltonian without smoothing, s = ∞, which strictly follows the minimum energy surface of the end states. The accuracy and efficiency of this method is tested on aspartic acid, lysine, and glutamic acid, which have two protonation states, a histidine with three states, a four-residue peptide with four states, and snake cardiotoxin with eight states. The pKa values estimated with the EDS-HREX method agree well with the experimental pKa values. The mean absolute errors of small benchmark systems range from 0.03 to 0.17 pKa units, and those of three titratable groups of snake cardiotoxin range from 0.2 to 1.6 pKa units. This study demonstrates that EDS-HREX is a potent theoretical framework, which gives the correct description of multiple protonation states and good calculated pKa values.

Constant pH simulations with the coarse-grained MARTINI model — Application to oleic acid aggregates

Long chain fatty acids are biologically important molecules with complex and pH sensitive aggregation behavior. The carboxylic head group of oleic acid is ionizable, with the pKa shifting to larger values, even above a value of 7, in certain aggregate states. While experiments have determined the macroscopic phase behavior, we have yet to understand the molecular level details for this complex behavior. This level of detail is likely required to fully appreciate the role of fatty acids in biology and for nanoscale biotechnological and industrial applications. Here, we introduce the use of constant pH molecular dynamics (MD) simulations with the coarse-grained MARTINI model and apply the method to oleic acid aggregates and a model lipid bilayer. By running simulations at different constant pH values, we determined titration curves and the resulting pKa for oleic acid in different environments. The coarse-grained model predicts positive pKa shifts, with a shift from 4.8 in water to 6.5 in a small micelle, and 6.6 in a dioleoylphosphatidylcholine (DOPC) bilayer, similar to experimental estimates. The size of the micelles increased as the pH increased, and correlated with the fraction of deprotonated oleic acid. We show this combination of constant pH MD and the coarse-grained MARTINI model can be used to model pH-dependent surfactant phase behavior. This suggests a large number of potential new applications of large-scale MARTINI simulations in other biological systems with ionizable molecules.

Resolving Differences in Substrate Specificities between Human and Parasite Phosphoribosyltransferases via Analysis of Functional Groups of Substrates and Receptors

We herein review experimental and theoretical approaches widely applied to delineation of the differences in substrate specificities between human and parasite phosphoribosyltransferases (PRTases), the latter of which are key targets for treatment of diseases caused by parasites. Standard Molecular Dynamics (MD) simulations have been applied to determine why the human PRTase prefers guanine over xanthine, whereas the Tritrichomonas foetus enzyme exhibits only a slight preference. We analyze this problem with the aid of standard MD simulations, as well as constant-pH MD simulations. Comparison of results of the two approaches reveals substantial differences, e.g. several Asp and Glu residues in the parasite enzyme, and one Glu residue in the human enzyme, are predicted to be permanently or frequently protonated during constant-pH simulations, whereas standard MD simulations assume that these residues are always ionized. Most interesting is the observation of a large conformational change, leading to tighter binding of the ligand, observed in constant-pH MD simulations of the parasite PRTase complexed with XMP, and lack of such a change in the human enzyme complexed with XMP.

Capturing molten globule state of α-lactalbumin through constant pH molecular dynamics simulations

The recently developed methods of constant pH molecular dynamics directly captures the correlation between protonation and conformation to probe protein structure, function, and dynamics. In this work, we investigate the effect of pH on the conformational properties of the protein human α-lactalbumin. Constant pH simulations at both acidic and alkaline medium indicate the formation of the molten globule state, which is in accordance with the previous experimental observations (especially, in acidic medium). The size of the protein measured by its radius of gyration (RG) exhibits a marked increase in both acidic and alkaline medium, which matches with the corresponding experimentally observed value of RG found in the molten globule. The probability of native contacts is also considerably reduced at acidic and basic pH as compared to that of native structure crystallized at neutral pH. The mean fractal dimension D2 of the protein records a sharp increase in basic medium as compared to those in neutral and acidic solutions implying a significant pH induced conformational change. The mean square fluctuations of all residues of the entire protein are found to increase by several folds in both acidic and basic medium, which may be correlated with the normalized solvent accessibility of the residues indicating role of solvent accessible surface area on protein internal dynamics. The helices comprising the α-domain of the protein are moderately preserved in the acidic and alkaline pH. However, the β-sheet structures present in the β-domain are completely disrupted in both acidic as well as basic pH.

PH-dependent dynamics of complex RNA macromolecules.

The role of pH-dependent protonation equilibrium in modulating RNA dynamics and function is one of the key unanswered questions in RNA biology. Molecular dynamics (MD) simulations can provide insight into the mechanistic roles of protonated nucleotides, but it is only capable of modeling fixed protonation states and requires prior knowledge of the key residue's protonation state. Recently, we developed a framework for constant pH molecular dynamics simulations (CPHMDMSλD) of nucleic acids, where the nucleotides' protonation states are modeled as dynamic variables that are coupled to the structural dynamics of the RNA. In the present study, we demonstrate the application of CPHMDMSλD to the lead-dependent ribozyme; establishing the validity of this approach for modeling complex RNA structures. We show that CPHMDMSλD accurately predicts the direction of the pKa shifts and reproduces experimentally-measured microscopic pKa values with an average unsigned error of 1.3 pKa units. The effects of coupled titration states in RNA structures are modeled, and the importance of conformation sampling is highlighted. The general accuracy of CPHMDMSλD simulations in reproducing pH-dependent observables reported in this work demonstrates that constant pH simulations provides a powerful tool to investigate pH-dependent processes in nucleic acids.

Fatty Acid Aggregates Simulated using Constant pH Molecular Dynamics with a Coarse-Grained Model

Fatty acids are crucial biomolecules, important for lipid metabolism, signaling, models for protocell membranes, soaps, industrial applications, and drug delivery. Oleic acid has complex phase behavior with respect to the protonation state of the carboxylic head group, which depends on the pH of the solution. Oils form at low pHs, vesicles at intermediate pHs, and micelles at high pHs. We use constant pH molecular dynamics with the MARTINI coarse-grained model to investigate oleic acid aggregates at different pH conditions. We determine titration curves for the oleic acid monomers in different aggregates, and observe a shift in the microscopic pKa. In agreement with experimental results, the pKa of a monomer in bulk water is ca. 4 and shifts to ca. 5.5 in a small micelle and ca. 8-9 in a fatty acid bilayer. There is strong anti-cooperative protonation behavior between monomers within an aggregate. This work presents a proof of concept for using constant pH simulations with the MARTINI model, and provides a physiochemical basis for the phase behavior of fatty acids in different aggregate environments.

Efficient Dynamic Protonation and Constant pH Simulations with Explicit Solvent: Calculation of Apparent pKa Values in Proteins

The pKa's of the ionizable amino acids are crucial for the function of many proteins as they are key factors that determine their electrostatic potential and its spatial distribution, often controlling and optimizing enzymatic catalysis. Further, during conformational motions pKa's and protonation states particularly of histidines may change. In established force field simulation, however, this effect is typically not included, and protonation states must therefore be either guessed or derived from experiment. There have been a number of approaches to include protonation effects within simulations, mainly based on continuum electrostatics or implicit solvent molecular dynamics [1--3]. However, these methods lack the effect of the hydrogen bonding and the entropy contribution that comes from the solvent. Here we present the implementation and application of a dynamic protonation atomistic simulation method with explicit solvent, which also allows for explicit solvent constant pH MD simulations, previously developed also in our group [4]. This method is used here to calculate the pKa's of the ionzable groups in proteins. In order to validate it, we selected a number of prototypic proteins and calculated titration curves and pKa values from constant pH simulations at a range of different pH values. The results compare favorably with measured values, and explain atomistically the strong deviations of some of the calculated pKa values from the solution ones. 1. Lee, M. S., Salsbury, F. R., Jr., and Brooks, C. L., III (2004), Proteins 56, 738–752. 2. Khandogin J, Brooks C. L., 3rd., Biophys J. 2005 Jul;89(1):141-57. Epub 2005 Apr, 29. 3. Mongan, J., Case, D. A., McCammon, J. A. J., Comput. Chem. 2004, 25, 2038–2048 4. Donnini S, Tegeler F, Groenhof G, Grubmüller H., J. Chem Theory and Comp 7: 1962–1978 (2011).

PH Replica-Exchange Method Based on Discrete Protonation States

We propose a new algorithm for obtaining proton titration curves of ionizable residues. The algorithm is a pH replica-exchange method (PHREM), which is based on the constant pH algorithm of Mongan et al. (J Comput Chem 2004;25:2038-2048). In the original replica-exchange method, simulations of different replicas are performed at different temperatures, and the temperatures are exchanged between the replicas. In our PHREM, simulations of different replicas are performed at different pH values, and the pHs are exchanged between the replicas. The PHREM was applied to a blocked amino acid and to two protein systems (snake cardiotoxin and turkey ovomucoid third domain), in conjunction with a generalized Born implicit solvent. The performance and accuracy of this algorithm and the original constant pH method (PHMD) were compared. For a single set of simulations at different pHs, the use of PHREM yields more accurate Hill coefficients of titratable residues. By performing multiple sets of constant pH simulations started with different initial states, the accuracy of predicted pKa values and Hill coefficients obtained with PHREM and PHMD methods becomes comparable. However, the PHREM algorithm exhibits better samplings of the protonation states of titratable residues and less scatter of the titration points and thus better precision of measured pKa values and Hill coefficients. In addition, PHREM exhibits faster convergence of individual simulations than the original constant pH algorithm.

PH replica‐exchange method based on discrete protonation states

We propose a new algorithm for obtaining proton titration curves of ionizable residues. The algorithm is a pH replica-exchange method (PHREM), which is based on the constant pH algorithm of Mongan et al. (J Comput Chem 2004;25:2038-2048). In the original replica-exchange method, simulations of different replicas are performed at different temperatures, and the temperatures are exchanged between the replicas. In our PHREM, simulations of different replicas are performed at different pH values, and the pHs are exchanged between the replicas. The PHREM was applied to a blocked amino acid and to two protein systems (snake cardiotoxin and turkey ovomucoid third domain), in conjunction with a generalized Born implicit solvent. The performance and accuracy of this algorithm and the original constant pH method (PHMD) were compared. For a single set of simulations at different pHs, the use of PHREM yields more accurate Hill coefficients of titratable residues. By performing multiple sets of constant pH simulations started with different initial states, the accuracy of predicted pK(a) values and Hill coefficients obtained with PHREM and PHMD methods becomes comparable. However, the PHREM algorithm exhibits better samplings of the protonation states of titratable residues and less scatter of the titration points and thus better precision of measured pK(a) values and Hill coefficients. In addition, PHREM exhibits faster convergence of individual simulations than the original constant pH algorithm.

Constant pH Molecular Dynamics in Explicit Solvent with λ-Dynamics.

pH is an important parameter in condensed-phase systems, because it determines the protonation state of titratable groups and thus influences the structure, dynamics, and function of molecules in solution. In most force field simulation protocols, however, the protonation state of a system (rather than its pH) is kept fixed and cannot adapt to changes of the local environment. Here, we present a method, implemented within the MD package GROMACS, for constant pH molecular dynamics simulations in explicit solvent that is based on the λ-dynamics approach. In the latter, the dynamics of the titration coordinate λ, which interpolates between the protonated and deprotonated states, is driven by generalized forces between the protonated and deprotonated states. The hydration free energy, as a function of pH, is included to facilitate constant pH simulations. The protonation states of titratable groups are allowed to change dynamically during a simulation, thus reproducing average protonation probabilities at a certain pH. The accuracy of the method is tested against titration curves of single amino acids and a dipeptide in explicit solvent.

In Silico Titration of Biomolecules: Explicit Solvent Constant pH Molecular Dynamics

The pH is an important parameter in macromolecular systems as it determines the protonation state of ionizable groups and consequently influences the structure, dynamics and function of molecules in solution. In most force field simulation protocols, however, the protonation state of a system (rather than its pH) is kept fix and cannot adapt to changes of the local environment. Here, we present a method to perform molecular dynamics simulations in explicit solvent at constant pH. During the simulation the protonation states of titratable groups are allowed to change dynamically, and the titration curves agree with experiment. Our method is based on the lambda-dynamics approach, in which the dynamics of the titration coordinate lambda is driven by generalized forces between the protonated and deprotonated states. Constant pH simulations can be achieved by accounting for the pH dependence of the hydration free energy. As a benchmark, titration curves of amino acid analogues and a di-peptide, as well as of turkey ovomucoid inhibitor protein were calculated.

Effects of pH on Conformational Equilibria of Intrinsically Disordered Proteins

Intrinsically disordered proteins (IDPs) adopt an ensemble of conformations under native physiological conditions. Despite their lack of folded structure, they perform important physiological functions and are predicted to constitute around 30% of the eukaryotic proteome. The success of disorder prediction using a protein's primary structure suggests that the propensity for disorder is encoded in the amino acid sequence. Recently, we found that net charge per residue segregates IDP sequences along a globule-to-coil transition, enabling construction of a sequence-space phase diagram that subdivides IDPs based on the polymeric character of their conformational ensembles. Here, we explore the effects of two perturbations that leave sequence composition unchanged: sequence permutation and pH titration. Using atomistic Monte Carlo simulations in ABSINTH implicit solvent, we find that permutants of an arginine-rich IDP sequence exhibit nearly identical polymeric attributes, yet differ in the details of their conformational ensembles. In contrast, fluorescence correlation spectroscopy shows a decrease in translational diffusion times with increasing pH, suggesting that electrostatic and conformational modulation of protonation equilibria is significant even for solvent-exposed titratable groups on flexible protein backbones. The effect is sufficient to induce a collapse transition in a polyarginine polypeptide at a pH well below the pKa of the arginine sidechain. We conclude that the polymeric character of an IDP is largely determined by its sequence composition, but emphasize the importance of accurate constant pH simulation technology in investigating details of its conformational equilibria.Supported by grants 5T32GM008802 from the NIH National Institute of General Medical Sciences and MCB 0718924 from the National Science Foundation.

Alchemical free energy simulations for biological complexes: powerful but temperamental …

Free energy simulations compare multiple ligand:receptor complexes by "alchemically" transforming one into another, yielding binding free energy differences. Since their introduction in the 1980s, many technical and theoretical obstacles were surmounted, and the method ("MDFE," since molecular dynamics are often used) has matured into a powerful tool. We describe its current status, its effectiveness, and the challenges it faces. MDFE has provided chemical accuracy for many systems but remains expensive, with significant human overhead costs. The bottlenecks have shifted, partly due to increased computer power. To study diverse sets of ligands, force field availability and accuracy can be a major difficulty. Another difficulty is the frequent need to consider multiple states, related to sidechain protonation or buried waters, for example. Sophisticated, automated methods to sample these states are maturing, such as constant pH simulations. Meanwhile, combinations of MDFE and simpler approaches, like continuum dielectric models, can be very effective. As illustrations, we show how, with careful force field parameterization, MDFE accurately predicts binding specificities between complex tetracycline ligands and their targets. We describe substrate binding to the aspartyl-tRNA synthetase enzyme, where many distinct electrostatic states play a role, and a histidine and a Mg(2+) ion act as coupled switches that help enforce a strict preference for the aspartate substrate, relative to several analogs. Overall, MDFE has achieved a predictive status, where novel ligands can be studied and molecular recognition elucidated in depth. It should play an increasing role in the analysis of complex cellular processes and biomolecular engineering.

Constant pH Simulations In Explicit Solvent Using The Lambda-Dynamics Approach

pH is an important parameter in condensed phase systems as it determines the protonation state of ionizable groups and consequently influences the structure, dynamics and function of molecules in solution. In the past ten years, few approaches have been applied to model the pH of a solution in the framework of Molecular Dynamics (MD) and Monte Carlo (MC) simulation methods. These include stochastic and mean field approximation methods to model the (de)protonation events and methods based on the lambda-dynamics approach, where the dynamics of the titration coordinate lambda is driven by the free energy gradient between the protonated and deprotonated states. In particular, the latter approach was so far limited to implicit solvent. We present here a method for constant pH simulations in explicit solvent that is based on the lambda-dynamics approach. The method has been implemented in the MD package Gromacs. The titration curve of single amino acids and small peptides and the shift in the pKa of an active site glutamic acid in the enzyme triosephosphate isomerase were correctly predicted. This preliminary tests suggest that the approach can be applied to simulate (bio)molecules with multiple titrating sites at constant pH.

Molecular simulation with variable protonation states at constant pH

A new method is presented for performing molecular simulations at constant pH. The method is a Monte Carlo procedure where trial moves consist of relatively short molecular dynamics trajectories, using a time-dependent potential energy that interpolates between the old and new protonation states. Conformations and protonation states are sampled from the correct statistical ensemble independent of the trial-move trajectory length, which may be adjusted to optimize efficiency. Because moves are not instantaneous, the method does not require the use of a continuum solvation model. It should also be useful in describing buried titratable groups that are not directly exposed to solvent, but have strong interactions with neighboring hydrogen bond partners. The feasibility of the method is demonstrated for a simple test case: constant-pH simulations of acetic acid in aqueous solution with an explicit representation of water molecules.

Effects of Solute−Solvent Proton Exchange on Polypeptide Chain Dynamics: A Constant-pH Molecular Dynamics Study

A method for performing implicit-solvent molecular dynamics simulations at constant pH was applied to a pentapeptide acetyl-Ala-Asp-Ala-Lys-Ala-amide at pH 4. As a reference, molecular dynamics simulations were done for the same peptide with two variants of its fixed protonation patterns expected to dominate at pH 4, i.e., with a protonated and a deprotonated side chain of the Asp residue and the protonated Lys residue in both cases. The dynamic trajectories of the peptide were used to discuss the problem of the significance of the solute-solvent proton exchange phenomena for the dynamics and structural distributions of the polypeptide chain. The Asp-Lys distance was used as a probe of the overall molecular structure of the investigated pentapeptide. To characterize the dynamics, distributions of the "waiting" times for a transition from a "short" distance conformation to a "long" distance conformation were constructed, based on the generated molecular dynamics trajectories. We show that the relaxation time for the transitions, derived from the constant-pH simulations, is very close to the relaxation time characterizing a permanently protonated molecule, although the average protonation probability of the short-distance conformation is close to zero. However, the distribution of the Asp-Lys distances obtained from constant-pH simulations cannot be reproduced as a linear combination of the distributions resulting from the simulations with fixed protonation states.

PH-Dependent Stability of a Decalysine α-Helix Studied by Explicit-Solvent Molecular Dynamics Simulations at Constant pH

The acidostat method previously developed for performing explicit-solvent molecular dynamics simulations at constant pH (J. Chem. Phys. 2001, 114, 9706) is applied to polyfunctional compounds, namely 1,4-diaminobutane and a decalysine peptide. The titration behavior of 1,4-diaminobutane is investigated by performing a series of simulations at different pH, using the acidostat method. The method accounts at least to some extent for site−site coupling and reproduces the experimental pKa values of the compound within half a pK unit, although the simulations reveal insufficient sampling of the protonation-state variables. In a second step, the ability of the acidostat method to account for correlations between the solution pH and the structure and dynamics of a biomolecule is tested by studying the pH-dependent stability of an α-helical decalysine peptide. To this end, four 32-ns constant-pH simulations at different pH values are performed. The results are compared to those of standard molecular dynamics simu...