Constant pH Molecular Dynamics Simulations Research Articles

Thermodynamic Coupling of Protonation and Conformational Equilibria in Proteins: Theory and Simulation

Ionization-coupled conformational phenomena are ubiquitous in biology. However, quantitative characterization of the underlying thermodynamic cycle comprised of protonation and conformational equilibria has remained an elusive goal. Here we use theory and continuous constant pH molecular dynamics (CpHMD) simulations to provide a thermodynamic description for the coupling of proton titration and conformational exchange between two distinct states of a protein. CpHMD simulations with a hybrid-solvent scheme and the pH-based replica-exchange (REX) protocol are applied to obtain the equilibrium constants and atomic details of the ionization-coupled conformational exchange between open and closed states of an engineered mutant of staphylococcal nuclease. Although the coupling of protonation and conformational equilibria is not exact in the simulation, the results are encouraging. They demonstrate that REX-CpHMD simulations can be used to study thermodynamics of ionization-coupled conformational processes—which has not possible using present experimental techniques or traditional simulations based on fixed protonation states.

Structural and Functional Difference of Pheromone Binding Proteins in Discriminating Chemicals in the Gypsy Moth, Lymantria Dispar

Pheromone-binding proteins (PBPs) of the gypsy moth, Lymantria dispar L., play an important role in olfaction. Here structures of PBPs were first built by Homology Modeling, and each model of PBPs had seven α-helices and a large hydrophobic cavity including 25 residues for PBP1 and 30 residues for PBP2. Three potential semiochemicals were first screened by CDOCKER program based on the PBP models and chemical database. These chemicals were Palmitic acid n-butyl ester (Pal), Bis(3,4-epoxycyclohexylmethyl) adipate (Bis), L-trans-epoxysuccinyl-isoleucyl-proline methyl ester propylamide (CA-074). The analysis of chemicals docking the proteins showed one hydrogen bond was established between the residues Lys94 and (+)-Disparlure ((+)-D), and л-л interactions were present between Phe36 of PBP1 and (+)-D. The Lys94 of PBP1 formed two and three hydrogen bonds with Bis and CA-074, respectively. There was no residue of PBP2 interacting with these four chemicals except Bis forming one hydrogen bond with Lys121. After simulating the conformational changes of LdisPBPs at pH7.3 and 5.5 by constant pH molecular dynamics simulation in implicit solvent, the N-terminal sequences of PBPs was unfolded, only having five α-helices, and PBP2 had larger binding pocket at 7.3 than PBP1. To investigate the changes of α-helices at different pH, far-UV and near-UV circular dichroism showed PBPs consist of α-helices, and the tertiary structures of PBP1 and PBP2 were influenced at pH7.3 and 5.5. The fluorescence binding assay indicated that PBP1 and PBP2 have similarly binding affinity to (+)-D at pH 5.5 and 7.3, respectively. At pH 5.5, the dissociation constant of the complex between PBP1 and 2-decyl-1-oxaspiro [2.2] pentane (OXP1) was 0.68±0.01μM, for (+)-D was 5.32±0.11μM, while PBP2 with OXP1 and (+)-D were 1.88±0.02μM and 5.54±0.04μM, respectively. Three chemicals screened had higher affinity to PBP1 than (+)-D except Pal at pH5.5, and had lower affinity than (+)-D at pH7.3. To PBP2, these chemicals had lower affinity than the sex pheromone except Bis at pH 5.5 and pH 7.3. Only PBP1 had higher affinity with Sal than the sex pheromone at pH 5.5. Therefore, the structures of PBP1 and PBP2 had different changes at pH5.5 and 7.3, showing different affinity to chemicals. This study helps understanding the role of PBPs as well as in developing more efficient chemicals for pest control.

Probing pH-dependent dissociation of HdeA dimers.

HdeA protein is a small, ATP-independent, acid stress chaperone that undergoes a dimer-to-monomer transition in acidic environments. The HdeA monomer binds a broad range of proteins to prevent their acid-induced aggregation. To understand better HdeA's function and mechanism, we perform constant-pH molecular dynamics simulations (CPHMD) to elucidate the details of the HdeA dimer dissociation process. First the pK(a) values of all the acidic titratable groups in HdeA are obtained and reveal a large pK(a) shift only for Glu(37). However, the pH-dependent monomer charge exhibits a large shift from -4 at pH > 6 to +6 at pH = 2.5, suggesting that the dramatic change in charge on each monomer may drive dissociation. By combining the CPHMD approach with umbrella sampling, we demonstrate a significant stability decrease of the HdeA dimer when the environmental pH changes from 4.0 to 3.5 and identify the key acidic residue-lysine interactions responsible for the observed pH sensing in HdeA chaperon activity function.

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.

Histidine-Mediated pH-Sensitive Regulation of M-Ficolin:GlcNAc Binding Activity in Innate Immunity Examined by Molecular Dynamics Simulations

BackgroundM-ficolin, a pathogen recognition molecule in the innate immune system, binds sugar residues including N-acetyl-D-glucosamine (GlcNAc), which is displayed on invading microbes and on apoptotic cells. The cis and trans Asp282-Cys283 peptide bond in the M-ficolin, which was found to occur at neutral and acidic pH in crystal structures, has been suggested to represent binding and non-binding activity, respectively. A detailed understanding of the pH-dependent conformational changes in M-ficolin and pH-mediated discrimination mechanism of GlcNAc-binding activity are crucial to both immune-surveillance and clearance of apoptotic cells.Methodology/Principal FindingsBy immunodetection analysis, we found that the pH-sensitive binding of GlcNAc is regulated by a conformational equilibrium between the active and inactive states of M-ficolin. We performed constant pH molecular dynamics (MD) simulation at a series of pH values to explore the pH effect on the cis-trans isomerization of the Asp282-Cys283 peptide bond in the M-ficolin fibrinogen-like domain (FBG). Analysis of the hydrogen bond occupancy of wild type FBG compared with three His mutants (H251A, H284A and H297A) corroborates that His284 is indispensible for pH-dependent binding. H251A formed new but weaker hydrogen bonds with GlcNAc. His297, unlike the other two His mutants, is more dependent on the solution pH and also contributes to cis-trans isomerization of the Asp282-Cys283 peptide bond in weak basic solution.Conclusions/SignificanceConstant pH MD simulation indicated that the cis active isomer of Asp282-Cys283 peptide bond was predominant around neutral pH while the trans bond gradually prevailed towards acidic environment. The protonation of His284 was found to be associated with the trans-to-cis isomerization of Asp282-Cys283 peptide bond which dominantly regulates the GlcNAc binding. Our MD simulation approach provides an insight into the pH-sensitive proteins and hence, ligand binding activity.

Molecular dynamics method for proteins with ionization-conformation coupling and equilibrium titration

A new realization of the constant-pH molecular dynamics simulation method is proposed. Molecular dynamics simulation is performed for a protein in the most probable ionization microstate of the current conformation taking into account the potential of mean force of protein molecule in the water-proton bath in equilibrium titration conditions (MD-pH-ET). It is shown that: 1) the optimal one is the simulation of the protein in the most probable ionization for a given conformation, taking into account the correction ionization potential of mean force, which results in the librium ensemble of ionization states; 2) new method MD-pH-ET allow one to carry out an optimization of protein structure and the total free energy of a protein in the aqueous solution at constant pH, and to calculate the pH-dependent properties. Method MD-pH-ET possesses the unique features: 1) it uses the most precise and computational-effective realization of calculation of the electrostatic energy of a protein in water solution, the model of continuous dielectric media with Poisson equation and the generalized Born method with "ideal" Born atomic radii; 2) it uses the same model of the potential energy surface in the ionization-conformational phase space, both for calculating the potential energy of the protein and atomic forces and for determining the most probable ionization states; 3) it calculates the total free energy of the protein in the aqueous solution in proton reservoir under the conditions of equilibrium titration. The workability of the new method MD-pH-ET is demonstrated for the molecule of protein BPTI.

Atomistic Simulation and Measurement of pH Dependent Cancer Therapeutic Interactions with Nanodiamond Carrier

In this work, we have combined constant-pH molecular dynamics simulations and experiments to provide a quantitative analysis of pH dependent interactions between doxorubicin hydrochloride (DOX) cancer therapeutic and faceted nanodiamond (ND) nanoparticle carriers. Our study suggests that when a mixture of faceted ND and DOX is dissolved in a solvent, the pH of this solvent plays a controlling role in the adsorption of DOX molecules on the ND. We find that the binding of DOX molecules on ND occurs only at high pH and requires at least ∼10% of ND surface area to be fully titrated for binding to occur. As such, this study reveals important mechanistic insight underlying an ND-based pH-controlled therapeutic platform.

Peptide Isomerization Induced by pH Change Regulates the S1 Binding Site in Ficolins

Infection-inflammation mediated interactions between human ficolin and the pathogen GlcNAc is associated with local acidosis, leading to antimicrobial action. Therefore, revealing the precise molecular conformation induced by pH-shift is crucial in understanding the immune response. Here, we performed constant-pH molecular dynamics simulations on the L-ficolin fibrinogen-like domain over pH 4.5–9. An unusual cis-Asn244-Cys245 peptide bond prevailed over the pH range in the S1 binding site. Analysis of the hydrogen-bond network at S1 suggested Asn244 to be indispensible for maintaining the cis form of Asn244-Cys245, and the absence of the hydroxyl group on Phe262 accounts for the lack of GlcNAc binding.

Constant-pH Molecular Dynamics Simulations Reveal a β-Rich Form of the Human Prion Protein

The misfolding of the prion protein (PrP) into a pathogenic β-rich form (PrP(Sc)) has been suggested to occur in the endocytic pathway, triggered by low pH. In this work we performed several constant-pH molecular dynamics simulations of human PrP 90-231 in the pH range 2-7, totaling more than 2 μs. We observed a strong conformational pH dependence where on average the helix content decreased and the β content increased toward acidic pH. Unlike some proposed models, the flexible N-terminus region did not gain stable structure at low pH. Rather, the main structural changes occurred on the helix-rich C-terminus core, as proposed in other models, namely, in the regions around 135-155 and 185-200. The protonation of His187 is found to be associated with a loss of interaction between two PrP subdomains, potentially playing a major role in the misfolding process. In one of the simulations at pH 2, a stable β-rich structure was formed that may be an intermediate of PrP(Sc) formation, indicating that misfolding may precede dimerization.

Simulations of Surfactant and Lipid Assemblies with Generalized Born Implicit Solvent Models

In recent years, all-atom and coarse-grained models have been developed for theoretical studies of surfactant and lipid assemblies. Here we describe the development of hybrid models incorporating atomic-level representation of surfactant molecules and continuum description for solvent. Specifically, we show that atomistic simulations of neutral, anionic and cationic surfactants can be performed using the generalized Born implicit solvent model by careful parameterization of atomic input radii. We applied the new models to constant pH molecular dynamics simulations of acid-base titration of fatty acid solubilized in dodecyl triehyleneglycol either (DE3), dodecylsulfaphate (SDS) and dodecyl tetramethylammonium (DTA) micelles. We show that simulation results are able to reproduce the experimental data and offer atomically detailed explanation for the abnormal titration behavior of fatty acid in cationic micelles. The combined atomistic solute and continuum solvent models significantly reduce computational and represent an attractive alternative approach for theoretical studies of interfacial phenomena involving surfactant and lipid assemblies.

Cofilin is a pH sensor for actin free barbed end formation: role of phosphoinositide binding

Newly generated actin free barbed ends at the front of motile cells provide sites for actin filament assembly driving membrane protrusion. Growth factors induce a rapid biphasic increase in actin free barbed ends, and we found both phases absent in fibroblasts lacking H+ efflux by the Na-H exchanger NHE1. The first phase is restored by expression of mutant cofilin-H133A but not unphosphorylated cofilin-S3A. Constant pH molecular dynamics simulations and nuclear magnetic resonance (NMR) reveal pH-sensitive structural changes in the cofilin C-terminal filamentous actin binding site dependent on His133. However, cofilin-H133A retains pH-sensitive changes in NMR spectra and severing activity in vitro, which suggests that it has a more complex behavior in cells. Cofilin activity is inhibited by phosphoinositide binding, and we found that phosphoinositide binding is pH-dependent for wild-type cofilin, with decreased binding at a higher pH. In contrast, phosphoinositide binding by cofilin-H133A is attenuated and pH insensitive. These data suggest a molecular mechanism whereby cofilin acts as a pH sensor to mediate a pH-dependent actin filament dynamics.

Prediction of Protein Solubility from Calculation of Transfer Free Energy

Solubility plays a major role in protein purification, and has serious implications in many diseases. We studied the effects of pH and mutations on protein solubility by calculating the transfer free energy from the condensed phase to the solution phase. The condensed phase was modeled as an implicit solvent, with a dielectric constant lower than that of water. To account for the effects of pH, the protonation states of titratable side chains were sampled by running constant-pH molecular dynamics simulations. Conformations were then selected for calculations of the electrostatic solvation energy: once for the condensed phase, and once for the solution phase. The average transfer free energy from the condensed phase to the solution phase was found to predict reasonably well the variations in solubility of ribonuclease Sa and insulin with pH. This treatment of electrostatic contributions combined with a similar approach for nonelectrostatic contributions led to a quantitative rationalization of the effects of point mutations on the solubility of ribonuclease Sa. This study provides valuable insights into the physical basis of protein solubility.

Evidence for a Catalytic Dyad in the Active Site of Homocitrate Synthase from Saccharomyces cerevisiae

Homocitrate synthase (acetyl-coenzyme A: 2-ketoglutarate C-transferase; E.C. 2.3.3.14) (HCS) catalyzes the condensation of acetyl-CoA (AcCoA) and alpha-ketoglutarate (alpha-KG) to give homocitrate and CoA. Although the structure of an HCS has not been solved, the structure of isopropylmalate synthase (IPMS), a homologue, has been solved (Koon, N., Squire, C. J., and Baker, E. N. (2004) Proc. Natl. Acad. Sci. U.S.A. 101, 8295-8300). Three active site residues in IPMS, Glu-218, His-379, and Tyr-410, were proposed as candidates for catalytic residues involved in deprotonation of the methyl group of AcCoA prior to the Claisen condensation to give homocitrylCoA. All three of the active site residues in IPMS are conserved in the HCS from Saccharomyces cerevisiae. Site-directed mutagenesis has been carried out to probe the role of the homologous residues, Glu-155, His-309, and Tyr-320, in the S. cerevisiae HCS. No detectable activity was observed for the H309A and H309N mutant enzyme, but a slight increase in activity was observed for H309A in the presence of 300 mM imidazole, which is still 1000-fold lower than that of wild type (wt). The E155Q and E155A mutant enzymes exhibited 1000-fold lower activity than wt. The activity of E155A, but not of E155Q, could be partially rescued by formate; a K act of 60 mM with a modest 4-fold maximum activation was observed. In the presence of formate, E155A gives k cat, K AcCoA, and K alpha-KG values of 0.0031 s (-1), 13 muM, and 39 microM, respectively, while a primary kinetic deuterium isotope effect of about 1.4 was obtained on V, with deuterium in the methyl of AcCoA. The pH dependence of k cat for E155A in the presence of formate gave a p K a of 7.9 for a group that must be protonated for optimum activity, similar to that observed for the wt enzyme. However, a partial change was observed on the acid side of the profile, compared to the all or none change observed for wt giving a p K a of about 6.7. The k cat for E155Q decreased at high pH, similar to the wt enzyme, but was pH independent at low pH. The Y320F mutant enzyme only lost 25-fold activity compared to that of the wt, giving k cat, K AcCoA, and K alpha-KG values of 0.039 s (-1), 33 microM, and 140 microM, respectively, and a primary kinetic deuterium isotope effect of 1.3 and 1.8 on V/ K AcCoA and V, respectively; the pH dependence of k cat was similar to that of the wt. These data, combined with a constant pH molecular dynamics simulation study, suggest that a catalytic dyad comprising Glu-155 and His-309 acts to deprotonate the methyl group of AcCoA, while Tyr320 is likely not directly involved in catalysis, but may aid in orienting the reactant and/or the catalytic dyad.

Calculating pH and Salt Dependence of Protein-Protein Binding

Ionic strength- (or salt-) effects on the protein-protein binding free energy has been included in many computational studies, while comparatively fewer computational studies have looked at the corresponding effect of pH. The pH dependence can be very complex if several groups change protonation state, while the ionic strength dependence usually scales as ln(I), and the main challenge is to predict the magnitude of the correlation. However, there is now very strong indication that pH effects due to binding induced changes in protonation states make a non-negligible contribution to the binding energy of most protein-protein complexes. This observation, together with more efficient pKa prediction methods and the emergence of constant pH molecular dynamics simulations to model the protonation-dependent structural changes will spark more experimental and theoretical work in pH effects on protein-protein binding.

Linking folding with aggregation in Alzheimer's β-amyloid peptides

Growing evidence suggests that the beta-amyloid (Abeta) peptides of Alzheimer's disease are generated in early endosomes and that small oligomers are the principal toxic species. We sought to understand whether and how the solution pH, which is more acidic in endosomes than the extracellular environment, affects the conformational processes of Abeta. Using constant pH molecular dynamics simulations of two model peptides, Abeta(1-28) and Abeta(10-42), we found that the folding landscape of Abeta is strongly modulated by pH and is most favorable for hydrophobically driven aggregation at pH 6. Thus, our theoretical findings substantiate the possibility that Abeta oligomers develop intracellularly before secretion into the extracellular milieu, where they may disrupt synaptic activity or act as seeds for plaque formation.

Simulating botulinum neurotoxin with constant pH molecular dynamics in Generalized Born implicit solvent

A new method was proposed by Mongan et al. for constant pH molecular dynamics simulation and was implemented in AMBER 8 package. Protonation states are modeled with different charge sets, and titrating residues are sampled from a Boltzmann distribution of protonation states. The simulation periodically adopts Monte Carlo sampling based on Generalized Born (GB) derived energies. However, when this approach was applied to a bio-toxin, Botulinum Neurotoxin Type A (BoNT/A) at pH 4.4, 4.7, 5.0, 6.8 and 7.2, the p K a predictions yielded by the method were inconsistent with the experimental values. The systems being simulated were divergent. Furthermore, the system behaviors in a very weak acidic solution (pH 6.8) and in a very weak basic solution (pH 7.2) were significantly different from the neutral case (pH 7.0). Hence, we speculate this method may require further study for modeling large biomolecule.

Folding Intermediate in the Villin Headpiece Domain Arises from Disruption of a N-Terminal Hydrogen-Bonded Network

The critical role of partially folded intermediates in protein misfolding and amyloid formation has sparked interest in exploring factors that control the formation of these metastable species. Recent NMR experiments reported a sparsely populated intermediate of the villin headpiece domain, in which the N-terminal subdomain is random coil-like under native conditions. Here we use continuous constant pH molecular dynamics simulations with replica-exchange sampling protocol to test the hypothesis that the conformational state obtained by simulations derived from a solution NMR structure represents a putative intermediate in which the N-terminal subdomain is partially unfolded. On the basis of the unique titration behavior of His41 as well as the structural and dynamic properties of this state, we propose that the loss of a hydrogen bond between Nδ of His41 and the backbone carbonyl oxygen of E14 in the NMR structure leads to the loss of inter-subdomain contacts as well as partial disruption of the N-terminal hydrophobic cluster. Thus, the loss of the hydrogen-bonded network and not the protonation of His41 is a prerequisite for unfolding.

Toward the Accurate First-Principles Prediction of Ionization Equilibria in Proteins

The calculation of pK(a) values for ionizable sites in proteins has been traditionally based on numerical solutions of the Poisson-Boltzmann equation carried out using a high-resolution protein structure. In this paper, we present a method based on continuous constant pH molecular dynamics (CPHMD) simulations, which allows the first-principles description of protein ionization equilibria. Our method utilizes an improved generalized Born implicit solvent model with an approximate Debye-Hückel screening function to account for salt effects and the replica-exchange (REX) protocol for enhanced conformational and protonation state sampling. The accuracy and robustness of the present method are demonstrated by 1 ns REX-CPHMD titration simulations of 10 proteins, which exhibit anomalously large pK(a) shifts for the carboxylate and histidine side chains. The experimental pK(a) values of these proteins are reliably reproduced with a root-mean-square error ranging from 0.6 unit for proteins containing few buried ionizable side chains to 1.0 unit or slightly higher for proteins containing ionizable side chains deeply buried in the core and experiencing strong charge-charge interactions. This unprecedented level of agreement with experimental benchmarks for the de novo calculation of pK(a) values suggests that the CPHMD method is maturing into a practical tool for the quantitative prediction of protein ionization equilibria, and this, in turn, opens a door to atomistic simulations of a wide variety of pH-coupled conformational phenomena in biological macromolecules such as protein folding or misfolding, aggregation, ligand binding, membrane interaction, and catalysis.

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 with Proton Tautomerism

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 12 and His 48 and the conformational states of His 48 and His 119 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.