Simulating proteins at constant pH: An approach combining molecular dynamics and Monte Carlo simulation

For the structure and function of proteins, the pH of the solution is one of the determining parameters. Current molecular dynamics (MD) simulations account for the solution pH only in a limited way by keeping each titratable site in a chosen protonation state. We present an algorithm that generates trajectories at a Boltzmann distributed ensemble of protonation states by a combination of MD and Monte Carlo (MC) simulation. The algorithm is useful for pH-dependent structural studies and to investigate in detail the titration behavior of proteins. The method is tested on the acidic residues of the protein hen egg white lysozyme. It is shown that small structural changes may have a big effect on the pK(A) values of titratable residues.     

Effects of pH on kinetics of binding of mRNA-cap analogs by translation initiation factor eIF4E

Stopped-flow spectrofluorimetry and a theoretical method for predicting protonation equilibria in polyelectrolytes were combined in an analysis of the pH dependence of the kinetics of binding of analogues of the 5'-mRNA cap to the cap binding protein eIF4E. The computer simulations and available experimental data indicate that there are two titratable groups in the binding site of the protein and two titratable groups on the ligands directly involved in the binding, in addition to stacking interactions described by other groups. The observed pH dependencies of the rate constants obtained from the stopped-flow experiments are consistent with this finding. In particular, it is concluded that binding of both forms of the cap analogs regarding protonation at the N1 position of the guanine ring is efficient, and the shift to a predominantly protonated form of the ring takes place after formation of the complex.     

Multiple pH regime molecular dynamics simulation for pK calculations

Ionisation equilibria in proteins are influenced by conformational flexibility, which can in principle be accounted for by molecular dynamics simulation. One problem in this method is the bias arising from the fixed protonation state during the simulation. Its effect is mostly exhibited when the ionisation behaviour of the titratable groups is extrapolated to pH regions where the predetermined protonation state of the protein may not be statistically relevant, leading to conformational sampling that is not representative of the true state. In this work we consider a simple approach which can essentially reduce this problem. Three molecular dynamics structure sets are generated, each with a different protonation state of the protein molecule expected to be relevant at three pH regions, and pK calculations from the three sets are combined to predict pK over the entire pH range of interest. This multiple pH molecular dynamics approach was tested on the GCN4 leucine zipper, a protein for which a full data set of experimental data is available. The pK values were predicted with a mean deviation from the experimental data of 0.29 pH units, and with a precision of 0.13 pH units, evaluated on the basis of equivalent sites in the dimeric GCN4 leucine zipper.     

The influence of a transmembrane pH gradient on protonation probabilities of bacteriorhodopsin: the structural basis of the back-pressure effect

Bacteriorhodopsin pumps protons across a membrane using the energy of light. The proton pumping is inhibited when the transmembrane proton gradient that the protein generates becomes larger than four pH units. This phenomenon is known as the back-pressure effect. Here, we investigate the structural basis of this effect by predicting the influence of a transmembrane pH gradient on the titration behavior of bacteriorhodopsin. For this purpose we introduce a method that accounts for a pH gradient in protonation probability calculations. The method considers that in a transmembrane protein, which is exposed to two different aqueous phases, each titratable residue is accessible for protons from one side of the membrane depending on its hydrogen-bond pattern. This method is applied to several ground-state structures of bacteriorhodopsin, which residues already present complicated titration behaviors in the absence of a proton gradient. Our calculations show that a pH gradient across the membrane influences in a non-trivial manner the protonation probabilities of six titratable residues which are known to participate in the proton transfer: D85, D96, D115, E194, E204, and the Schiff base. The residues connected to one side of the membrane are influenced by the pH on the other side because of their long-range electrostatic interactions within the protein. In particular, D115 senses the pH at the cytoplasmic side of the membrane and transmits this information to D85 and the Schiff base. We propose that the strong electrostatic interactions found between D85, D115, and the Schiff base as well as the interplay of their respective protonation states under the influence of a transmembrane pH gradient are responsible for the back-pressure effect on bacteriorhodopsin.     

Assessing the acid-base and conformational properties of histidine residues in human prion protein (125-228) by means of pK(a) calculations and molecular dynamics simulations

A thorough study of the acid-base behavior of the four histidines and the other titratable residues of the structured domain of human prion protein (125-228) is presented. By using multi-tautomer electrostatic calculations, average titration curves have been built for all titratable residues, using the whole bundles of NMR structures determined at pH 4.5 and 7.0. According to our results, (1) only histidine residues are likely to be involved in the first steps of the pH-driven conformational transition of prion protein; (2) the pK(a)'s of His140 and His177 are approximately 7.0, whereas those of His155 and His187 are < 5.5. 10-ns long molecular dynamics simulations have been performed on five different models, corresponding to the most significant combinations of histidine protonation states. A critical comparison between the available NMR structures and our computational results (1) confirms that His155 and His187 are the residues whose protonation is involved in the conformational rearrangement of huPrP in mildly acidic condition, and (2) shows how their protonation leads to the destructuration of the C-terminal part of HB and to the loss of the last turn of HA that represent the crucial microscopic steps of the rearrangement.     

Computational investigation of proton transfer,pKa shifts and pH‐optimum of protein–DNA and protein–RNA complexes

Protein–nucleic acid interactions play a crucial role in many biological processes. This work investigates the changes of pKa values and protonation states of ionizable groups (including nucleic acid bases) that may occur at protein–nucleic acid binding. Taking advantage of the recently developed pKa calculation tool DelphiPka, we utilize the large protein–nucleic acid interaction database (NPIDB database) to model pKa shifts caused by binding. It has been found that the protein's interfacial basic residues experience favorable electrostatic interactions while the protein acidic residues undergo proton uptake to reduce the energy cost upon the binding. This is in contrast with observations made for protein–protein complexes. In terms of DNA/RNA, both base groups and phosphate groups of nucleotides are found to participate in binding. Some DNA/RNA bases undergo pKa shifts at complex formation, with the binding process tending to suppress charged states of nucleic acid bases. In addition, a weak correlation is found between the pH‐optimum of protein–DNA/RNA binding free energy and the pH‐optimum of protein folding free energy. Overall, the pH‐dependence of protein–nucleic acid binding is not predicted to be as significant as that of protein–protein association. Proteins 2017; 85:282–295. © 2016 Wiley Periodicals, Inc.     

Toward accurate prediction of pKa values for internal protein residues: the importance of conformational relaxation and desolvation energy

Proton uptake or release controls many important biological processes, such as energy transduction, virus replication, and catalysis. Accurate pK(a) prediction informs about proton pathways, thereby revealing detailed acid-base mechanisms. Physics-based methods in the framework of molecular dynamics simulations not only offer pK(a) predictions but also inform about the physical origins of pK(a) shifts and provide details of ionization-induced conformational relaxation and large-scale transitions. One such method is the recently developed continuous constant pH molecular dynamics (CPHMD) method, which has been shown to be an accurate and robust pK(a) prediction tool for naturally occurring titratable residues. To further examine the accuracy and limitations of CPHMD, we blindly predicted the pK(a) values for 87 titratable residues introduced in various hydrophobic regions of staphylococcal nuclease and variants. The predictions gave a root-mean-square deviation of 1.69 pK units from experiment, and there were only two pK(a)'s with errors greater than 3.5 pK units. Analysis of the conformational fluctuation of titrating side-chains in the context of the errors of calculated pK(a) values indicate that explicit treatment of conformational flexibility and the associated dielectric relaxation gives CPHMD a distinct advantage. Analysis of the sources of errors suggests that more accurate pK(a) predictions can be obtained for the most deeply buried residues by improving the accuracy in calculating desolvation energies. Furthermore, it is found that the generalized Born implicit-solvent model underlying the current CPHMD implementation slightly distorts the local conformational environment such that the inclusion of an explicit-solvent representation may offer improvement of accuracy.     

Ionization states of residues in OmpF and mutants: effects of dielectric constant and interactions between residues

To understand ion permeation, one must assign correct ionization states to titratable amino acid residues in protein channels. We report on the effects of physical and methodological assumptions in calculating the protonation states at neutral bulk pH of titratable residues lining the lumen of the native Escherichia coli OmpF channel, and five mutants. We systematically considered a wide range of assumed protein dielectric constants and all plausible combinations of protonation states for electrostatically interacting side chains, and three different levels of accounting for solute shielding: 1), full nonlinear Poisson-Boltzmann; 2), linearized Poisson-Boltzmann; and 3), neglect of solute shielding. For this system we found it acceptable to neglect solute shielding, a result we postulate to be generalizable to narrow lumens of other protein channels. For the large majority of residues, the protonation state at neutral bulk pH was found to be independent of the assumed dielectric constant of the protein, and unambiguously determined by the calculation; for native OmpF only Asp-127 has a protonation state that is sensitive to the assumed protein dielectric constant. Our results are significant for understanding two published experimental observations: the structure of the narrow part of the channel, and the ionic selectivity of OmpF mutants.     

Introducing Titratable Water to All-Atom Molecular Dynamics at Constant pH

Recent development of titratable coions has paved the way for realizing all-atom molecular dynamics at constant pH. To further improve physical realism, here we describe a technique in which proton titration of the solute is directly coupled to the interconversion between water and hydroxide or hydronium. We test the new method in replica-exchange continuous constant pH molecular dynamics simulations of three proteins, HP36, BBL, and HEWL. The calculated pK_a values based on 10-ns sampling per replica have the average absolute and root-mean-square errors of 0.7 and 0.9 pH units, respectively. Introducing titratable water in molecular dynamics offers a means to model proton exchange between solute and solvent, thus opening a door to gaining new insights into the intricate details of biological phenomena involving proton translocation.     

The pH-induced release of iron from transferrin investigated with a continuum electrostatic model.

A reduction in pH induces the release of iron from transferrin in a process that involves a conformational change in the protein from a closed to an open form. Experimental evidence suggests that there must be changes in the protonation states of certain, as yet not clearly identified, residues in the protein accompanying this conformational change. Such changes in protonation states of residues and the consequent changes in electrostatic interactions are assumed to play a large part in the mechanism of release of iron from transferrin. Using the x-ray crystal structures of human ferri- and apo-lactoferrin, we calculated the pKa values of the titratable residues in both the closed (iron-loaded) and open (iron-free) conformations with a continuum electrostatic model. With the knowledge of a residue's pKa value, its most probable protonation state at any specified pH may be determined. The preliminary results presented here are in good agreement with the experimental observation that the binding of ferric iron and the synergistic anion bicarbonate/carbonate results in the release of approximately three H+ ions. It is suggested that the release of these three H+ ions may be accounted for, in most part, by the deprotonation of the bicarbonate and residues Tyr-92, Lys-243, Lys-282, and Lys-285 together with the protonation of residues Asp-217 and Lys-277.     

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.     

Protonation linked equilibria and apparent affinity constants: the thermodynamic profile of the α-chymotrypsin–proflavin interaction

Protonation/deprotonation equilibria are frequently linked to binding processes involving proteins. The presence of these thermodynamically linked equilibria affects the observable thermodynamic parameters of the interaction (K _obs, ΔH _obs⁰). In order to try and elucidate the energetic factors that govern these binding processes, a complete thermodynamic characterisation of each intrinsic equilibrium linked to the complexation event is needed and should furthermore be correlated to structural information. We present here a detailed study, using NMR and ITC, of the interaction between α-chymotrypsin and one of its competitive inhibitors, proflavin. By performing proflavin titrations of the enzyme, at different pH values, we were able to highlight by NMR the effect of the complexation of the inhibitor on the ionisable residues of the catalytic triad of the enzyme. Using ITC we determined the intrinsic thermodynamic parameters of the different equilibria linked to the binding process. The possible driving forces of the interaction between α-chymotrypsin and proflavin are discussed in the light of the experimental data and on the basis of a model of the complex. This study emphasises the complementarities between ITC and NMR for the study of binding processes involving protonation/deprotonation equilibria. Electronic supplementary material The online version of this article (doi:) contains supplementary material, which is available to authorized users.     

Long-range nature of the interactions between titratable groups in Bacillus agaradhaerens family 11 xylanase: pH titration of B. agaradhaerens xylanase

Xylanase from Bacillus agaradhaerens belongs to a large group of glycosyl hydrolases which catalyze the degradation of xylan. The protonation behavior of titratable groups of the uniformly (15)N- and (13)C-labeled xylanase was investigated by multinuclear NMR spectroscopy. A total of 224 chemical shift titration curves corresponding to (1)H, (13)C, and (15)N resonances revealed pK(a) values for all aspartic and glutamic acid residues, as well as for the C-terminal carboxylate and histidine residues. Most of the titratable groups exhibit a complex titration behavior, which is most likely due to the mutual interactions with other neighboring groups or due to an unusual local microenvironment. Subsite -1 containing the catalytic dyad shows a long-range interaction over 9 A with Asp21 via two hydrogen bonds with Asn45 as the mediator. This result illuminates the pivotal role of the conserved position 45 among family 11 endoxylanases, determining an alkaline pH optimum by asparagine residues or an acidic pH optimum by an aspartate. The asymmetric interactions of neighboring tryptophan side chains with respect to the catalytic dyad can be comprehended as a result of hydrogen bonding and aromatic stacking. Most of the chemical shift-pH profiles of the backbone amides exhibit biphasic behavior with two distinct inflection points, which correspond to the pK(a) values of the nearby acidic side chains. However, the alternation of both positive and negative slopes of individual amide titration curves is interpreted as a consequence of a simultaneous reorganization of side chain conformational space at pH approximately 6 and/or an overall change in the hydrogen network in the substrate binding cleft.     

Acid-base equilibria in rhodopsin: dependence of the protonation state of glu134 on its environment

Glutamic acid E134 in rhodopsin is part of a highly conserved triad, D(E)RY, located near the cytoplasmic lipid/water interface in transmembrane helix 3 of G protein-coupled receptors (GPCRs). A large body of experimental evidence suggests that the protonation of E134 plays a role in the mechanism of activation of rhodopsin and other GPCRs as well. For E134 to change its protonation state, its pK(a) value must shift from values below physiological pH to higher values. Because of the proximity of the triad to the lipid/water interface, it was hypothesized that a change in solvent around E134 from water to lipid could induce such a shift in pK(a). To test this hypothesis, the pK(a) values of the titratable amino acid residues in rhodopsin have been calculated and the change in solvent around E134 was modeled by shifting the position of the lipid/water interface. The approach used to carry out the pK(a) calculations takes into account the partial immersion of transmembrane proteins in lipid. Qualitative experimental evidence is available for several residues regarding their likely protonation state in rhodopsin at or near physiological pH. Comparison of the calculated pK(a) values with these experimental findings shows good agreement between the two. Notably, glutamic acids E122 and E181 were found to be protonated. The pK(a) values were then calculated for a range of lipid/water interface positions. Although the surrounding solvent of several titratable residues changed from water to lipid in this range, leading to pK(a) shifts in most cases, only for E134 would the shift lead to a change in protonation state at physiological pH. Thus, our results show that the protonation state of E134 is particularly sensitive to its environment. This sensitivity together with the location of E134 near the actual position of the lipid/water interface could be a strategic element in the mechanism of activation of rhodopsin.     

Modeling and simulation of biological systems with stochasticity

Mathematical modeling is a powerful approach for understanding the complexity of biological systems. Recently, several successful attempts have been made for simulating complex biological processes like metabolic pathways, gene regulatory networks and cell signaling pathways. The pathway models have not only generated experimentally verifiable hypothesis but have also provided valuable insights into the behavior of complex biological systems. Many recent studies have confirmed the phenotypic variability of organisms to an inherent stochasticity that operates at a basal level of gene expression. Due to this reason, development of novel mathematical representations and simulations algorithms are critical for successful modeling efforts in biological systems. The key is to find a biologically relevant representation for each representation. Although mathematically rigorous and physically consistent, stochastic algorithms are computationally expensive, they have been successfully used to model probabilistic events in the cell. This paper offers an overview of various mathematical and computational approaches for modeling stochastic phenomena in cellular systems.     

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.     

Exploring the early stages of the pH‐induced conformational change of influenza hemagglutinin

Hemagglutinin (HA) mediates the membrane fusion process of influenza virus through its pH‐induced conformational change. However, it remains challenging to study its structure reorganization pathways in atomic details. Here, we first applied continuous constant pH molecular dynamics approach to predict the pK_a values of titratable residues in H2 subtype HA. The calculated net‐charges in HA1 globular heads increase from 0e (pH 7.5) to +14e (pH 4.5), indicating that the charge repulsion drives the detrimerization of HA globular domains. In HA2 stem regions, critical pH sensors, such as Glu103₂, His18₁, and Glu89₁, are identified to facilitate the essential structural reorganizations in the fusing pathways, including fusion peptide release and interhelical loop transition. To probe the contribution of identified pH sensors and unveil the early steps of pH‐induced conformational change, we carried out conventional molecular dynamics simulations in explicit water with determined protonation state for each titratable residue in different environmental pH conditions. Particularly, energy barriers involving previously uncharacterized hydrogen bonds and hydrophobic interactions are identified in the fusion peptide release pathway. Nevertheless, comprehensive comparisons across HA family members indicate that different HA subtypes might employ diverse pH sensor groups along with different fusion pathways. Finally, we explored the fusion inhibition mechanism of antibody CR6261 and small molecular inhibitor TBHQ, and discovered a novel druggable pocket in H2 and H5 subtypes. Our results provide the underlying mechanism for the pH‐driven conformational changes and also novel insight for anti‐flu drug development. Proteins 2014; 82:2412–2428. © 2014 Wiley Periodicals, Inc.     

Oxidation of the non-heme iron complex in photosystem II

In photosystem II (PSII), the redox properties of the non-heme iron complex (Fe complex) are sensitive to the redox state of quinones (Q(A/)(B)), which may relate to the electron/proton transfer. We calculated the redox potentials for one-electron oxidation of the Fe complex in PSII [E(m)(Fe)] based on the reference value E(m)(Fe) = +400 mV at pH 7 in the Q(A)(0)Q(B)(0) state, considering the protein environment in atomic detail and the associated changes in protonation pattern. Our model yields the pH dependence of E(m)(Fe) with -60 mV/pH as observed in experimental redox titration. We observed significant deprotonation at D1-Glu244 in the hydrophilic loop region upon Fe complex oxidation. The calculated pK(a) value for D1-Glu244 depends on the Fe complex redox state, yielding a pK(a) of 7.5 and 5.5 for Fe(2+) and Fe(3+), respectively. To account for the pH dependence of E(m)(Fe), a model involving not only D1-Glu244 but also the other titratable residues (five Glu in the D-de loops and six basic residues near the Fe complex) seems to be needed, implying the existence of a network of residues serving as an internal proton reservoir. Reduction of Q(A/B) yields +302 mV and +268 mV for E(m)(Fe) in the Q(A)(-)Q(B)(0) and Q(A)(0)Q(B)(-) states, respectively. Upon formation of the Q(A)(0)Q(B)(-) state, D1-His252 becomes protonated. Forming Fe(3+)Q(B)H(2) by a proton-coupled electron transfer process from the initial state Fe(2+)Q(B)(-) results in deprotonation of D1-His252. The two EPR signals observed at g = 1.82 and g = 1.9 in the Fe(2+)Q(A)(-) state of PSII may be attributed to D1-His252 with variable and fixed protonation, respectively.