Predicting pKa values with continuous constant pH molecular dynamics

Anions and Cations Affect Amino Acid Dissociation Equilibria via Distinct Mechanisms

Salts reduce the pKa of weak acids by a mechanism sensitive to ion identity and concentration via charge screening of the deprotonated state. In this study, we utilize constant pH molecular dynamics simulations to understand the molecular mechanism behind the salt-dependent dissociation of aspartic acid (Asp). We calculate the pKa of Asp in the presence of a monovalent salt and investigate Hofmeister ion effects by systematically varying the ionic radii. We observe that increasing the anion size leads to a monotonic decrease in Asp pKa. Conversely, the cation size affects the pKa nonmonotonically, interpretable in the context of the law of matching water affinity. The net effect of salt on Asp acidity is governed by an interplay of solvation and competing ion interactions. The proposed mechanism is rather general and can be applicable to several problems in Hofmeister ion chemistry, such as pH effects on protein stability and soft matter interfaces.

Understanding the Key Roles of pH Buffer in Accelerating Lignin Degradation by Lignin Peroxidase

pH buffer plays versatile roles in both biology and chemistry. In this study, we unravel the critical role of pH buffer in accelerating degradation of the lignin substrate in lignin peroxidase (LiP) using QM/MM MD simulations and the nonadiabatic electron transfer (ET) and proton-coupled electron transfer (PCET) theories. As a key enzyme involved in lignin degradation, LiP accomplishes the oxidation of lignin via two consecutive ET reactions and the subsequent C-C cleavage of the lignin cation radical. The first one involves ET from Trp171 to the active species of Compound I, while the second one involves ET from the lignin substrate to the Trp171 radical. Differing from the common view that pH = 3 may enhance the oxidizing power of Cpd I via protonation of the protein environment, our study shows that the intrinsic electric fields have minor effects on the first ET step. Instead, our study shows that the pH buffer of tartaric acid plays key roles during the second ET step. Our study shows that the pH buffer of tartaric acid can form a strong H-bond with Glu250, which can prevent the proton transfer from the Trp171-H^•+ cation radical to Glu250, thereby stabilizing the Trp171-H^•+ cation radical for the lignin oxidation. In addition, the pH buffer of tartaric acid can enhance the oxidizing power of the Trp171-H^•+ cation radical via both the protonation of the proximal Asp264 and the second-sphere H-bond with Glu250. Such synergistic effects of pH buffer facilitate the thermodynamics of the second ET step and reduce the overall barrier of lignin degradation by ∼4.3 kcal/mol, which corresponds to a rate acceleration of 10³-fold that agrees with experiments. These findings not only expand our understanding on pH-dependent redox reactions in both biology and chemistry but also provide valuable insights into tryptophan-mediated biological ET reactions.

Biologically relevant small variations of intra-cellular pH can have significant effect on stability of protein-DNA complexes, including the nucleosome

Stability of a protein-ligand complex may be sensitive to pH of its environment. Here we explore, computationally, stability of a set of protein-nucleic acid complexes using fundamental thermodynamic linkage relationship. The nucleosome, as well as an essentially random selection of 20 protein complexes with DNA or RNA, are included in the analysis. An increase in intra-cellular/intra-nuclear pH destabilizes most complexes, including the nucleosome. We propose to quantify the effect by ΔΔG^0.3-the change in the binding free energy due to pH increase of 0.3 units, corresponding to doubling of the H ⁺ activity; variations of pH of this amplitude can occur in living cells, including in the course of the cell cycle, and in cancer cells relative to normal ones. We suggest, based on relevant experimental findings, a threshold of biological significance of 1 2 k B T ( ∼ 0.3 k c a l / m o l ) for changes of stability of chromatin-related protein-DNA complexes: a change in the binding affinity above the threshold may have biological consequences. We find that for 70% of the examined complexes, Δ Δ G 0.3 > 1 2 k B T (for 10%, ΔΔG^0.3 is between 3 and 4 k _B T). Thus, small but relevant variations of intra-nuclear pH of 0.3 may have biological consequences for many protein-nucleic acid complexes. The binding affinity between the histone octamer and its DNA, which directly affects the DNA accessibility in the nucleosome, is predicted to be highly sensitive to intra-nuclear pH. A variation of 0.3 units results in ΔΔG^0.3 ∼ 10k _B T ( ∼ 6 k c a l / m o l ) ; for spontaneous unwrapping of 20 bp long entry/exit fragments of the nucleosomal DNA, ΔΔG^0.3 = 2.2k _B T; partial disassembly of the nucleosome into the tetrasome is characterized by ΔΔG^0.3 = 5.2k _B T. The predicted pH -induced modulations of the nucleosome stability are significant enough to suggest that they may have consequences relevant to the biological function of the nucleosome. Accessibility of the nucleosomal DNA is predicted to positively correlate with pH variations during the cell cycle; an increase in intra-cellular pH seen in cancer cells is predicted to lead to a more accessible nucleosomal DNA; a drop in pH associated with apoptosis is predicted to make nucleosomal DNA less accessible. We speculate that processes that depend on accessibility to the DNA in the nucleosomes, such as transcription or DNA replication, might become upregulated due to relatively small, but nevertheless realistic increases of intra-nuclear pH.

Galvani Offset Potential and Constant-pH Simulations of Membrane Proteins

A central problem in computational biophysics is the treatment of titratable residues in molecular dynamics simulations of large biological macromolecular systems. Conventional simulation methods ascribe a fixed ionization state to titratable residues in accordance with their pK_a and the pH of the system, assuming that an effective average model will be able to capture the predominant behavior of the system. While this assumption may be justifiable in many cases, it is certainly limited, and it is important to design alternative methodologies allowing a more realistic treatment. Constant-pH simulation methods provide powerful approaches to handle titratable residues more realistically by allowing the ionization state to vary statistically during the simulation. Extending the molecular mechanical (MM) potential energy function to a family of potential functions accounting for different ionization states, constant-pH simulations are designed to sample all accessible configurations and ionization states, properly weighted according to their Boltzmann factor. Because protonation and deprotonation events correspond to a change in the total charge, difficulties arise when the long-range Coulomb interaction is treated on the basis of an idealized infinite simulation model and periodic boundary conditions with particle-mesh Ewald lattice sums. Charging free-energy calculations performed under these conditions in aqueous solution depend on the Galvani potential of the bulk water phase. This has important implications for the equilibrium and nonequilibrium constant-pH simulation methods grounded in the relative free-energy difference corresponding to the protonated and unprotonated residues. Here, the effect of the Galvani potential is clarified, and a simple practical solution is introduced to address this issue in constant-pH simulations of the acid-sensing ion channel (ASIC).

Characterizing Protein Protonation Microstates Using Monte Carlo Sampling

Proteins are polyelectrolytes with acidic and basic amino acids Asp, Glu, Arg, Lys, and His, making up ≈25% of the residues. The protonation state of residues, cofactors, and ligands defines a "protonation microstate". In an ensemble of proteins some residues will be ionized and others neutral, leading to a mixture of protonation microstates rather than in a single one as is often assumed. The microstate distribution changes with pH. The protein environment also modifies residue proton affinity so microstate distributions change in different reaction intermediates or as ligands are bound. Particular protonation microstates may be required for function, while others exist simply because there are many states with similar energy. Here, the protonation microstates generated in Monte Carlo sampling in MCCE are characterized in HEW lysozyme as a function of pH and bacterial photosynthetic reaction centers (RCs) in different reaction intermediates. The lowest energy and highest probability microstates are compared. The ΔG, ΔH, and ΔS between the four protonation states of Glu35 and Asp52 in lysozyme are shown to be calculated with reasonable precision. At pH 7 the lysozyme charge ranges from 6 to 10, with 24 accepted protonation microstates, while RCs have ≈50,000. A weighted Pearson correlation analysis shows coupling between residue protonation states in RCs and how they change when the quinone in the Q_B site is reduced. Protonation microstates can be used to define input MD parameters and provide insight into the motion of protons coupled to reactions.

Ribose and Non-Ribose A2A Adenosine Receptor Agonists: Do They Share the Same Receptor Recognition Mechanism?

Adenosine receptors have been a promising class of targets for the development of new therapies for several diseases. In recent years, a renewed interest in this field has risen, thanks to the implementation of a novel class of agonists that lack the ribose moiety, once considered essential for the agonistic profile. Recently, an X-ray crystal structure of the A_2A adenosine receptor has been solved, providing insights about the receptor activation from this novel class of agonists. Starting from this structural information, we have performed supervised molecular dynamics (SuMD) simulations to investigate the binding pathway of a non-nucleoside adenosine receptor agonist as well as one of three classic agonists. Furthermore, we analyzed the possible role of water molecules in receptor activation.

Hybrid MC/MD for protein design

Computational protein design relies on simulations of a protein structure, where selected amino acids can mutate randomly, and mutations are selected to enhance a target property, such as stability. Often, the protein backbone is held fixed and its degrees of freedom are modeled implicitly to reduce the complexity of the conformational space. We present a hybrid method where short molecular dynamics (MD) segments are used to explore conformations and alternate with Monte Carlo (MC) moves that apply mutations to side chains. The backbone is fully flexible during MD. As a test, we computed side chain acid/base constants or pK_a's in five proteins. This problem can be considered a special case of protein design, with protonation/deprotonation playing the role of mutations. The solvent was modeled as a dielectric continuum. Due to cost, in each protein we allowed just one side chain position to change its protonation state and the other position to change its type or mutate. The pK_a's were computed with a standard method that scans a range of pH values and with a new method that uses adaptive landscape flattening (ALF) to sample all protonation states in a single simulation. The hybrid method gave notably better accuracy than standard, fixed-backbone MC. ALF decreased the computational cost a factor of 13.

Protonation Dynamics in the K-Channel of Cytochrome c Oxidase Estimated from Molecular Dynamics Simulations

Proton transfer reactions are one of the most fundamental processes in biochemistry. We present a simplistic approach for estimating proton transfer probabilities in a membrane protein, cytochrome c oxidase. We combine short molecular dynamics simulations at discrete protonation states with a Monte Carlo approach to exchange between those states. Requesting for a proton transfer the existence of a hydrogen-bonded connection between the two source and target residues of the exchange, restricts the acceptance of transfers to only those in which a proton-relay is possible. Together with an analysis of the hydrogen-bonded connectivity in one of the proton-conducting channels of cytochrome c oxidase, this approach gives insight into the protonation dynamics of the hydrogen-bonded networks. The connectivity and directionality of the networks are coupled to the conformation of an important protein residue in the channel, K362, rendering proton transfer in the entire channel feasible in only one of the two major conformations. Proton transport in the channel can thus be regulated by K362 not only through its possible role as a proton carrier itself, but also by allowing or preventing proton transport via water residues.

Continuous Constant pH Molecular Dynamics Simulations of Transmembrane Proteins

Many membrane channels, transporters, and receptors utilize a pH gradient or proton coupling to drive functionally relevant conformational transitions. Conventional molecular dynamics simulations employ fixed protonation states, thus neglecting the coupling between protonation and conformational equilibria. Here we describe the membrane-enabled hybrid-solvent continuous constant pH molecular dynamics method for capturing atomic details of proton-coupled conformational dynamics of transmembrane proteins. Example protocols from our recent application studies of proton channels and ion/substrate transporters are discussed.

Predicting Reactive Cysteines with Implicit-Solvent-Based Continuous Constant pH Molecular Dynamics in Amber

Cysteines existing in the deprotonated thiolate form or having a tendency to become deprotonated are important players in enzymatic and cellular redox functions and frequently exploited in covalent drug design; however, most computational studies assume cysteines as protonated. Thus, developing an efficient tool that can make accurate and reliable predictions of cysteine protonation states is timely needed. We recently implemented a generalized Born (GB) based continuous constant pH molecular dynamics (CpHMD) method in Amber for protein pKa calculations on CPUs and GPUs. Here we benchmark the performance of GB-CpHMD for predictions of cysteine pKa's and reactivities using a data set of 24 proteins with both down- and upshifted cysteine pKa's. We found that 10 ns single-pH or 4 ns replica-exchange CpHMD titrations gave root-mean-square errors of 1.2-1.3 and correlation coefficients of 0.8-0.9 with respect to experiment. The accuracy of predicting thiolates or reactive cysteines at physiological pH with single-pH titrations is 86 or 81% with a precision of 100 or 90%, respectively. This performance well surpasses the traditional structure-based methods, particularly a widely used empirical pKa tool that gives an accuracy less than 50%. We discuss simulation convergence, dependence on starting structures, common determinants of the pKa downshifts and upshifts, and the origin of the discrepancies from the structure-based calculations. Our work suggests that CpHMD titrations can be performed on a desktop computer equipped with a single GPU card to predict cysteine protonation states for a variety of applications, from understanding biological functions to covalent drug design.

GPU-Accelerated Implementation of Continuous Constant pH Molecular Dynamics in Amber: pKa Predictions with Single-pH Simulations

We present a GPU implementation of the continuous constant pH molecular dynamics (CpHMD) based on the most recent generalized Born implicit-solvent model in the pmemd engine of the Amber molecular dynamics package. To test the accuracy of the tool for rapid pK_a predictions, a series of 2 ns single-pH simulations were performed for over 120 titratable residues in 10 benchmark proteins that were previously used to test the various continuous CpHMD methods. The calculated pK_a's showed a root-mean-square deviation of 0.80 and correlation coefficient of 0.83 with respect to experiment. Also, 90% of the pK_a's were converged with estimated errors below 0.1 pH units. Surprisingly, this level of accuracy is similar to our previous replica-exchange simulations with 2 ns per replica and an exchange attempt frequency of 2 ps^-1 (Huang, Harris, and Shen J. Chem. Inf. Model. 2018 , 58 , 1372 - 1383 ). Interestingly, for the linked titration sites in two enzymes, although residue-specific protonation state sampling in the single-pH simulations was not converged within 2 ns, the protonation fraction of the linked residues appeared to be largely converged, and the experimental macroscopic pK_a values were reproduced to within 1 pH unit. Comparison with replica-exchange simulations with different exchange attempt frequencies showed that the splitting between the two macroscopic pK_a's is underestimated with frequent exchange attempts such as 2 ps^-1, while single-pH simulations overestimate the splitting. The same trend is seen for the single-pH vs replica-exchange simulations of a hydrogen-bonded aspartyl dyad in a much larger protein. A 2 ns single-pH simulation of a 400-residue protein takes about 1 h on a single NVIDIA GeForce RTX 2080 graphics card, which is over 1000 times faster than a CpHMD run on a single CPU core of a high-performance computing cluster node. Thus, we envision that GPU-accelerated continuous CpHMD may be used in routine pK_a predictions for a variety of applications, from assisting MD simulations with protonation state assignment to offering pH-dependent corrections of binding free energies and identifying reactive hot spots for covalent drug design.

Computational chemistry methods to investigate the effects caused by DNA variants linked with disease

Computational chemistry offers variety of tools to study properties of biological macromolecules. These tools vary in terms of levels of details from quantum mechanical treatment to numerous macroscopic approaches. Here, we provide a review of computational chemistry algorithms and tools for modeling the effects of genetic variations and their association with diseases. Particular emphasis is given on modeling the effects of missense mutations on stability, conformational dynamics, binding, hydrogen bond network, salt bridges, and pH-dependent properties of the corresponding macromolecules. It is outlined that the disease may be caused by alteration of one or several of above-mentioned biophysical characteristics, and a successful prediction of pathogenicity requires detailed analysis of how the alterations affect the function of involved macromolecules. The review provides a short list of most commonly used algorithms to predict the molecular effects of mutations as well.

Modeling and Simulation of Oligonucleotide Hybrids: Outlining a Strategy

Molecular dynamics simulations with a state-of-the-art force field provide an atomistic detailed description of the structural and thermodynamic features of biomolecules. Effects of chemical modifications and of the environment such as sequence, solvent, and ionic strength can explicitly be taken into account. Molecular simulation techniques can also provide insight in change in binding affinity, in protonation (pK_a shift) and tautomeric propensity due to changes in the environment or in the molecular system. The quality and reliability of a simulation depend on the quality of the force field and on the reproducibility of the data, and validation depends on the availability of suitable experimental data. Here, we describe the workflow to investigate oligonucleotide hybrids using molecular simulation including hardware and software information.

Reservoir pH replica exchange

We present the reservoir pH replica exchange (R-pH-REM) method for constant pH simulations. The R-pH-REM method consists of a two-step procedure; the first step involves generation of one or more reservoirs of conformations. Each reservoir is obtained from a standard or enhanced molecular dynamics simulation with a constrained (fixed) protonation state. In the second step, fixed charge constraints are relaxed, as the structures from one or more reservoirs are periodically injected into a constant pH or a pH-replica exchange (pH-REM) simulation. The benefit of this two-step process is that the computationally intensive part of conformational search can be decoupled from constant pH simulations, and various techniques for enhanced conformational sampling can be applied without the need to integrate such techniques into the pH-REM framework. Simulations on blocked Lys, KK, and KAAE peptides were used to demonstrate an agreement between pH-REM and R-pH-REM simulations. While the reservoir simulations are not needed for these small test systems, the real need arises in cases when ionizable molecules can sample two or more conformations separated by a large energy barrier, such that adequate sampling is not achieved on a time scale of standard constant pH simulations. Such problems might be encountered in protein systems that exploit conformational transitions for function. A hypothetical case is studied, a small molecule with a large torsional barrier; while results of pH-REM simulations depend on the starting structure, R-pH-REM calculations on this model system are in excellent agreement with a theoretical model.

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 K_a's of the buried and hydrogen bonded catalytic dyads in five enzymes. We asked the question what determines the p K_a order, i.e., what makes a residue proton donor vs a nucleophile. The continuous constant pH molecular dynamics simulations captured the experimental p K_a 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.

Computational Study of Uracil Tautomeric Forms in the Ribosome: The Case of Uracil and 5-Oxyacetic Acid Uracil in the First Anticodon Position of tRNA

Tautomerism is important in many biomolecular interactions, not least in RNA biology. Crystallographic studies show the possible presence of minor tautomer forms of transfer-RNA (tRNA) anticodon bases in the ribosome. The hydrogen positions are not resolved in the X-ray studies, and we have used ab initio calculations and molecular dynamics simulations to understand if and how the minor enol form of uracil (U), or the modified uracil 5-oxyacetic acid (cmo⁵U), can be accommodated in the tRNA-messenger-RNA interactions in the ribosome decoding center. Ab initio calculations on isolated bases show that the modification affects the keto-enol equilibrium of the uracil base only slightly; the keto form is dominant (>99.99%) in both U and cmo⁵U. Other factors such as interactions with the surrounding nucleotides or ions would be required to shift the equilibrium toward the enol tautomer. Classical molecular simulations show a better agreement with the X-ray structures for the enol form, but free energy calculations indicate that the most stable form is the keto. In the ribosome, the enol tautomers of U and cmo⁵U pair with a guanine forming two hydrogen bonds, which do not involve the enol group. The oxyacetic acid modification has a minor effect on the keto-enol equilibrium.

Constant pH Molecular Dynamics Reveals How Proton Release Drives the Conformational Transition of a Transmembrane Efflux Pump

AcrB is the inner-membrane transporter of an E. coli AcrAB-TolC tripartite efflux complex, which plays a major role in the intrinsic resistance to clinically important antibiotics. AcrB pumps a wide range of toxic substrates by utilizing the proton gradient between periplasm and cytoplasm. Crystal structures of AcrB revealed three distinct conformational states of the transport cycle, substrate access, binding, and extrusion or loose (L), tight (T), and open (O) states. However, the specific residue(s) responsible for proton binding/release and the mechanism of proton-coupled conformational cycling remain controversial. Here we use the newly developed membrane hybrid-solvent continuous constant pH molecular dynamics technique to explore the protonation states and conformational dynamics of the transmembrane domain of AcrB. Simulations show that both Asp407 and Asp408 are deprotonated in the L/T states, while only Asp408 is protonated in the O state. Remarkably, release of a proton from Asp408 in the O state results in large conformational changes, such as the lateral and vertical movement of transmembrane helices as well as the salt-bridge formation between Asp408 and Lys940 and other side chain rearrangements among essential residues. Consistent with the crystallographic differences between the O and L protomers, simulations offer dynamic details of how proton release drives the O-to-L transition in AcrB and address the controversy regarding the proton/drug stoichiometry. This work offers a significant step toward characterizing the complete cycle of proton-coupled drug transport in AcrB and further validates the membrane hybrid-solvent CpHMD technique for studies of proton-coupled transmembrane proteins which are currently poorly understood.

Development of constant-pH simulation methods in implicit solvent and applications in biomolecular systems

pH is a critical parameter for biological and technological systems directly related with electrical charges. It can give rise to peculiar electrostatic phenomena, which also makes them more challenging. Due to the quantum nature of the process, involving the forming and breaking of chemical bonds, quantum methods should ideally by employed. Nevertheless, due to the very large number of ionizable sites, different macromolecular conformations, salt conditions, and all other charged species, the CPU time cost simply becomes prohibitive for computer simulations, making this a quite complex problem. Simplified methods based on Monte Carlo sampling have been devised and will be reviewed here, highlighting the updated state-of-the-art of this field, advantages, and limitations of different theoretical protocols for biomolecular systems (proteins and nucleic acids). Following a historical perspective, the discussion will be associated with the applications to protein interactions with other proteins, polyelectrolytes, and nanoparticles.

Comparing pairwise-additive and many-body generalized Born models for acid/base calculations and protein design

Generalized Born (GB) solvent models are common in acid/base calculations and protein design. With GB, the interaction between a pair of solute atoms depends on the shape of the protein/solvent boundary and, therefore, the positions of all solute atoms, so that GB is a many-body potential. For compute-intensive applications, the model is often simplified further, by introducing a mean, native-like protein/solvent boundary, which removes the many-body property. We investigate a method for both acid/base calculations and protein design that uses Monte Carlo simulations in which side chains can explore rotamers, bind/release protons, or mutate. The fluctuating protein/solvent dielectric boundary is treated in a way that is numerically exact (within the GB framework), in contrast to a mean boundary. Its originality is that it captures the many-body character while retaining the residue-pairwise complexity given by a fixed boundary. The method is implemented in the Proteus protein design software. It yields a slight but systematic improvement for acid/base constants in nine proteins and a significant improvement for the computational design of three PDZ domains. It eliminates a source of model uncertainty, which will facilitate the analysis of other model limitations. © 2017 Wiley Periodicals, Inc.

Role of Pseudoisocytidine Tautomerization in Triplex-Forming Oligonucleotides: In Silico and in Vitro Studies

Pseudoisocytidine (^ΨC) is a synthetic cytidine analogue that can target DNA duplex to form parallel triplex at neutral pH. Pseudoisocytidine has mainly two tautomers, of which only one is favorable for triplex formation. In this study, we investigated the effect of sequence on ^ΨC tautomerization using λ-dynamics simulation, which takes into account transitions between states. We also performed in vitro binding experiments with sequences containing ^ΨC and furthermore characterized the structure of the formed triplex using molecular dynamics simulation. We found that the neighboring methylated or protonated cytidine promotes the formation of the favorable tautomer, whereas the neighboring thymine or locked nucleic acid has a poor effect, and consecutive ^ΨC has a negative influence. The deleterious effect of consecutive ^ΨC in a triplex formation was confirmed using in vitro binding experiments. Our findings contribute to improving the design of ^ΨC-containing triplex-forming oligonucleotides directed to target G-rich DNA sequences.

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.

All-Atom Continuous Constant pH Molecular Dynamics With Particle Mesh Ewald and Titratable Water

Development of a pH stat to properly control solution pH in biomolecular simulations has been a long-standing goal in the community. Toward this goal recent years have witnessed the emergence of the so-called constant pH molecular dynamics methods. However, the accuracy and generality of these methods have been hampered by the use of implicit-solvent models or truncation-based electrostatic schemes. Here we report the implementation of the particle mesh Ewald (PME) scheme into the all-atom continuous constant pH molecular dynamics (CpHMD) method, enabling CpHMD to be performed with a standard MD engine at a fractional added computational cost. We demonstrate the performance using pH replica-exchange CpHMD simulations with titratable water for a stringent test set of proteins, HP36, BBL, HEWL, and SNase. With the sampling time of 10 ns per replica, most pK_a's are converged, yielding the average absolute and root-mean-square deviations of 0.61 and 0.77, respectively, from experiment. Linear regression of the calculated vs experimental pK_a shifts gives a correlation coefficient of 0.79, a slope of 1, and an intercept near 0. Analysis reveals inadequate sampling of structure relaxation accompanying a protonation-state switch as a major source of the remaining errors, which are reduced as simulation prolongs. These data suggest PME-based CpHMD can be used as a general tool for pH-controlled simulations of macromolecular systems in various environments, enabling atomic insights into pH-dependent phenomena involving not only soluble proteins but also transmembrane proteins, nucleic acids, surfactants, and polysaccharides.

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.

Finite-size effect on the charging free energy of protein in explicit solvent

The finite-size effect in periodic system is examined for the charging free energy of protein in explicit solvent over a variety of charged states. The key to the finite-size correction is the self-energy, which is defined as the interaction energy of the solute with its own periodic images and the neutralizing background. By employing the thermodynamic-integration method with systematically varied sizes of the unit cell of molecular dynamics (MD) simulations, we show for ubiquitin that the self-energy corrects the finite-size effect on the charging free energy within 1 kcal/mol at total charges of -5e, -1e, neutral, and +1e and within 5 kcal/mol even for a highly charged state with +8e. We then sought the additional correction from the solvation effect using the numerical solution to the Poisson equation of the protein with implicit solvent. This correction reduces the cell-size dependence of the charging free energy at +8e to 3 kcal/mol and is well expressed as the self-energy divided by the dielectric constant of solvent water.

The AIBLHiCoS Method: Predicting Aqueous pKa Values from Gas-Phase Equilibrium Bond Lengths

The proposed AIBLHiCoS method predicts a given compound's pKa in aqueous solution from a single ab initio bond length only, after geometry optimization in the gas phase. Here we provide simple and predictive equations for naphthols and chemically similar biomolecules. Each linear equation corresponds to a High-Correlation Subset (HiCoS) that expresses the novel type of linear free energy relationship discovered here. The naphthol family exhibits a clear and strong relationship with the phenol family, with the "active" C-O bond always producing the highest correlations. The proposed method can isolate erroneous experiments and operate in non-aqueous solution and at different temperatures. Moreover, the existence of "active fragments" is demonstrated in a variety of sizable biomolecules for which the pKa is successfully predicted.

Dehydration of multilamellar fatty acid membranes: towards a computational model of the stratum corneum

The level of hydration controls the cohesion between apposed lamellae of saturated free fatty acids found in the lipid matrix of stratum corneum, the outermost layer of mammalian skin. This multilamellar lipid matrix is highly impermeable to water and ions, so that the local hydration shell of its fatty acids may not always be in equilibrium with the acidity and relative humidity, which significantly change over a course of days during skin growth. The homeostasis of the stratum corneum at each moment of its growth likely requires a balance between two factors, which affect in opposite ways the diffusion of hydrophilic species through the stratum corneum: (i) an increase in water order as the lipid lamellae come in closer contact, and (ii) a decrease in water order as the fraction of charged fatty acids is lowered by pH. Herein molecular dynamics simulations are employed to estimate the impact of both effects on water molecules confined between lamellae of fatty acids. Under conditions where membrane undulations are energetically favorable, the charged fatty acids are able to sequester cations around points of contact between lamellae that are fully dehydrated, while essentially maintaining a multilamellar structure for the entire system. This observation suggests that the undulations of the fatty acid lamellae control the diffusion of hydrophilic species through the water phase by altering the positional and rotational order of water molecules in the embedded/occluded "droplets."

pKa values in proteins determined by electrostatics applied to molecular dynamics trajectories

For a benchmark set of 194 measured pKa values in 13 proteins, electrostatic energy computations are performed in which pKa values are computed by solving the Poisson-Boltzmann equation. In contrast to the previous approach of Karlsberg(+) (KB(+)) that essentially used protein crystal structures with variations in their side chain conformations, the present approach (KB2(+)MD) uses protein conformations from four molecular dynamics (MD) simulations of 10 ns each. These MD simulations are performed with different specific but fixed protonation patterns, selected to sample the conformational space for the different protonation patterns faithfully. The root-mean-square deviation between computed and measured pKa values (pKa RMSD) is shown to be reduced from 1.17 pH units using KB(+) to 0.96 pH units using KB2(+)MD. The pKa RMSD can be further reduced to 0.79 pH units, if each conformation is energy-minimized with a dielectric constant of εmin = 4 prior to calculating the electrostatic energy. The electrostatic energy expressions upon which the computations are based have been reformulated such that they do not involve terms that mix protein and solvent environment contributions and no thermodynamic cycle is needed. As a consequence, conformations of the titratable residues can be treated independently in the protein and solvent environments. In addition, the energy terms used here avoid the so-called intrinsic pKa and can therefore be interpreted without reference to arbitrary protonation states and conformations.

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

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

Recent development and application of constant pH molecular dynamics

Solution pH is a critical environmental factor for chemical and biological processes. Over the last decade, significant efforts have been made in the development of constant pH molecular dynamics (pHMD) techniques for gaining detailed insights into pH-coupled dynamical phenomena. In this article we review the advancement of this field in the past five years, placing a special emphasis on the development of the all-atom continuous pHMD technique. We discuss various applications, including the prediction of pK_a shifts for proteins, nucleic acids and surfactant assemblies, elucidation of pH-dependent population shifts, protein-protein and protein-RNA binding, as well as the mechanisms of pH-dependent self-assembly and phase transitions of surfactants and peptides. We also discuss future directions for the further improvement of the pHMD techniques.

Study of pK values and effective dielectric constants of ionizable residues in pentapeptides and in staphylococcal nuclease (SNase) using a mean-field approach

The determination of pK values of amino acid residues as a function of temperature and ionic concentration is crucial to understanding the dynamics of various biological processes such as adsorption of peptides and their interactions with active sites of enzymes. In this study we developed a mean-field model to calculate the position-dependent dielectric constants of ionizable groups and the mean electrostatic potential on the surface. Such potential, which takes into account the contributions exerted by neighboring groups and ions in solution, is responsible for the fine-tuning of the pK value of each residue. The proposed model was applied to the amino acids Asp, Glu, Lys, His, Tyr, and Cys, and since the results were consistent with experimentally obtained values, the model was extended and applied to computation of pK values of Gly and Ala pentapeptides and of ionizable residues of the enzyme staphylococcal nuclease (SNase). In this latter case, we used an approach similar to a first-neighbors approximation, and the results turned out to be in good agreement with previously reported data when considering only the interactions of charged groups located at distances of maximally 20 Å. These considerations and the little computational cost involved turn the suggested approach into a promising tool for the modeling of force fields in computational simulations.

The transition from the native to the acid-state characterized by multi-spectroscopy approach: study for the holo-form of bovine α-lactalbumin

The transition of the holo-form of bovine α-lactalbumin from the native (N) to the pH-generated acidic-state (A-state) was analyzed by probing its tertiary and secondary structure using a concerted spectroscopic approach combining near- and far-UV circular dichroism (CD), electrospray ionization ion mobility mass spectrometry (ESI-IM-MS), vibrational circular dichroism (VCD), and Fourier transform infrared spectroscopy (FTIR) in the attenuated total reflection (ATR) and transmission (TR) modes. The spectroscopic results, which relied on the interaction of an electromagnetic field with different molecular targets, confirmed the decay of extensive rigid side-chain packing interactions during the pH-induced N→A-state transition and revealed the targets' dependence on secondary structural changes. Independent analyses of the spectral changes using two methods of multivariate analysis, such as principal component analysis and two-dimensional correlation spectroscopy, revealed small but significant differences in the secondary structure as a result of the all-or-none transition. The cooperativity of the transition was quantitatively described using values corresponding to the mid-point (tm) and width of the transition (Δtm). The averages of the two parameters, calculated using the data collected by the different probes, were equal to 3.5±0.2 and 0.6±0.1(SE), respectively. The variable two-state nature of the cooperative N→A-state transition confirmed that the protonation of the side chain carboxyl groups on the Asp and Glu residues and that the release of a Ca(2+) ion induced structural changes on both the secondary and tertiary levels. The changes have been confirmed by results obtained from the concerted spectroscopic approach.

BcsA and BcsB form the catalytically active core of bacterial cellulose synthase sufficient for in vitro cellulose synthesis

Cellulose is a linear extracellular polysaccharide. It is synthesized by membrane-embedded glycosyltransferases that processively polymerize UDP-activated glucose. Polymer synthesis is coupled to membrane translocation through a channel formed by the cellulose synthase. Although eukaryotic cellulose synthases function in macromolecular complexes containing several different enzyme isoforms, prokaryotic synthases associate with additional subunits to bridge the periplasm and the outer membrane. In bacteria, cellulose synthesis and translocation is catalyzed by the inner membrane-associated bacterial cellulose synthase (Bcs)A and BcsB subunits. Similar to alginate and poly-β-1,6 N-acetylglucosamine, bacterial cellulose is implicated in the formation of sessile bacterial communities, termed biofilms, and its synthesis is likewise stimulated by cyclic-di-GMP. Biochemical studies of exopolysaccharide synthesis are hampered by difficulties in purifying and reconstituting functional enzymes. We demonstrate robust in vitro cellulose synthesis reconstituted from purified BcsA and BcsB proteins from Rhodobacter sphaeroides. Although BcsA is the catalytically active subunit, the membrane-anchored BcsB subunit is essential for catalysis. The purified BcsA-B complex produces cellulose chains of a degree of polymerization in the range 200-300. Catalytic activity critically depends on the presence of the allosteric activator cyclic-di-GMP, but is independent of lipid-linked reactants. Our data reveal feedback inhibition of cellulose synthase by UDP but not by the accumulating cellulose polymer and highlight the strict substrate specificity of cellulose synthase for UDP-glucose. A truncation analysis of BcsB localizes the region required for activity of BcsA within its C-terminal membrane-associated domain. The reconstituted reaction provides a foundation for the synthesis of biofilm exopolysaccharides, as well as its activation by cyclic-di-GMP.

Monte Carlo simulations of proteins at constant pH with generalized Born solvent, flexible sidechains, and an effective dielectric boundary

Titratable residues determine the acid/base behavior of proteins, strongly influencing their function; in addition, proton binding is a valuable reporter on electrostatic interactions. We describe a method for pK(a) calculations, using constant-pH Monte Carlo (MC) simulations to explore the space of sidechain conformations and protonation states, with an efficient and accurate generalized Born model (GB) for the solvent effects. To overcome the many-body dependency of the GB model, we use a "Native Environment" approximation, whose accuracy is shown to be good. It allows the precalculation and storage of interactions between all sidechain pairs, a strategy borrowed from computational protein design, which makes the MC simulations themselves very fast. The method is tested for 12 proteins and 167 titratable sidechains. It gives an rms error of 1.1 pH units, similar to the trivial "Null" model. The only adjustable parameter is the protein dielectric constant. The best accuracy is achieved for values between 4 and 8, a range that is physically plausible for a protein interior. For sidechains with large pKa shifts, ≥2, the rms error is 1.6, compared to 2.5 with the Null model and 1.5 with the empirical PROPKA method.

Protonation of trimethylamine N-oxide (TMAO) is required for stabilization of RNA tertiary structure

The osmolyte trimethylamine N-oxide (TMAO) stabilizes the tertiary but not the secondary structures of RNA. However, molecular dynamics simulations performed on the PreQ1 riboswitch showed that TMAO destabilizes the tertiary riboswitch structure, leading us to hypothesize that the presence of RNA could result in enhanced population of the protonated form, TMAOP. Constant pH replica exchange simulations showed that a percentage of TMAO is indeed protonated, thus contributing to the stability of the tertiary but not the secondary structure of PreQ1. TMAOP results in an unfavorable dehydration of phosphodiester backbone, which is compensated by electrostatic attraction between TMAOP and the phosphate groups. In addition, TMAOP interacts with specific sites in the tertiary RNA structure, mimicking the behavior of positively charged ions and of the PreQ1 ligand in stabilizing RNA. Finally, we predict that TMAO-induced stabilization of RNA tertiary structures should be strongly pH dependent.

Rational modification of protein stability by targeting surface sites leads to complicated results

The rational modification of protein stability is an important goal of protein design. Protein surface electrostatic interactions are not evolutionarily optimized for stability and are an attractive target for the rational redesign of proteins. We show that surface charge mutants can exert stabilizing effects in distinct and unanticipated ways, including ones that are not predicted by existing methods, even when only solvent-exposed sites are targeted. Individual mutation of three solvent-exposed lysines in the villin headpiece subdomain significantly stabilizes the protein, but the mechanism of stabilization is very different in each case. One mutation destabilizes native-state electrostatic interactions but has a larger destabilizing effect on the denatured state, a second removes the desolvation penalty paid by the charged residue, whereas the third introduces unanticipated native-state interactions but does not alter electrostatics. Our results show that even seemingly intuitive mutations can exert their effects through unforeseen and complex interactions.

Charge-leveling and proper treatment of long-range electrostatics in all-atom molecular dynamics at constant pH

Recent development of constant pH molecular dynamics (CpHMD) methods has offered promise for adding pH-stat in molecular dynamics simulations. However, until now the working pH molecular dynamics (pHMD) implementations are dependent in part or whole on implicit-solvent models. Here we show that proper treatment of long-range electrostatics and maintaining charge neutrality of the system are critical for extending the continuous pHMD framework to the all-atom representation. The former is achieved here by adding forces to titration coordinates due to long-range electrostatics based on the generalized reaction field method, while the latter is made possible by a charge-leveling technique that couples proton titration with simultaneous ionization or neutralization of a co-ion in solution. We test the new method using the pH-replica-exchange CpHMD simulations of a series of aliphatic dicarboxylic acids with varying carbon chain length. The average absolute deviation from the experimental pK(a) values is merely 0.18 units. The results show that accounting for the forces due to extended electrostatics removes the large random noise in propagating titration coordinates, while maintaining charge neutrality of the system improves the accuracy in the calculated electrostatic interaction between ionizable sites. Thus, we believe that the way is paved for realizing pH-controlled all-atom molecular dynamics in the near future.

pH-sensitive residues in the p19 RNA silencing suppressor protein from carnation Italian ringspot virus affect siRNA binding stability

Tombusviruses, such as Carnation Italian ringspot virus (CIRV), encode a protein homodimer called p19 that is capable of suppressing RNA silencing in their infected hosts by binding to and sequestering short-interfering RNA (siRNA) away from the RNA silencing pathway. P19 binding stability has been shown to be sensitive to changes in pH but the specific amino acid residues involved have remained unclear. Using constant pH molecular dynamics simulations, we have identified key pH-dependent residues that affect CIRV p19-siRNA binding stability at various pH ranges based on calculated changes in the free energy contribution from each titratable residue. At high pH, the deprotonation of Lys60, Lys67, Lys71, and Cys134 has the largest effect on the binding stability. Similarly, deprotonation of several acidic residues (Asp9, Glu12, Asp20, Glu35, and/or Glu41) at low pH results in a decrease in binding stability. At neutral pH, residues Glu17 and His132 provide a small increase in the binding stability and we find that the optimal pH range for siRNA binding is between 7.0 and 10.0. Overall, our findings further inform recent experiments and are in excellent agreement with data on the pH-dependent binding profile.

Atomistic simulations of pH-dependent self-assembly of micelle and bilayer from fatty acids

Detailed knowledge of the self-assembly and phase behavior of pH-sensitive surfactants has implications in areas such as targeted drug delivery. Here we present a study of the formation of micelle and bilayer from lauric acids using a state-of-the-art simulation technique, continuous constant pH molecular dynamics (CpHMD) with conformational sampling in explicit solvent and the pH-based replica-exchange protocol. We find that at high pH conditions a spherical micelle is formed, while at low pH conditions a bilayer is formed with a considerable degree of interdigitation. The mid-point of the phase transition is in good agreement with experiment. Preliminary investigation also reveals that the effect of counterions and salt screening shifts the transition mid-point and does not change the structure of the surfactant assembly. Based on these data we suggest that CpHMD simulations may be applied to computational design of surfactant-based nano devices in the future.

Chaperone activation by unfolding

Conditionally disordered proteins can alternate between highly ordered and less ordered configurations under physiological conditions. Whereas protein function is often associated with the ordered conformation, for some of these conditionally unstructured proteins, the opposite applies: Their activation is associated with their unfolding. An example is the small periplasmic chaperone HdeA, which is critical for the ability of enteric bacterial pathogens like Escherichia coli to survive passage through extremely acidic environments, such as the human stomach. At neutral pH, HdeA is a chaperone-inactive dimer. On a shift to low pH, however, HdeA monomerizes, partially unfolds, and becomes rapidly active in preventing the aggregation of substrate proteins. By mutating two aspartic acid residues predicted to be responsible for the pH-dependent monomerization of HdeA, we have succeeded in isolating an HdeA mutant that is active at neutral pH. We find this HdeA mutant to be substantially destabilized, partially unfolded, and mainly monomeric at near-neutral pH at a concentration at which it prevents aggregation of a substrate protein. These results provide convincing evidence for direct activation of a protein by partial unfolding.

Considering protonation as a posttranslational modification regulating protein structure and function

Posttranslational modification is an evolutionarily conserved mechanism for regulating protein activity, binding affinity, and stability. Compared with established posttranslational modifications such as phosphorylation or ubiquitination, posttranslational modification by protons within physiological pH ranges is a less recognized mechanism for regulating protein function. By changing the charge of amino acid side chains, posttranslational modification by protons can drive dynamic changes in protein conformation and function. Addition and removal of a proton is rapid and reversible and, in contrast to most other posttranslational modifications, does not require an enzyme. Signaling specificity is achieved by only a minority of sites in proteins titrating within the physiological pH range. Here, we examine the structural mechanisms and functional consequences of proton posttranslational modification of pH-sensing proteins regulating different cellular processes.

pH-dependent dynamics of complex RNA macromolecules

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

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.

Potential of mean force of water-proton bath and molecular dynamic simulation of proteins at constant pH

An advanced implicit solvent model of water-proton bath for protein simulations at constant pH is presented. The implicit water-proton bath model approximates the potential of mean force of a protein in water solvent in a presence of hydrogen ions. Accurate and fast computational implementation of the implicit water-proton bath model is developed using the continuum electrostatic Poisson equation model for calculation of ionization equilibrium and the corrected MSR6 generalized Born model for calculation of the electrostatic atom-atom interactions and forces. Molecular dynamics (MD) method for protein simulation in the potential of mean force of water-proton bath is developed and tested on three proteins. The model allows to run MD simulations of proteins at constant pH, to calculate pH-dependent properties and free energies of protein conformations. The obtained results indicate that the developed implicit model of water-proton bath provides an efficient way to study thermodynamics of biomolecular systems as a function of pH, pH-dependent ionization-conformation coupling, and proton transfer events.

Constant pH Molecular Dynamics Simulations of Nucleic Acids in Explicit Solvent

The nucleosides of adenine and cytosine have pKa values of 3.50 and 4.08, respectively, and are assumed to be unprotonated under physiological conditions. However, evidence from recent NMR and X-Ray crystallography studies has revealed the prevalence of protonated adenine and cytosine in RNA macromolecules. Such nucleotides with elevated pKa values may play a role in stabilizing RNA structure and participate in the mechanism of ribozyme catalysis. With the work presented here, we establish the framework and demonstrate the first constant pH MD simulations (CPHMD) for nucleic acids in explicit solvent in which the protonation state is coupled to the dynamical evolution of the RNA system via λ-dynamics. We adopt the new functional form λ(Nexp) for λ that was recently developed for Multi-Site λ-Dynamics (MSλD) and demonstrate good sampling characteristics in which rapid and frequent transitions between the protonated and unprotonated states at pH = pKa are achieved. Our calculated pKa values of simple nucleotides are in a good agreement with experimentally measured values, with a mean absolute error of 0.24 pKa units. This work demonstrates that CPHMD can be used as a powerful tool to investigate pH-dependent biological properties of RNA macromolecules.