Constant pH Molecular Dynamics Research Articles

Constant-pH Simulations with the Polarizable Atomic Multipole AMOEBA Force Field.

Accurately predicting protein behavior across diverse pH environments remains a significant challenge in biomolecular simulations. Existing constant-pH molecular dynamics (CpHMD) algorithms are limited to fixed-charge force fields, hindering their application to biomolecular systems described by permanent atomic multipoles or induced dipoles. This work overcomes these limitations by introducing the first polarizable CpHMD algorithm in the context of the Atomic Multipole Optimized Energetics for Biomolecular Applications (AMOEBA) force field. Additionally, our implementation in the open-source Force Field X (FFX) software has the unique ability to handle titration state changes for crystalline systems including flexible support for all 230 space groups. The evaluation of constant-pH molecular dynamics (CpHMD) with the AMOEBA force field was performed on 11 crystalline peptide systems that span the titrating amino acids (Asp, Glu, His, Lys, and Cys). Titration states were correctly predicted for 15 out of the 16 amino acids present in the 11 systems, including for the coordination of Zn2+ by cysteines. The lone exception was for a HIS-ALA peptide where CpHMD predicted both neutral histidine tautomers to be equally populated, whereas the experimental model did not consider multiple conformers and diffraction data are unavailable for rerefinement. This work demonstrates the promise polarizable CpHMD simulations for pKa predictions, the study of biochemical mechanisms such as the catalytic triad of proteases, and for improved protein-ligand binding affinity accuracy in the context of pharmaceutical lead optimization.

Resolving coupled pH titrations using alchemical free energy calculations.

In a protein, nearby titratable sites can be coupled: the (de)protonation of one may affect the other. The degree of this interaction depends on several factors and can influence the measured . Here, we derive a formalism based on double free energy differences ( ) for quantifying the individual site values of coupled residues. As values can be obtained by means of alchemical free energy calculations, the presented approach allows for a convenient estimation of coupled residue s in practice. We demonstrate that our approach and a previously proposed microscopic formalism, can be combined with alchemical free energy calculations to resolve pH-dependent protein values. Toy models and both, regular and constant-pH molecular dynamics simulations, alongside experimental data, are used to validate this approach. Our results highlight the insights gleaned when coupling and microstate probabilities are analyzed and suggest extensions to more complex enzymatic contexts. Furthermore, we find that naïvely computed values that ignore coupling, can be significantly improved when coupling is accounted for, in some cases reducing the error by half. In short, alchemical free energy methods can resolve the values of both uncoupled and coupled residues.

Constant-pH molecular dynamics simulations of progressive activation in a proton-gated ion channel

Constant-pH molecular dynamics simulations of progressive activation in a proton-gated ion channel

How is Membrane Permeation of Small Ionizable Molecules Affected by Protonation Kinetics?

According to the pH-partition hypothesis, the aqueous solution adjacent to a membrane is a mixture of the ionization states of the permeating molecule at fixed Henderson-Hasselbalch concentrations, such that each state passes through the membrane in parallel with its own specific permeability. An alternative view, based on the assumption that the rate of switching ionization states is instantaneous, represents the permeation of ionizable molecules via an effective Boltzmann-weighted average potential (BWAP). Such an assumption is used in constant-pH molecular dynamics simulations. The inhomogeneous solubility-diffusion framework can be used to compute the pH-dependent membrane permeability for each of these two limiting treatments. With biased WTM-eABF molecular dynamics simulations, we computed the potential of mean force and diffusivity of each ionization state of two weakly basic small molecules: nicotine, an addictive drug, and varenicline, a therapeutic for treating nicotine addiction. At pH = 7, the BWAP effective permeability is greater than that determined by pH-partitioning by a factor of 2.5 for nicotine and 5 for varenicline. To assess the importance of ionization kinetics, we present a Smoluchowski master equation that includes explicitly the protonation and deprotonation processes coupled with the diffusive motion across the membrane. At pH = 7, the increase in permeability due to the explicit ionization kinetics is negligible for both nicotine and varenicline. This finding is reaffirmed by combined Brownian dynamics and Markov state model simulations for estimating the permeability of nicotine while allowing changes in its ionization state. We conclude that for these molecules the pH-partition hypothesis correctly captures the physics of the permeation process. The small free energy barriers for the permeation of nicotine and varenicline in their deprotonated neutral forms play a crucial role in establishing the validity of the pH-partitioning mechanism. Essentially, BWAP fails because ionization kinetics are too slow on the time scale of membrane crossing to affect the permeation of small ionizable molecules such as nicotine and varenicline. For the singly protonated state of nicotine, the computational results agree well with experimental measurements (P1 = 1.29 × 10-7 cm/s), but the agreement for neutral (P0 = 6.12 cm/s) and doubly protonated nicotine (P2 = 3.70 × 10-13 cm/s) is slightly worse, likely due to factors associated with the aqueous boundary layer (neutral form) or leaks through paracellular pathways (doubly protonated form).

Phbuilder: A Tool for Efficiently Setting up Constant pH Molecular Dynamics Simulations in GROMACS.

Constant pH molecular dynamics (MD) is a powerful technique that allows the protonation state of residues to change dynamically, thereby enabling the study of pH dependence in a manner that has not been possible before. Recently, a constant pH implementation was incorporated into the GROMACS MD package. Although this implementation provides good accuracy and performance, manual modification and the preparation of simulation input files are required, which can be complicated, tedious, and prone to errors. To simplify and automate the setup process, we present phbuilder, a tool that automatically prepares constant pH MD simulations for GROMACS by modifying the input structure and topology as well as generating the necessary parameter files. phbuilder can prepare constant pH simulations from both initial structures and existing simulation systems, and it also provides functionality for performing titrations and single-site parametrizations of new titratable group types. The tool is freely available at www.gitlab.com/gromacs-constantph. We anticipate that phbuilder will make constant pH simulations easier to set up, thereby making them more accessible to the GROMACS user community.

Modeling pH-Dependent Biomolecular Photochemistry.

The tuning mechanism of pH can be extremely challenging to model computationally in complex biological systems, especially with respect to the photochemical properties. This article reports a protocol aimed at modeling pH-dependent photodynamics using a combination of constant-pH molecular dynamics and semiclassical nonadiabatic molecular dynamics simulations. With retinal photoisomerization in Anabaena sensory rhodopsin (ASR) as a testbed, we show that our protocol produces pH-dependent photochemical properties, such as the isomerization quantum yield or decay rates. We decompose our results into single-titrated residue contributions, identifying some key tuning amino acids. Additionally, we assess the validity of the single protonation state picture to represent the system at a given pH and propose the most populated protein charge state as a compromise between cost and accuracy.

Probing the pH Sensitivity of OprM: Insights into Metastable States and Semi-Open Conformation.

Efflux pumps are specialized transport proteins that play a key role in the bacterial defense against a wide spectrum of antibiotics. Hence, understanding the biophysical mechanism associated with this complex system of drug expulsion becomes crucial. This work deals with some vital aspects of the outer membrane factor (OMF) of MexAB-OprM. After being passed through MexB and MexA, efflux substrates have to go through OprM for their final judgment. Thus, it is very important to understand the periplasmic pore opening mechanism and the associated biophysical changes during this process. Our study captures a detailed analysis of the pore opening mechanism involving OprM. With powerful molecular dynamics (MD) techniques such as well-tempered metadynamics, the presence of metastable states in between open and closed states was confirmed. Also, upon mutating R376, the energy barrier for the conversion of the close to open conformation decreases, indicating an important role played by the residue. Further, constant pH MD was performed to capture the effect of pH in both conformations. OprM exhibits distinct conformational states at pH values greater than 5.5 and lower than 5.5, suggesting its pH-responsive characteristics. Overall, our study elucidates a crucial undertaking toward discovering potential inhibitors for MexAB-OprM efflux pumps.

Understanding Cysteine Reactivity in Protein Environments with Electric Fields.

The role cysteine residues play in proteins is mediated by their protonation state, whereby the thiolate form of the side chain is highly reactive while the thiol form is more inert. However, the pKa of cysteine residues is hard to predict as it can differ widely from its reference value in solution, an effect that is accentuated by local effects in the heterogeneous protein environment. Here, we present a new approach to the prediction of cysteine reactivity based on electric field calculations at the thiol/thiolate group. We validated our approach by predicting the protonation state of cysteine residues in different protein environments (in the active site, at the protein surface, and buried within the protein interior), including Cys-25 in papaya protease omega, which was proven problematic for the more traditional constant pH molecular dynamics (MD) technique. We predict pKa shifts consistent with experimental observations, and the decomposition of the electric fields into contributions from molecular fragments provides a direct handle to rationalize local pH and pKa effects in proteins without introducing parameters other than those of the force field used for MD simulations.

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.

Enhanced Uptake of Anticancer C6-Pep Dimer in a Model Membrane Caused by Differential pKa in Acidic Microenvironment.

Acidic tumor microenvironment (TME) presents a challenge for the action of antitumor drugs by acting as an additional barrier for the passive crossing of the cell membrane by chemotherapic agents playing a critical role in the proliferation of tumor cells. Anticancer lipopeptide C6-Pep dimer containing the leucine zipper motif shows an increased uptake into the model tumor membrane in TME, and application of external heat might lead to the uncoiling of the zipper, which could result in cell lysis. This work investigated the cause of this increased uptake of C6-Pep dimer into the bilayer model in TME. Accurate protonation states of all the titratable residues of the C6-Pep dimer in TME were determined using constant pH molecular dynamics. In TME, except for two terminal Glu5 residues, all other Glu residues in the C6-Pep dimer were permanently protonated. The remaining Glu5 residues had differential pKa values, leading to the construction of four possible dimers with different fixed protonation states, and molecular dynamics was used to study their interaction with the anionic bilayer. Except for the dimer at a physiological pH, the other dimers were positively charged and could readily adsorb on the membrane surface. The free energy of insertion of these dimers in the bilayer was lower for single and double protonated Glu5-containing dimers than for the others. After the insertion of the lipopeptides into the membrane, thinning of the bilayer in the vicinity of dimers and an increase in area per lipid of the bilayer were observed for all systems, indicating destabilization of the bilayer due to this intercalation. This study shows that the anticancer lipopeptide C6-Pep utilizes the TME around a tumor cell for insertion into the membrane.

Histidine Protonation and Conformational Switching in Diphtheria Toxin Translocation Domain.

Protonation of key histidine residues has been long implicated in the acid-mediated cellular action of the diphtheria toxin translocation (T-) domain, responsible for the delivery of the catalytic domain into the cell. Here, we use a combination of computational (constant-pH Molecular Dynamics simulations) and experimental (NMR, circular dichroism, and fluorescence spectroscopy along with the X-ray crystallography) approaches to characterize the initial stages of conformational change happening in solution in the wild-type T-domain and in the H223Q/H257Q double mutant. This replacement suppresses the acid-induced transition, resulting in the retention of a more stable protein structure in solutions at pH 5.5 and, consequently, in reduced membrane-disrupting activity. Here, for the first time, we report the pKa values of the histidine residues of the T-domain, measured by NMR-monitored pH titrations. Most peaks in the histidine side chain spectral region are titrated with pKas ranging from 6.2 to 6.8. However, the two most up-field peaks display little change down to pH 6, which is a limiting pH for this protein in solution at concentrations required for NMR. These peaks are absent in the double mutant, suggesting they belong to H223 and H257. The constant-pH simulations indicate that for the T-domain in solution, the pKa values for histidine residues range from 3.0 to 6.5, with those most difficult to protonate being H251 and H257. Taken together, our experimental and computational data demonstrate that previously suggested cooperative protonation of all six histidines in the T-domain does not occur.

Integrative Approach to Dissect the Drug Resistance Mechanism of the H172Y Mutation of SARS-CoV-2 Main Protease.

Nirmatrelvir is an orally available inhibitor of SARS-CoV-2 main protease (Mpro) and the main ingredient of Paxlovid, a drug approved by the U.S. Food and Drug Administration for high-risk COVID-19 patients. Recently, a rare natural mutation, H172Y, was found to significantly reduce nirmatrelvir's inhibitory activity. As the COVID-19 cases skyrocket in China and the selective pressure of antiviral therapy builds in the US, there is an urgent need to characterize and understand how the H172Y mutation confers drug resistance. Here, we investigated the H172Y Mpro's conformational dynamics, folding stability, catalytic efficiency, and inhibitory activity using all-atom constant pH and fixed-charge molecular dynamics simulations, alchemical and empirical free energy calculations, artificial neural networks, and biochemical experiments. Our data suggest that the mutation significantly weakens the S1 pocket interactions with the N-terminus and perturbs the conformation of the oxyanion loop, leading to a decrease in the thermal stability and catalytic efficiency. Importantly, the perturbed S1 pocket dynamics weaken the nirmatrelvir binding in the P1 position, which explains the decreased inhibitory activity of nirmatrelvir. Our work demonstrates the predictive power of the combined simulation and artificial intelligence approaches, and together with biochemical experiments, they can be used to actively surveil continually emerging mutations of SARS-CoV-2 Mpro and assist the optimization of antiviral drugs. The presented approach, in general, can be applied to characterize mutation effects on any protein drug targets.

Basis for Accurate Protein pKa Prediction with Machine Learning.

pH regulates protein structures and the associated functions in many biological processes via protonation and deprotonation of ionizable side chains where the titration equilibria are determined by pKa's. To accelerate pH-dependent molecular mechanism research in the life sciences or industrial protein and drug designs, fast and accurate pKa prediction is crucial. Here we present a theoretical pKa data set PHMD549, which was successfully applied to four distinct machine learning methods, including DeepKa, which was proposed in our previous work. To reach a valid comparison, EXP67S was selected as the test set. Encouragingly, DeepKa was improved significantly and outperforms other state-of-the-art methods, except for the constant-pH molecular dynamics, which was utilized to create PHMD549. More importantly, DeepKa reproduced experimental pKa orders of acidic dyads in five enzyme catalytic sites. Apart from structural proteins, DeepKa was found applicable to intrinsically disordered peptides. Further, in combination with solvent exposures, it is revealed that DeepKa offers the most accurate prediction under the challenging circumstance that hydrogen bonding or salt bridge interaction is partly compensated by desolvation for a buried side chain. Finally, our benchmark data qualify PHMD549 and EXP67S as the basis for future developments of protein pKa prediction tools driven by artificial intelligence. In addition, DeepKa built on PHMD549 has been proven an efficient protein pKa predictor and thus can be applied immediately to, for example, pKa database construction, protein design, drug discovery, and so on.

Analysis of the ERK Pathway Cysteinome for Targeted Covalent Inhibition of RAF and MEK Kinases

The ERK pathway is one of the most important signaling cascades involved in tumorigenesis. So far, eight noncovalent inhibitors of RAF and MEK kinases in the ERK pathway have been approved by the FDA for the treatment of cancers; however, their efficacies are limited due to various resistance mechanisms. There is an urgent need to develop novel targeted covalent inhibitors. Here we report a systematic study of the covalent ligandabilities of the ERK pathway kinases (ARAF, BRAF, CRAF, KSR1, KSR2, MEK1, MEK2, ERK1, and ERK2) using constant pH molecular dynamics titration and pocket analysis. Our data revealed that the hinge GK (gate keeper)+3 cysteine in RAF family kinases (ARAF, BRAF, CRAF, KSR1, and KSR2) and the back loop cysteine in MEK1 and MEK2 are reactive and ligandable. Structure analysis suggests that the type II inhibitors belvarafenib and GW5074 may be used as scaffolds for designing pan-RAF or CRAF-selective covalent inhibitors directed at the GK+3 cysteine, while the type III inhibitor cobimetinib may be modified to label the back loop cysteine in MEK1/2. The reactivities and ligandabilities of the remote cysteine in MEK1/2 and the DFG-1 cysteine in MEK1/2 and ERK1/2 are also discussed. Our work provides a starting point for medicinal chemists to design novel covalent inhibitors of the ERK pathway kinases. The computational protocol is general and can be applied to the systematic evaluation of covalent ligandabilities of the human cysteinome.

How the pH Controls Photoprotection in the Light-Harvesting Complex of Mosses.

In response to varying light conditions, light-harvesting complexes (LHCs) switch from a light-harvesting state to a quenched state to protect the photosynthetic organism from excessive light irradiation in a strategy known as nonphotochemical quenching (NPQ). NPQ is activated by an acidification of the thylakoid lumen, which is sensed directly or indirectly by the LHC, resulting in a conformational change of the complex that leads to the quenched state. The conformational changes responsible for NPQ activation and their connection to specific quenching mechanisms are still unknown. Here, we investigate the pH-triggered conformational changes in the light-harvesting complex stress-related (LHCSR) of mosses. By combining constant-pH molecular dynamics and enhanced sampling techniques, we find that the pH sensitivity of the complex is driven by the coupled protonation of three residues modulating the conformation of the short amphipathic helix placed at the lumen side of the embedding membrane. Combining these results with quantum mechanics/molecular mechanics calculations, we show that the quenching mechanism sensitive to the pH goes through a charge-transfer between a carotenoid and an excited chlorophyll, which is controlled by the protein conformation.

Constant pH molecular dynamics simulations of membrane proteins.

Environment-driven dynamic protonation of titratable residues in a protein plays a key role in membrane transport and pH response in various biological processes. The dominant protonation state of a titratable amino-acid residue depends on environmental conditions like pH, dielectric environment (ionic strength, proximity to water-membrane interfaces), other amino-acid residues in its vicinity, etc. Constant pH Molecular Dynamics (CpHMD) techniques enable the study of dynamic protonation events in a molecular simulation and can be regarded as an in-silico titration process. In CpHMD techniques, the protonation state of a titratable residue is allowed to propagate between the different ionization states. As the protonation states of different residues change, the charge of the total system fluctuates. In hybrid solvent-based schemes, where the changes in the protonation states are handled using an implicit solvent representation, these charge fluctuations do not pose a problem. However, in case of an explicit-solvent representation, where the electrostatics are handled using a particle mesh ewald (PME) scheme, the varying charge poses a problem in reliable evaluation of energetics and thereby the pKa of the titratable residues. A large aqueous fraction in the simulation box is known to mitigate the charge artifacts. However, compared to soluble proteins, membrane protein simulation systems have a large non-aqueous volume fraction, owing to the lipid membrane. In this study, we compare different CpHMD methodologies for estimation of pKa in membrane proteins. These methods include both an explicit solvent-based CpHMD as well as hybrid solvent-based CpHMD. We present the pKa values of the acidic residues (aspartic acid, glutamic acid) from these simulations. In general, the hybrid solvent-based scheme showed a faster convergence compared to the explicit solvent-based scheme. However, in case of buried residues, the hybrid solvent-based scheme took a longer time to converge.

The electronic structure of genome editors from the first principles

Ab-initio molecular dynamics enables following the dynamics of biological systems from the first principles, describing the electronic structure and offering the opportunity to “watch” the evolution of biochemical processes with unique resolution, beyond the capabilities of state-of-the-art experimental techniques. This article reports the role of first-principles (ab-initio) molecular dynamics (MD) in the CRISPR-Cas9 genome editing revolution, achieving a profound understanding of the enzymatic function and offering valuable insights for enzyme engineering. We introduce the methodologies and explain the use of ab-initio MD simulations to establish the two-metal dependent mechanism of DNA cleavage in the RuvC domain of the Cas9 enzyme, and how a second catalytic domain, HNH, cleaves the target DNA with the aid of a single metal ion. A detailed description of how ab-initio MD is combined with free-energy methods—i.e., thermodynamic integration and metadynamics—to break and form chemical bonds is given, explaining the use of these methods to determine the chemical landscape and establish the catalytic mechanism in CRISPR-Cas9. The critical role of classical methods is also discussed, explaining theory and application of constant pH MD simulations, used to accurately predict the catalytic residues’ protonation states. Overall, first-principles methods are shown to unravel the electronic structure and reveal the catalytic mechanism of the Cas9 enzyme, providing valuable insights that can serve for the design of genome editing tools with improved catalytic efficiency or controllable activity.

Different pKa Shifts of Internal GLU8 in Human β-Endorphin Amyloid Revealing a Coupling of Internal Ionization and Stepwise Fibril Disassembly.

As a functional amyloid, human β-endorphin amyloid fibril features a β-solenoid conformation and store peptide hormones within acidic secretory granules, which would be released into the blood through fibril disassembly when the cellular milieu pH increases from acidic to neutral level on exocytosis. To gain detailed atomic mechanism of β-endorphin amyloid fibrils' pH-responsive disassembly, we conduct constant pH molecular dynamics simulations to investigate the structural and dynamical properties of β-endorphin amyloid fibrils in experiencing the environmental pH changes. Our results demonstrate a clear pKa shift of the internal ionizable residue of GLU8, and this shift becomes even more pronounced when it is buried more deeply in the amyloid fibrils. The unusual pKa of GLU8 reveals that its protonation state changes from the protonated state in the acidic secretory granule to the deprotonated state in the neutral pH conditions in the blood, where the deprotonation of GLU8 leads to unfavorable interactions within the hydrophobic core of the amyloid and subsequent fibril disassembly. The different pKa shifts of GLU8 relative to its positions in the amyloid fibril indicate that the β-endorphin amyloid fibril disassembly is a stepwise process, accounting for the experimental observation that the disassembly always initiates from the outermost layer. This study reveals the critical role of the protonation state of GLU8 in amyloid fibrils' pH-responsive disassembly, and provides clear insights for understanding the structural transitions of amyloids in hormone secretion. This study also provides theoretical basis for designing pH-sensitive biological tools for specific use with precise positioning of ionizable residues into the hydrophobic interior of proteins.

The role of conformational change and key glutamic acid residues in the ClC-ec1 antiporter

The triple glutamine (Q) mutant (QQQ) structure of a Cl−/H+ antiporter from Escherichia coli (ClC-ec1) displaying a novel backbone arrangement has been used to challenge the long-held notion that Cl−/H+ antiporters do not operate through large conformational motions. The QQQ mutant substitutes the glutamine residue for an external glutamate E148, an internal glutamate E203, and a third glutamate E113 that hydrogen-bonds with E203. However, it is unknown if QQQ represents a physiologically relevant state, as well as how the protonation of the wild-type glutamates relates to the global dynamics. We herein apply continuous constant-pH molecular dynamics to investigate the H+-coupled dynamics of ClC-ec1. Although any large-scale conformational rearrangement upon acidification would be due to the accumulation of excess charge within the protein, protonation of the glutamates significantly impacts mainly the local structure and dynamics. Despite the fact that the extracellular pore enlarges at acidic pHs, an occluded ClC-ec1 within the active pH range of 3.5–7.5 requires a protonated E148 to facilitate extracellular Cl− release. E203 is also involved in the intracellular H+ transfer as an H+ acceptor. The water wire connection of E148 with the intracellular solution is regulated by the charge states of the E113/E203 dyad with coupled proton titration. However, the dynamics extracted from our simulations are not QQQ-like, indicating that the QQQ mutant does not represent the behavior of the wild-type ClC-ec1. These findings reinforce the necessity of having a protonatable residue at the E203 position in ClC-ec1 and suggest that a higher level of complexity exists for the intracellular H+ transfer in Cl−/H+ antiporters.

GPU-Accelerated All-Atom Particle-Mesh Ewald Continuous Constant pH Molecular Dynamics in Amber.

Constant pH molecular dynamics (MD) simulations sample protonation states on the fly according to the conformational environment and user specified pH conditions; however, the current accuracy is limited due to the use of implicit-solvent models or a hybrid solvent scheme. Here, we report the first GPU-accelerated implementation, parametrization, and validation of the all-atom continuous constant pH MD (CpHMD) method with particle-mesh Ewald (PME) electrostatics in the Amber22 pmemd.cuda engine. The titration parameters for Asp, Glu, His, Cys, and Lys were derived for the CHARMM c22 and Amber ff14sb and ff19sb force fields. We then evaluated the PME-CpHMD method using the asynchronous pH replica-exchange titration simulations with the c22 force field for six benchmark proteins, including BBL, hen egg white lysozyme (HEWL), staphylococcal nuclease (SNase), thioredoxin, ribonuclease A (RNaseA), and human muscle creatine kinase (HMCK). The root-mean-square deviation from the experimental pKa's of Asp, Glu, His, and Cys is 0.76 pH units, and the Pearson's correlation coefficient for the pKa shifts with respect to model values is 0.80. We demonstrated that a finite-size correction or much enlarged simulation box size can remove a systematic error of the calculated pKa's and improve agreement with experiment. Importantly, the simulations captured the relevant biology in several challenging cases, e.g., the titration order of the catalytic dyad Glu35/Asp52 in HEWL and the coupled residues Asp19/Asp21 in SNase, the large pKa upshift of the deeply buried catalytic Asp26 in thioredoxin, and the large pKa downshift of the deeply buried catalytic Cys283 in HMCK. We anticipate that PME-CpHMD will offer proper pH control to improve the accuracies of MD simulations and enable mechanistic studies of proton-coupled dynamical processes that are ubiquitous in biology but remain poorly understood due to the lack of experimental tools and limitation of current MD simulations.