Constant pH Molecular Dynamics Simulations Research Articles

How Reactive are Druggable Cysteines in Protein Kinases?

Targeted covalent inhibitors (TCIs) have been successfully developed as high-affinity and selective inhibitors of enzymes of the protein kinase family. These drugs typically act by undergoing an electrophilic addition with an active-site cysteine residue, so design of a TCI begins with the identification of a "druggable" cysteine. These electrophilic additions generally require deprotonation of the thiol to form a reactive anionic thiolate, so the acidity of the residue is a critical factor. Few experimental measurements of the p Ka's of druggable cysteines have been reported, so computational prediction could prove to be very important in selecting reactive cysteine targets. Here we report the computed p Ka's of druggable cysteines in selected protein kinases that are of clinical relevance for targeted therapies. The p Ka's of the cysteines were calculated using advanced computational methods based on all-atom replica-exchange thermodynamic integration molecular dynamics simulations in explicit solvent. We found that the acidities of druggable cysteines within protein kinases are diverse and elevated, indicating enormous differences in their reactivity. Constant-pH molecular dynamics simulations were also performed on selected protein kinases, and the results confirmed this varied range in the acidities of druggable cysteines. Many of these active-site cysteines have low exposure to solvent molecules, elevating their p Ka values. Electrostatic interactions with nearby anionic residues also elevate the p Ka's of cysteine residues in the active site. The results suggest that some cysteine residues within kinase binding sites will be slow to react with a TCI because of their low acidity. Several oncogenic kinase mutations were also modeled and found to have p Ka's similar to that of the wild-type kinase.

Understanding the pH-Dependent Reaction Mechanism of a Glycoside Hydrolase Using High-Resolution X-ray and Neutron Crystallography

Glycoside hydrolases (GHs) commonly use the retaining or inverting mechanisms to hydrolyze carbohydrates, and the rates of catalysis are usually pH dependent. Deeper understanding of these pH-dependent reaction mechanisms is of great importance for protein engineering and drug design. We used high-resolution X-ray crystallography to analyze the sugar ring configurations of an oligosaccharide ligand during hydrolysis for the family 11 GH, and the results support the 1S3 → 4H3 → 4C1 conformational itinerary. These results indicate that sugar ring flexibility may help to distort and break the glycosidic bond. Constant pH molecular dynamics simulations and neutron crystallography demonstrate that the catalytic glutamate residue (E177) has alternate conformational changes to transfer a proton to cleave the glycosidic bond. Furthermore, a neutron crystallography analysis shows that the H-bond length between E177 and its nearby tyrosine residue (Y88) is shortened when the pH increases, preventing E177 from rotat...

Membrane-Induced p Ka Shifts in wt-pHLIP and Its L16H Variant.

The pH (low) insertion peptides (pHLIPs) is a family of peptides that are able to insert into a lipid bilayer at acidic pH. The molecular mechanism of pHLIPs insertion, folding, and stability in the membrane at low pH is based on multiple protonation events, which are challenging to study at the molecular level. More specifically, the relation between the experimental p K of insertion (p Kexp) of pHLIPs and the p Ka of the key residues is yet to be clarified. We carried out a computational study, complemented with new experimental data, and established the influence of (de)protonation of titrable residues on the stability of the peptide membrane-inserted state. Constant-pH molecular dynamics simulations were employed to calculate the p Ka values of these residues along the membrane normal. In the wt-pHLIP, we identified Asp14 as the key residue for the stability of the membrane-inserted state, and its p Ka value is strongly correlated with the experimental p Kexp measured in thermodynamics studies. Also, in order to narrow down the pH range at which pHLIP is stable in the membrane, we designed a new pHLIP variant, L16H, where Leu in the 16th position was replaced by a titrable His residue. Our results showed that the L16H variant undergoes two transitions. The calculated p Ka and experimentally observed p Kexp values are in good agreement. Two distinct p Kexp values delimit a pH range where the L16H peptide is stably inserted in the membrane, while, outside this range, the membrane-inserted state is destabilized and the peptide exits from the bilayer. pHLIP peptides have been successfully used to target cancer cells for the delivery of diagnostics and therapeutic agents to acidic tumors. The fine-tuning of the stability of the pHLIP inserted state and its restriction to a narrow well-defined pH range might allow the design of new peptides, able to discriminate between tissues with different extracellular pH values.

Surface Charge Effects on Fe(II) Sorption and Oxidation at (110) Goethite Surfaces

Iron(III) oxides and oxyhydroxides are among the most reactive minerals in the environment, with surfaces that become charged when immersed in water. The governing role of surface charge over interfacial processes such as metal sorption is well understood. However, its role in interfacial redox reactions, such as when metal sorption is coupled to interfacial electron transfer (ET), is not. This is mainly because surface charge affects not only the types and densities of surface complexes formed but also their respective driving forces for ET. An important case is Fe(II)-catalyzed recrystallization of Fe(III)-oxyhydroxides, in which Fe(II) sorption and interfacial ET are closely linked. In this study, we used replica-exchange constant-pH molecular dynamics simulations (Zarzycki, P.; Smith, D. M.; Rosso, K. M. J. Chem. Theory Comput. 2015, 11, 1715–1724) to calculate the distance-dependent electrostatic potential at charged (110) surfaces of goethite particles, assessing its effect on previously computed Fe...

Predicting Catalytic Proton Donors and Nucleophiles in Enzymes: How Adding Dynamics Helps Elucidate the Structure-Function Relationships.

Despite the relevance of understanding structure-function relationships, robust prediction of proton donors and nucleophiles in enzyme active sites remains challenging. Here we tested three types of state-of-the-art computational methods to calculate the p Ka's of the buried and hydrogen bonded catalytic dyads in five enzymes. We asked the question what determines the p Ka order, i.e., what makes a residue proton donor vs a nucleophile. The continuous constant pH molecular dynamics simulations captured the experimental p Ka orders and revealed that the negative nucleophile is stabilized by increased hydrogen bonding and solvent exposure as compared to the proton donor. Surprisingly, this simple trend is not apparent from crystal structures and the static structure-based calculations. While the generality of the findings awaits further testing via a larger set of data, they underscore the role of dynamics in bridging enzyme structures and functions.

Role of Counterions in Constant-pH Molecular Dynamics Simulations of PAMAM Dendrimers.

Electrostatic interactions play a pivotal role in the structure and mechanism of action of most biomolecules. There are several conceptually different methods to deal with electrostatics in molecular dynamics simulations. Ionic strength effects are usually introduced using such methodologies and can have a significant impact on the quality of the final conformation space obtained. We have previously shown that full system neutralization can lead to wrong lipidic phases in the 25% PA/PC bilayer (J. Chem. Theory Comput. 2014,10, 5483–5492). In this work, we investigate how two limit approaches to the ionic strength treatment (implicitly with GRF or using full system neutralization with either GRF or PME) can influence the conformational space of the second-generation PAMAM dendrimer. Constant-pH MD simulations were used to map PAMAM’s conformational space at its full pH range (from 2.5 to 12.5). Our simulations clearly captured the coupling between protonation and conformation in PAMAM. Interestingly, the dendrimer conformational distribution was almost independent of the ionic strength treatment methods, which is in contrast to what we have observed in charged lipid bilayers. Overall, our results confirm that both GRF with implicit ionic strength and a fully neutralized system with PME are valid approaches to model charged globular systems, using the GROMOS 54A7 force field.

Constant-pH Molecular Dynamics Simulations for Large Biomolecular Systems.

An increasingly important endeavor is to develop computational strategies that enable molecular dynamics (MD) simulations of biomolecular systems with spontaneous changes in protonation states under conditions of constant pH. The present work describes our efforts to implement the powerful constant-pH MD simulation method, based on a hybrid nonequilibrium MD/Monte Carlo (neMD/MC) technique within the highly scalable program NAMD. The constant-pH hybrid neMD/MC method has several appealing features; it samples the correct semigrand canonical ensemble rigorously, the computational cost increases linearly with the number of titratable sites, and it is applicable to explicit solvent simulations. The present implementation of the constant-pH hybrid neMD/MC in NAMD is designed to handle a wide range of biomolecular systems with no constraints on the choice of force field. Furthermore, the sampling efficiency can be adaptively improved on-the-fly by adjusting algorithmic parameters during the simulation. Illustrative examples emphasizing medium- and large-scale applications on next-generation supercomputing architectures are provided.

Surfactant pKa Calculations using Molecular Dynamics Simulations

The protonation state for ionizable surfactants is crucial for self-assembly, phase behavior, and morphology. For example, oleic acid forms vesicles at intermediate pH conditions (∼7-9), micelles at high pH, and oil phases at low pH. The pKa for oleic acid shifts from ∼4.8 as monomers in water to ∼7.5 in model phospholipid vesicles. It remains challenging to accurately model surfactant self-assembly and protonation behavior with molecular dynamics simulations. We use constant pH molecular dynamics simulations and thermodynamic integration calculations for protonating a single fatty acid in different chemical environments.

Replica exchange enveloping distribution sampling (RE-EDS): A robust method to estimate multiple free-energy differences from a single simulation.

In molecular dynamics (MD) simulations, free-energy differences are often calculated using free energy perturbation or thermodynamic integration (TI) methods. However, both techniques are only suited to calculate free-energy differences between two end states. Enveloping distribution sampling (EDS) presents an attractive alternative that allows to calculate multiple free-energy differences in a single simulation. In EDS, a reference state is simulated which "envelopes" the end states. The challenge of this methodology is the determination of optimal reference-state parameters to ensure equal sampling of all end states. Currently, the automatic determination of the reference-state parameters for multiple end states is an unsolved issue that limits the application of the methodology. To resolve this, we have generalised the replica-exchange EDS (RE-EDS) approach, introduced by Lee et al. [J. Chem. Theory Comput. 10, 2738 (2014)] for constant-pH MD simulations. By exchanging configurations between replicas with different reference-state parameters, the complexity of the parameter-choice problem can be substantially reduced. A new robust scheme to estimate the reference-state parameters from a short initial RE-EDS simulation with default parameters was developed, which allowed the calculation of 36 free-energy differences between nine small-molecule inhibitors of phenylethanolamine N-methyltransferase from a single simulation. The resulting free-energy differences were in excellent agreement with values obtained previously by TI and two-state EDS simulations.

Molecular dynamics simulation of chitinase I from Thermomyces lanuginosus SSBP to ensure optimal activity

The fungal chitinase I obtained from Thermomyces lanuginosus SSBP, a thermophilic deuteromycete, has an optimum growth temperature and pH of 323.15 K and 6.5, respectively. This enzyme plays an important task in the defence mechanism of organisms against chitin-containing parasites by hydrolysing β-1, 4-linkages in chitin. It acts as both anti-fungal and biofouling agents, with some being thermostable and suitable for the industrial applications. Three-dimensional model of chitinase I enzyme was predicted and analysed using various bioinformatics tools. The structure of chitinase I exhibited a well-defined TIM barrel topology with an eight-stranded α/β domain. Structural analysis and folding studies at temperatures ranging from 300 to 375 K using 10 ns molecular dynamics simulations clearly showed the stability of the protein was evenly distributed even at higher temperatures, in accordance with the experimental results. We also carried out a number of 20 ns constant pH molecular dynamics simulations of chitinase I at a pH range 2–6 in a solvent. This work was aimed at establishing the optimum activity and stability profiles of chitinase I. We observed a strong conformational pH dependence of chitinase I and the enzyme retained their characteristic TIM barrel topology at low pH.

Structuring Peptide Dendrimers through pH Modulation and Substrate Binding.

Dendrimers are a family of ramified synthetic molecules. pH effects and electrostatic interactions are known to be crucial players to explain the conformational and functional behaviors observed in these systems. Nonetheless, to date, no computational study involving these systems has explicitly addressed the protonation equilibrium taking place at different pH values for dendrimers containing multiple ionizable sites. Herein, we present the results of constant-pH molecular dynamics simulations performed at several pH values for four peptide dendrimers of different generations (from one to four) composed of the same type of amino acids: histidines, serines, and diaminopropionic acid. These dendrimers are known to catalyze the hydrolysis of pyrene sulfonate esters. Constant-pH MD simulations in the presence of substrate molecules at the optimum pH for catalysis are also reported. The results show that first and second generation dendrimers are almost structurally unresponsive to pH variations. For third and fourth generation dendrimers, pH plays a structuring role, with markedly different behaviors being observed when passing from acidic to neutral pH. Protonation-conformation coupling effects influence several intramolecular interactions, which, in turn, modulate the shape and structure at the different pH values. The atypical and highly pH-dependent protonation profiles of some histidine residues are also investigated. The interactions between dendrimers and substrates restrict the conformational space available to the dendrimers and enforce conformational homogeneity. This structuring effect is a consequence of the dendrimer-substrate interactions which occur through stabilizing hydrogen bonds and ion pairs between the substrate sulfonate groups and the dendrimer residues. Our results provide original fundamental data contributing to the development of novel pH-modulated dendritic systems and the improvement of the existing ones.

Unraveling HIV protease flaps dynamics by Constant pH Molecular Dynamics simulations.

The active site of HIV protease (HIV-PR) is covered by two flaps. These flaps are known to be essential for the catalytic activity of the HIV-PR, but their exact conformations at the different stages of the enzymatic pathway remain subject to debate. Understanding the correct functional dynamics of the flaps might aid the development of new HIV-PR inhibitors. It is known that, the HIV-PR catalytic efficiency is pH-dependent, likely due to the influence of processes such as charge transfer and protonation/deprotonation of ionizable residues. Several Molecular Dynamics (MD) simulations have reported information about the HIV-PR flaps. However, in MD simulations the protonation of a residue is fixed and thus it is not possible to study the correlation between conformation and protonation state. To address this shortcoming, this work attempts to capture, through Constant pH Molecular Dynamics (CpHMD), the conformations of the apo, substrate-bound and inhibitor-bound HIV-PR, which differ drastically in their flap arrangements. The results show that the HIV-PR flaps conformations are defined by the protonation of the catalytic residues Asp25/Asp25' and that these residues are sensitive to pH changes. This study suggests that the catalytic aspartates can modulate the opening of the active site and substrate binding.

Acid–Base Properties and Surface Charge Distribution of the Water-Soluble Au102(pMBA)44 Nanocluster

The pKa of the p-mercaptobenzoic acid (pMBA) ligands in the Au102(pMBA)44 nanocluster was measured by using acid–base and IR titration. The observed macroscopic pKa = 6.18 ± 0.05 is significantly more basic than that of free pMBA (pKa = 4.16), and the protonation behavior is anticooperative according to the Hill coefficient n = 0.64 ± 0.04. The cluster is truly water-soluble when more than 22 and insoluble when fewer than 7 ligands are in the deprotonated state. In order to obtain more insight into the anticooperative character, the cluster was modeled at pH ∼6.2 using constant pH molecular dynamics simulations. The pKa values of the individual pMBAs are in the range of 5.18–7.58, depending on the positions of the ligands, generating a nonuniform charge distribution on the cluster surface. The cluster interactions with inorganic and biological environment depend on the surface charge distribution, stressing the importance of understanding the protonation state of the cluster.

Charge-Neutral Constant pH Molecular Dynamics Simulations Using a Parsimonious Proton Buffer.

In constant pH molecular dynamics simulations, the protonation states of titratable sites can respond to changes of the pH and of their electrostatic environment. Consequently, the number of protons bound to the biomolecule, and therefore the overall charge of the system, fluctuates during the simulation. To avoid artifacts associated with a non-neutral simulation system, we introduce an approach to maintain neutrality of the simulation box in constant pH molecular dynamics simulations, while maintaining an accurate description of all protonation fluctuations. Specifically, we introduce a proton buffer that, like a buffer in experiment, can exchange protons with the biomolecule enabling its charge to fluctuate. To keep the total charge of the system constant, the uptake and release of protons by the buffer are coupled to the titration of the biomolecule with a constraint. We find that, because the fluctuation of the total charge (number of protons) of a typical biomolecule is much smaller than the number of titratable sites of the biomolecule, the number of buffer sites required to maintain overall charge neutrality without compromising the charge fluctuations of the biomolecule, is typically much smaller than the number of titratable sites, implying markedly enhanced simulation and sampling efficiency.

Tracking Changes in Protonation and Conformation during Photoactivation of a Phytochrome Protein

Phytochromes are photosensor proteins in plants and bacteria. The biological response is mediated by structural changes that follow photon absorption in the protein complex. The initial step is the photoisomerization of the biliverdin chromophore. How this leads to large-scale structural changes of the whole complex is, however, poorly understood. In this work, we use molecular dynamics (MD) simulations to investigate the structural changes after isomerization. In particular, we perform MD simulations at constant pH, using a recently developed method, to explore the effect of chromophore isomerization on the protonation (pKa) of nearby residues. In addition, we use a hybrid quantum mechanics/molecular mechanics approach to investigate the effect of isomerization, protonation and protein conformational changes on the absorption spectrum of the protein, for which experimental data are available. Here, we will first describe the constant pH MD simulations, and then compare the calculated spectra to experiment, and discuss the implications of our results for the photo-switching mechanism.

Constant-pH MD Simulations of DMPA/DMPC Lipid Bilayers.

Current constant-pH molecular dynamics (CpHMD) simulations provide a proper treatment of pH effects on the structure and dynamics of soluble biomolecules like peptides and proteins. However, addressing such effects on lipid membrane assemblies has remained problematic until now, despite the important role played by lipid ionization at physiological pH in a plethora of biological processes. Modeling (de)protonation events in these systems requires a proper consideration of the physicochemical features of the membrane environment, including a sound treatment of solution ions. Here, we apply our recent CpHMD-L method to the study of pH effects on a 25% DMPA/DMPC bilayer membrane model, closely reproducing the correct lipid phases of this system, namely, gel-fluid coexistence at pH 4 and a fluid phase at pH 7. A significant transition is observed for the membrane ionization and mechanical properties at physiological pH, providing a molecular basis for the well-established role of phosphatidic acid (PA) as a key player in the regulation of many cellular events. Also, as reported experimentally, we observed pH-induced PA-PA lipid aggregation at acidic pH. By including the titration of anionic phospholipids, the current methodology makes possible to simulate lipid bilayers with increased realism. To the best of our knowledge, this is the first simulation study dealing with a continuous phospholipid bilayer with pH titration of all constituent lipids.

The Effect of Membrane Environment on Surfactant Protein C Stability Studied by Constant-pH Molecular Dynamics.

Pulmonary surfactant protein C (SP-C) is a small peptide with two covalently linked fatty acyl chains that plays a crucial role in the formation and stabilization of the pulmonary surfactant reservoirs during the compression and expansion steps of the respiratory cycle. Although its function is known to be tightly related to its highly hydrophobic character and key interactions maintained with specific lipid components, much is left to understand about its molecular mechanism of action. Also, although it adopts a mainly helical structure while associated with the membrane, factors as pH variation and deacylation have been shown to affect its stability and function. In this work, the conformational behavior of both the acylated and deacylated SP-C isoforms was studied in a DPPC bilayer under different pH conditions using constant-pH molecular dynamics simulations. Our findings show that both protein isoforms are remarkably stable over the studied pH range, even though the acylated isoform exhibits a labile helix-turn-helix motif rarely observed in the other isoform. We estimate similar tilt angles for the two isoforms over the studied pH range, with a generally higher degree of internalization of the basic N-terminal residues in the deacylated case, and observe and discuss some protonation-conformation coupling effects. Both isoforms establish contacts with the surrounding lipid molecules (preferentially with the sn-2 ester bonds) and have a local effect on the conformational behavior of the surrounding lipid molecules, the latter being more pronounced for acylated SP-C.

Proton Dynamics on Goethite Nanoparticles and Coupling to Electron Transport.

The surface chemistry of metal oxide particles is governed by the charge that develops at the interface with aqueous solution. Mineral transformation, biogeochemical reactions, remediation, and sorption dynamics are profoundly affected in response. Here we report implementation of replica-exchange constant-pH molecular dynamics simulations that use classical molecular dynamics for exploring configurational space and Metropolis Monte Carlo walking through protonation space with a simulated annealing escape route from metastable configurations. By examining the archetypal metal oxide, goethite (α-FeOOH), we find that electrostatic potential gradients spontaneously arise between intersecting low-index crystal faces and across explicitly treated oxide nanoparticles at a magnitude exceeding the Johnson-Nyquist voltage fluctuation. Fluctuations in adsorbed proton density continuously repolarize the surface potential bias between edge-sharing crystal faces, at a rate slower than the reported electron-polaron hopping rate in goethite interiors. This suggests that these spontaneous surface potential fluctuations will control the net movement of charge carriers in the lattice.

Interpretation of pH-activity profiles for acid-base catalysis from molecular simulations.

The measurement of reaction rate as a function of pH provides essential information about mechanism. These rates are sensitive to the pK(a) values of amino acids directly involved in catalysis that are often shifted by the enzyme active site environment. Experimentally observed pH-rate profiles are usually interpreted using simple kinetic models that allow estimation of "apparent pK(a)" values of presumed general acid and base catalysts. One of the underlying assumptions in these models is that the protonation states are uncorrelated. In this work, we introduce the use of constant pH molecular dynamics simulations in explicit solvent (CpHMD) with replica exchange in the pH-dimension (pH-REMD) as a tool to aid in the interpretation of pH-activity data of enzymes and to test the validity of different kinetic models. We apply the methods to RNase A, a prototype acid-base catalyst, to predict the macroscopic and microscopic pK(a) values, as well as the shape of the pH-rate profile. Results for apo and cCMP-bound RNase A agree well with available experimental data and suggest that deprotonation of the general acid and protonation of the general base are not strongly coupled in transphosphorylation and hydrolysis steps. Stronger coupling, however, is predicted for the Lys41 and His119 protonation states in apo RNase A, leading to the requirement for a microscopic kinetic model. This type of analysis may be important for other catalytic systems where the active forms of the implicated general acid and base are oppositely charged and more highly correlated. These results suggest a new way for CpHMD/pH-REMD simulations to bridge the gap with experiments to provide a molecular-level interpretation of pH-activity data in studies of enzyme mechanisms.

Multiscale Modeling of a Conditionally Disordered pH-Sensing Chaperone

The pH-sensing chaperone HdeA promotes the survival of enteropathogenic bacteria during transit through the harshly acidic environment of the mammalian stomach. At low pH, HdeA transitions from an inactive, folded, dimer to chaperone-active, disordered, monomers to protect against the acid-induced aggregation of periplasmic proteins. Toward achieving a detailed mechanistic understanding of the pH response of HdeA, we develop a multiscale modeling approach to capture its pH-dependent thermodynamics. Our approach combines pKa (logarithmic acid dissociation constant) calculations from all-atom constant pH molecular dynamics simulations with coarse-grained modeling and yields new, atomic-level, insights into HdeA chaperone function that can be directly tested by experiment. “pH triggers” that significantly destabilize the dimer are each located near the N-terminus of a helix, suggesting that their neutralization at low pH destabilizes the helix macrodipole as a mechanism of monomer disordering. Moreover, we observe a non-monotonic change in the pH-dependent stability of HdeA, with maximal stability of the dimer near pH5. This affect is attributed to the protonation Glu37, which exhibits an anomalously high pKa value and is located within the hydrophobic dimer interface. Finally, the pH-dependent binding pathway of HdeA comprises a partially unfolded, dimeric intermediate that becomes increasingly stable relative to the native dimer at lower pH values and displays key structural features for chaperone–substrate interaction. We anticipate that the insights from our model will help inform ongoing NMR and biochemical investigations.