Langevin dynamics of proteins at constant pH

Constant pH Simulation with FMM Electrostatics in GROMACS. (A) Design and Applications

The structural dynamics of biological macromolecules, such as proteins, DNA/RNA, or complexes thereof, are strongly influenced by protonation changes of their typically many titratable groups, which explains their sensitivity to pH changes. Conversely, conformational and environmental changes of the biomolecule affect the protonation state of these groups. With few exceptions, conventional force field-based molecular dynamics (MD) simulations neither account for these effects nor do they allow for coupling to a pH buffer. Here, we present design decisions and applications of a rigorous Hamiltonian interpolation λ-dynamics constant pH method in GROMACS, which rests on GPU-accelerated Fast Multipole Method (FMM) electrostatics. Our implementation supports both CHARMM36m and Amber99sb*-ILDN force fields and is largely automated to enable seamless switching from regular MD to constant pH MD, involving minimal changes to the input files. Here, the first of two companion papers describes the underlying constant pH protocol and sample applications to several prototypical benchmark systems such as cardiotoxin V, lysozyme, and staphylococcal nuclease. Enhanced convergence is achieved through a new dynamic barrier height optimization method, and high pKa accuracy is demonstrated. We use Functional Mode Analysis (FMA) and Mutual Information (MI) to explore the complex intra- and intermolecular couplings between the protonation states of titratable groups as well as those between protonation states and conformational dynamics. We identify striking conformation-dependent pKa variations and unexpected inter-residue couplings. Conformation-protonation coupling is identified as a primary cause of the slow protonation convergence notorious to constant pH simulations involving multiple titratable groups, suggesting enhanced sampling methods to accelerate convergence.

Constant pH Simulation with FMM Electrostatics in GROMACS. (B) GPU Accelerated Hamiltonian Interpolation

The structural dynamics of biological macromolecules, such as proteins, DNA/RNA, or their complexes, are strongly influenced by protonation changes of their typically many titratable groups, which explains their pH sensitivity. Conversely, conformational and environmental changes in the biomolecule affect the protonation state of these groups. With a few exceptions, conventional force field-based molecular dynamics (MD) simulations do not account for these effects, nor do they allow for coupling to a pH buffer. The λ-dynamics method implements this coupling and thus allows for MD simulations at constant pH. It uses separate Hamiltonians for the protonated and deprotonated states of each titratable group, with a dynamic λ variable that continuously interpolates between them. However, rigorous implementations of Hamiltonian Interpolation (HI) λ-dynamics are prohibitively slow for typical numbers of sites when used with particle mesh Ewald (PME). To circumvent this problem, it has recently been proposed to interpolate the charges (QI) instead of the Hamiltonians. Here, in the second of two companion papers, we propose a rigorous yet efficient Multipole-Accelerated Hamiltonian Interpolation (MAHI) method to perform λ-dynamics in GROMACS. Starting from a charge-scaled Hamiltonian, precomputed with the Fast Multipole Method (FMM), the correct HI forces are calculated with negligible computational overhead. However, other electrostatic solvers, such as PME, can also be used for the precomputation. We compare Hamiltonian interpolation with charge interpolation and show that HI leads to more frequent transitions between protonation states, resulting in better sampling and accuracy. Our accuracy and performance benchmarks show that introducing, e.g., 512 titratable sites to a one million atom MD system increases runtime by less than 20% compared to a regular FMM-based simulation. We have integrated the scheme into our GPU-accelerated FMM code for the simulation software GROMACS, allowing easy and effortless transitions from standard force field simulations to constant pH simulations.

Titration of Adenine in a GA Mismatch with Grand Canonical Simulations

Computational prediction of the pKa of ionizable groups remains a central challenge in biomolecular modeling. Although all-atom fixed-charge force fields could be accurate to describe the interaction network within the biomolecules, proper sampling techniques are required to obtain the thermodynamic information in the (de)protonation event. Sufficient sampling requires an ensemble of structures from simulations with proper treatments of the acid-base equilibria, and the grand canonical simulation technique could be used to model the growth/annihilation of hydrogen atoms by merging Hamiltonians of different protonation states into one simulation ensemble. The electrostatic feature of nucleotide systems is especially difficult to model, and the situation becomes more challenging when the ionizable site is highly perturbed. Although there are many successful predictions obtained from the grand canonical constant pH simulations, few reports focus on highly perturbed nucleotide systems with unconventional base-pair features. In this work, with the discrete constant pH method, we investigate the titration thermodynamics of an adenine in the catalytic triad in a 35-nucleotide single-stranded RNA hairpin, featuring an unconventional GA mismatch and a substantially shifted pKa value. Validation tests are performed with two system setups, both of which provide pKa predictions in good agreement with the experimental value. A single-configuration-based technique is used to calculate the pKa for comparison. The current success indicates the predictive power of the current nucleotide modeling framework.

Wang-Landau Reaction Ensemble Method: Simulation of Weak Polyelectrolytes and General Acid-Base Reactions

We present a novel method for the study of weak polyelectrolytes and general acid-base reactions in molecular dynamics and Monte Carlo simulations. The approach combines the advantages of the reaction ensemble and the Wang-Landau sampling method. Deprotonation and protonation reactions are simulated explicitly with the help of the reaction ensemble method, while the accurate sampling of the corresponding phase space is achieved by the Wang-Landau approach. The combination of both techniques provides a sufficient statistical accuracy such that meaningful estimates for the density of states and the partition sum can be obtained. With regard to these estimates, several thermodynamic observables like the heat capacity or reaction free energies can be calculated. We demonstrate that the computation times for the calculation of titration curves with a high statistical accuracy can be significantly decreased when compared to the original reaction ensemble method. The applicability of our approach is validated by the study of weak polyelectrolytes and their thermodynamic properties.

Origin of pKa Shifts of Internal Lysine Residues in SNase Studied Via Equal-Molar VMMS Simulations in Explicit Water

Protein internal ionizable groups can exhibit large shifts in pK_a values. Although the environment and interaction changes have been extensively studied both experimentally and computationally, direct calculation of pK_a values of these internal ionizable groups in explicit water is challenging due to energy barriers in solvent interaction and in conformational transition. The virtual mixture of multiple states (VMMS) method is a new approach designed to study chemical state equilibrium. This method constructs a virtual mixture of multiple chemical states in order to sample the conformational space of all states simultaneously and to avoid crossing energy barriers related to state transition. By applying VMMS to 25 variants of staphylococcal nuclease with lysine residues at internal positions, we obtained the pK_a values of these lysine residues and investigated the physics underlining the pK_a shifts. Our calculation results agree reasonably well with experimental measurements, validating the VMMS method for pK_a calculation and providing molecular details of the protonation equilibrium for protein internal ionizable groups. Based on our analyses of protein conformation relaxation, lysine side chain flexibility, water penetration, and the microenvironment, we conclude that the hydrophobicity of the microenvironment around the lysine side chain (which affects water penetration differently for different protonation states) plays an important role in the pK_a shifts.

Coupling between protonation and conformation in cytochrome c oxidase: Insights from constant-pH MD simulations

Cytochrome c oxidases (CcOs) are the terminal enzymes of the respiratory chain in mitochondria and most bacteria. These enzymes reduce dioxygen (O(2)) to water and, simultaneously, generate a transmembrane electrochemical proton gradient. Despite their importance in the aerobic metabolism and the large amount of structural and biochemical data available for the A1-type CcO family, there is still no consensually accepted description of the molecular mechanisms operating in this protein. A substantial number of questions about the CcO's working mechanism remain to be answered, including how the protonation behavior of some key residues is modulated during a reduction cycle and how is the conformation of the protein affected by protonation. The main objective of this work was to study the protonation-conformation coupling in CcOs and identify the molecular factors that control the protonation state of some key residues. In order to directly capture the interplay between protonation and conformational effects, we have performed constant-pH MD simulations of an A1-type CcO inserted into a lipid bilayer in two redox states (oxidized and reduced) at physiological pH. From the simulations, we were able to identify several groups with unusual titration behavior that are highly dependent on the protein redox state, including the A-propionate from heme a and the D-propionate from heme a3, two key groups possibly involved in proton pumping. The protonation state of these two groups is heavily influenced by subtle conformational changes in the protein (notably of R481(I) and R482(I)) and by small changes in the hydrogen bond network.

A Virtual Mixture Approach to the Study of Multistate Equilibrium: Application to Constant pH Simulation in Explicit Water

Chemical and thermodynamic equilibrium of multiple states is a fundamental phenomenon in biology systems and has been the focus of many experimental and computational studies. This work presents a simulation method to directly study the equilibrium of multiple states. This method constructs a virtual mixture of multiple states (VMMS) to sample the conformational space of all chemical states simultaneously. The VMMS system consists of multiple subsystems, one for each state. The subsystem contains a solute and a solvent environment. The solute molecules in all subsystems share the same conformation but have their own solvent environments. Transition between states is implicated by the change of their molar fractions. Simulation of a VMMS system allows efficient calculation of relative free energies of all states, which in turn determine their equilibrium molar fractions. For systems with a large number of state transition sites, an implicit site approximation is introduced to minimize the cost of simulation. A direct application of the VMMS method is for constant pH simulation to study protonation equilibrium. Applying the VMMS method to a heptapeptide of 3 ionizable residues, we calculated the pKas of those residues both with all explicit states and with implicit sites and obtained consistent results. For mouse epidermal growth factor of 9 ionizable groups, our VMMS simulations with implicit sites produced pKas of all 9 ionizable groups and the results agree qualitatively with NMR measurement. This example demonstrates the VMMS method can be applied to systems of a large number of ionizable groups and the computational cost scales linearly with the number of ionizable groups. For one of the most challenging systems in constant pH calculation, SNase Δ+PHS/V66K, our VMMS simulation shows that it is the state-dependent water penetration that causes the large deviation in lysine66’s pKa.

Computer simulation plays an important role to understand molecular systems and has been applied to problems of increasing complexity. Multistate equilibrium is a fundamental concept behind the structure and function of biological systems. Due to the limit in computing resources and lack of good alternative methods, computer simulation has been conducted for systems in a single state, sampling from one state to another to infer equilibrium properties. This sequential approach has been successful in many cases such as protonation equilibrium with implicit solvation model. However, state transition is difficult when explicit solvent is used for more accurate solvation description. Many efforts have been dedicated to overcome this difficulty. Analogous to real multistate systems, we proposed a virtual mixture of multiple states (VMMS) to directly simulate the equilibrium. State transitions are replaced by changes in state molar fractions. Mimicking a test tube environment, all states are simulated in parallel to equilibrate with each other. Application to constant pH simulation in explicit water demonstrates the capability of this method. It is expected that the VMMS method will find more applications in biological problems related to the equilibrium of competing states.

Constant pH Replica Exchange Molecular Dynamics in Explicit Solvent Using Discrete Protonation States: Implementation, Testing, and Validation

By utilizing Graphics Processing Units, we show that constant pH molecular dynamics simulations (CpHMD) run in Generalized Born (GB) implicit solvent for long time scales can yield poor pK_a predictions as a result of sampling unrealistic conformations. To address this shortcoming, we present a method for performing constant pH molecular dynamics simulations (CpHMD) in explicit solvent using a discrete protonation state model. The method involves standard molecular dynamics (MD) being propagated in explicit solvent followed by protonation state changes being attempted in GB implicit solvent at fixed intervals. Replica exchange along the pH-dimension (pH-REMD) helps to obtain acceptable titration behavior with the proposed method. We analyzed the effects of various parameters and settings on the titration behavior of CpHMD and pH-REMD in explicit solvent, including the size of the simulation unit cell and the length of the relaxation dynamics following protonation state changes. We tested the method with the amino acid model compounds, a small pentapeptide with two titratable sites, and hen egg white lysozyme (HEWL). The proposed method yields superior predicted pK_a values for HEWL over hundreds of nanoseconds of simulation relative to corresponding predicted values from simulations run in implicit solvent.

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.

Conformational changes of enzymes and DNA in molecular dynamics: influenced by pH, temperature, and ligand

Protein conformation, which has been a research hotspot for human diseases, is an important factor of protein properties. Recently, a series of approaches have been utilized to investigate the conformational changes under different conditions. Some of them have gained promising achievements, but it is still deficient in the detail researches at the atomic level. In this chapter, a series of computational examples of protein conformational changes under different pH environment, temperature, and ligand binding are described. We further show some useful methods, such as constant pH molecular dynamics simulations, molecular docking, and molecular mechanics Poisson-Boltzmann surface area/generalized Born surface area calculations. In comparison with the experimental results, the methods mentioned above are reasonable to detect and predict the interaction between residue and residue, residue and DNA, and residue and ligand. Additionally, some crucial interactions that cause protein conformational changes are discovered and discussed in this chapter. In summary, our work can give penetrating information to understand the pH-, temperature-, and ligand-induced conformational change mechanisms.

Charge parametrization of the DvH-c3 heme group: validation using constant-(pH,E) molecular dynamics simulations

We studied the effect of using different heme group charge parametrization methods and schemes (Merz-Kollman, CHELPG, and single- and multiconformational RESP) on the quality of the results produced by the constant-(pH,E) MD method, applied to the redox titration of Desulfovibrio vulgaris Hildenborough cytochrome c(3). These new and more accurate charge sets enabled us to overcome the previously reported dependence of the method's performance on the dielectric constant, ε, assigned to the protein region. In particular, we found a systematic, clear shift of the E(mod) toward more negative values than those previously reported, in agreement with an electrostatics based reasoning. The simulations showed strong coupling between protonating/redox sites. We were also able to capture significant direct and, especially, indirect interactions between hemes, such as those mediated by histidine 67. Our results highlight the importance of having a good quantum description of the system before deriving atomic partial charges for classic force fields.

Constant pH molecular dynamics (CpHMD) and molecular docking studies of CquiOBP1 pH-induced ligand releasing mechanism

The odorant binding protein of Culex quinquefasciatus (CquiOBP1), expressed on the insect antenna, is crucial for the investigation of trapping baited with oviposition semi-chemicals and controlling mosquito populations. The acidic titratable residues pKa prediction and the ligand binding poses investigation in two systems (pH 7 and pH 5) are studied by constant pH molecular dynamics (CpHMD) and molecular docking methods. Research results reveal that the change of the protonation states would disrupt some important H-bonds, such as Asp 66-Asp 70, Glu 105-Asn 102, etc. The cleavage of these H-bonds leads to the movement of the relative position of hydrophobic tunnel, N- and C- termini loops and pH-sensing triad (His23-Tyr54-Val125) in acid solution. Ligand MOP has lower affinity and shows different binding poses to protein CquiOBP1 at pH 5. This ligand may be released from another tunnel between helices α3 and α4 in acidic environment. However, it would bind to the protein with high affinity in neutral environment. This work could provide more penetrating understanding of the pH-induced ligand-releasing mechanism.

pH-Dependent conformational changes in proteins and their effect on experimental pK(a)s: the case of Nitrophorin 4

The acid-base behavior of amino acids is an important subject of study due to their prominent role in enzyme catalysis, substrate binding and protein structure. Due to interactions with the protein environment, their pK_as can be shifted from their solution values and, if a protein has two stable conformations, it is possible for a residue to have different “microscopic”, conformation-dependent pK_a values. In those cases, interpretation of experimental measurements of the pK_a is complicated by the coupling between pH, protonation state and protein conformation. We explored these issues using Nitrophorin 4 (NP4), a protein that releases NO in a pH sensitive manner. At pH 5.5 NP4 is in a closed conformation where NO is tightly bound, while at pH 7.5 Asp30 becomes deprotonated, causing the conformation to change to an open state from which NO can easily escape. Using constant pH molecular dynamics we found two distinct microscopic Asp30 pK_as: 8.5 in the closed structure and 4.3 in the open structure. Using a four-state model, we then related the obtained microscopic values to the experimentally observed “apparent” pK_a, obtaining a value of 6.5, in excellent agreement with experimental data. This value must be interpreted as the pH at which the closed to open population transition takes place. More generally, our results show that it is possible to relate microscopic structure dependent pKa values to experimentally observed ensemble dependent apparent pK_as and that the insight gained in the relatively simple case of NP4 can be useful in several more complex cases involving a pH dependent transition, of great biochemical interest.

The interaction of an amino acid with its protein environment can result in an acid-base behavior that is very different from what would be observed in solution. This environment can be greatly altered when the protein changes conformation. As a result, the amino acid will have two different “microscopic” pK_a values. Nitrophorin 4 is a good case study to explore this behavior, because it undergoes a pH-dependent conformational change that is well characterized experimentally. Using computer simulation tools, we found that the key titratable Aspartic acid 30, has two very different microscopic pK_as: 4.3 and 8.5, which are significantly different to the observed transition pK_a in solution. However, using a simple model, we were able to understand how this causes the conformational change to take place at pH∼6.5, as measured experimentally. The insight gained in this relatively simple case can be useful in other more complex cases where the apparent pK_a is also a result of the interplay of different conformations where some amino acids experience very different environments.

pH-replica exchange molecular dynamics in proteins using a discrete protonation method

Protonation equilibria in biological molecules modulates structure, dynamics, and function. A pH-replica exchange molecular dynamics (pH-REMD) method is described here to improve the coupling between conformational and protonation sampling. Under a Hamiltonian replica exchange setup, conformations are swapped between two neighboring replicas, which themselves are at different pHs. The method has been validated on a series of biological systems. We applied pH-REMD to a series of model compounds, to an terminally charged ADFDA pentapeptide, and to a heptapeptide derived from the ovomucoid third domain (OMTKY3). In all of those systems, the predicted pK(a) by pH-REMD is very close to the experimental value and almost identical to the ones obtained by constant pH molecular dynamics (CpH MD). The method presented here, pH-REMD, has the advantage of faster convergence properties due to enhanced sampling of both conformation and protonation spaces.

Constant pH replica exchange molecular dynamics in biomolecules using a discrete protonation model

A constant pH replica exchange molecular dynamics (REMD) method is proposed and implemented to improve coupled protonation and conformational state sampling. By mixing conformational sampling at constant pH (with discrete protonation states) with a temperature ladder, this method avoids conformational trapping. Our method was tested and applied to seven different biological systems. The constant pH REMD not only predicted pKa correctly for small, model compounds but also converged faster than constant pH molecular dynamics (MD). We further tested our constant pH REMD on a heptapeptide from ovomucoid third domain (OMTKY3). Although constant pH REMD and MD produced very close pKa values, the constant pH REMD showed its advantage in the efficiency of conformational and protonation state samplings.

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.

Protein-Folding Dynamics: Overview of Molecular Simulation Techniques

Molecular dynamics (MD) is an invaluable tool with which to study protein folding in silico. Although just a few years ago the dynamic behavior of a protein molecule could be simulated only in the neighborhood of the experimental conformation (or protein unfolding could be simulated at high temperature), the advent of distributed computing, new techniques such as replica-exchange MD, new approaches (based on, e.g., the stochastic difference equation), and physics-based reduced models of proteins now make it possible to study protein-folding pathways from completely unfolded structures. In this review, we present algorithms for MD and their extensions and applications to protein-folding studies, using all-atom models with explicit and implicit solvent as well as reduced models of polypeptide chains.

The pH-dependent conformational states of kyotorphin: a constant-pH molecular dynamics study

An extensive conformational study of the analgesic dipeptide kyotorphin (L-Tyr-L-Arg) at different pH values was performed using a constant-pH molecular dynamics method. This dipeptide showed a remarkable pH-dependent conformational variety. The protonation of the N-terminal amine was identified as a key element in the transition between the more extended and the more packed conformational states, as monitored by the dihedral angle defined by the atoms 1Cbeta-1Calpha-2Calpha-2Cbeta. The principal-component analysis of kyotorphin identified two major conformational populations (the extended trans and the packed cis) together with conformations that occur exclusively at extreme pH values. Other, less stable conformations were also identified, which help us to understand the transitions between the predominant populations. The fitting of kyotorphin's conformational space to the structure of morphine resulted in a set of conformers that were able to fulfill most of the constraints for the mu-receptor. These results suggest that there may be strong similarities between the kyotorphin receptor and the structural family of opioid receptors.

Assessment of two theoretical methods to estimate potentiometric titration curves of peptides: comparison with experiment

We compared the ability of two theoretical methods of pH-dependent conformational calculations to reproduce experimental potentiometric titration curves of two models of peptides: Ac-K5-NHMe in 95% methanol (MeOH)/5% water mixture and Ac-XX(A)7OO-NH2 (XAO) (where X is diaminobutyric acid, A is alanine, and O is ornithine) in water, methanol (MeOH), and dimethyl sulfoxide (DMSO), respectively. The titration curve of the former was taken from the literature, and the curve of the latter was determined in this work. The first theoretical method involves a conformational search using the electrostatically driven Monte Carlo (EDMC) method with a low-cost energy function (ECEPP/3 plus the SRFOPT surface-solvation model, assumming that all titratable groups are uncharged) and subsequent reevaluation of the free energy at a given pH with the Poisson-Boltzmann equation, considering variable protonation states. In the second procedure, molecular dynamics (MD) simulations are run with the AMBER force field and the generalized Born model of electrostatic solvation, and the protonation states are sampled during constant-pH MD runs. In all three solvents, the first pKa of XAO is strongly downshifted compared to the value for the reference compounds (ethylamine and propylamine, respectively); the water and methanol curves have one, and the DMSO curve has two jumps characteristic of remarkable differences in the dissociation constants of acidic groups. The predicted titration curves of Ac-K5-NHMe are in good agreement with the experimental ones; better agreement is achieved with the MD-based method. The titration curves of XAO in methanol and DMSO, calculated using the MD-based approach, trace the shape of the experimental curves, reproducing the pH jump, while those calculated with the EDMC-based approach and the titration curve in water calculated using the MD-based approach have smooth shapes characteristic of the titration of weak multifunctional acids with small differences between the dissociation constants. Nevertheless, quantitative agreement between theoretically predicted and experimental titration curves is not achieved in all three solvents even with the MD-based approach, which is manifested by a smaller pH range of the calculated titration curves with respect to the experimental curves. The poorer agreement obtained for water than for the nonaqueous solvents suggests a significant role of specific solvation in water, which cannot be accounted for by the mean-field solvation models.

Constant-pH Molecular Dynamics with Ionic Strength Effects: Protonation−Conformation Coupling in Decalysine

A new implementation of the stochastic titration method for constant-pH molecular dynamics is presented, which introduces ionic strength effects in the simulations. In addition, the new implementation uses a faster molecular dynamics algorithm and an improved treatment of protonation events and of their effect on force field parameters. This new methodology is applied to a decalysine peptide, yielding very good quantitative agreement with experiments, both in terms of titration and helix-coil transition. The results show a significant dependence on ionic strength, illustrating the importance of including this parameter in constant-pH molecular dynamics simulations. Overall, the method seems to properly capture the protonation-conformation coupling and its dependence on ionic strength.

Constant pH molecular dynamics in generalized Born implicit solvent

A new method is proposed for constant pH molecular dynamics (MD), employing generalized Born (GB) electrostatics. Protonation states are modeled with different charge sets, and titrating residues sample a Boltzmann distribution of protonation states as the simulation progresses, using Monte Carlo sampling based on GB-derived energies. The method is applied to four different crystal structures of hen egg-white lysozyme (HEWL). pK(a) predictions derived from the simulations have root-mean-square (RMS) error of 0.82 relative to experimental values. Similarity of results between the four crystal structures shows the method to be independent of starting crystal structure; this is in contrast to most electrostatics-only models. A strong correlation between conformation and protonation state is noted and quantitatively analyzed, emphasizing the importance of sampling protonation states in conjunction with dynamics.

Gibbs adsorption isotherm combined with Monte Carlo sampling to see action of cosolutes on protein folding

Driven by conditions set by smaller solutes, proteins fold and unfold. Experimentally, these conditions are stated as intensive variables--pH and other chemical potentials--as though small solutes were infinite resources that come at an externally varied free energy cost. Computationally, the finite spaces of simulation allow only fixed numbers of these solutes. By combining the analytic Gibbs adsorption isotherm with the computational Monte Carlo sampling of polymer configurations, we have been able to overcome an inherent limitation of computer simulation. The idea is to compute analytically the free energy changes wrought by solutes on each particular configuration. Then numerical computation is needed only to sample the set of configurations as efficiently as when no bathing solute is present. For illustration, the procedure is applied to an idealized two-dimensional heteropolymer to yield lessons about the effect of cosolutes on protein stability.

Conformational sampling for the impatient

Several new methods for sampling conformations of biomolecules have appeared recently. A brief review thereof is presented, with particular emphasis on applications that have been published, and suitability for different kinds of systems. Four methods (namely: RESPA, replica-exchange molecular dynamics, CONCOORD and Gaussian network method) are readily applicable for biomolecular systems.

Constant-pH molecular dynamics study of protonation-structure relationship in a heptapeptide derived from ovomucoid third domain

Molecular dynamics (MD) simulations with implicit solvent and variable protonation states for titratable residues at constant pH are performed for a short peptide derived from ovomucoid third domain (OMTKY3), acetyl-Ser-Asp-Asn-Lys-Thr-Tyr-Gly-amide (residues 26-32 of OMTKY3). Nuclear magnetic resonance (NMR) measurements indicate that the pK(a) for Asp is 3.6. However, if the charge on Lys is neutralized by acetylation, then the pK(a) for Asp is 4.0. These pK(a)'s, and therefore the Asp-Lys interaction, are insensitive to changes in ionic strength. The constant-pH MD simulations for both variants of the heptapeptide yield Asp pK(a) values that are 0.6-0.9 pH units greater than experimental values, but the difference between the variants that is observed in the NMR experiments is reproduced much better. Moreover, the simulations suggest that Asp-Lys interactions do not dominate the behavior of this heptapeptide, even for normal Lys residue where there is a possibility of forming a salt bridge between negatively charged Asp and positively charged Lys. This is consistent with the experimentally observed independence of Asp pK(a) values with respect to ionic strength. Another important result of the simulations with variable protonation states is that they lead to ensembles of the heptapeptide structures that are different from those derived from simulations with fixed protonation states. It should be stressed that these results are for structures generated entirely by computer simulations without any restrictions imposed by experimental data.