{'fraction_non_alphanumeric': 0.036122775269251714, 'fraction_numerical': 0.025616248221093874, 'mean_word_length': 4.516913831540176, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
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 condition; however, the accuracy of most current applications is limited due to the use of implicit-solvent models or hybrid solvent scheme. Here we report the first GPU-accelerated implementation, parameterization, and validation of the all-atom continuous constant pH MD (CpHMD) method with particle-mesh Ewald (PME) electrostatics in the Am-ber22 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 replicaexchange 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.


Introduction
Accurate and efficient molecular modeling of proton-coupled dynamic processes is important, as biological functions and material properties often depend on protonation and deprotonation. For example, many secondary active membrane transporters utilize pH gradient and proton coupling to drive the conformation transitions for function. 1 Many enzymes have pH-dependent catalytic activities, e.g., the active site of SARS-CoV-2 main protease collapses upon protonation of a conserved histidine residue. 2,3 Well known examples of pH-dependent materials include aminopolysaccharide chitosan which self assembles into hydrogels in response to a small increase in solution pH. 4 The ability to model proton-coupled dynamic processes is also important for studying protein-ligand binding, as upon protein-ligand association, the protonation state of the protein and the ligand may change. 5,6 Unlike conventional molecular dynamics (MD) that assumes fixed protonation states, constant pH MD allows protonation states to evolve with time according to the conformational environment and a preset solution pH. Currently, the most applied constant pH approaches are based on λ dynamics and the hybrid Monte-Carlo (MC)/MD scheme. 7 The former [8][9][10][11][12] uses continuous titration coordinates to propagate protonation states based on an extended Lagrange approach called λ dynamics, 13 while the latter [14][15][16][17] combines MD with periodic MC sampling of discrete protonated and deprotonated states. The details of these techniques can be found in the recent reviews. 7,18,19 The early implementations of the continuous λ dynamics based and the discrete MC based constant pH methods made use of the implicit-solvent (e.g., generalized Born) models. 8,9,15 The use of implicit-solvent models significantly reduces system size and allows faster sampling of solute conformational states relative to simulations with explicit water models. 18 However, for many biologically relevant systems, e.g., transmembrane proteins (with heterogeneous dielectric environment), nucleic acids (highly charged), and protein-ligand and protein-protein bound states, implicit-solvent models are not sufficiently accurate. This has motivated the development of explicit-solvent constant pH methods, which include the hybrid-solvent and all-atom approaches.
In a hybrid-solvent constant pH scheme, the MC based protonation state sampling or λ dynamics based propagation of titration coordinates is conducted in implicit solvent, while MD is conducted in explicit water. The first hybrid-solvent constant pH method was developed by Baptista and Soares, who combined MD in explicit solvent with MC sampling based on Poisson-Boltzmann (PB) calculations. 14 Later, a hybrid-solvent continuous constant pH method was implemented in CHARMM, 20 and a hybrid-solvent discrete constant pH approach was implemented in Amber. 16 The hybrid-solvent constant pH methods demonstrated improved accuracy for conformational sampling and consequently better agreement with experimental pK a 's; 16,20 however, a Hamiltonian cannot be expressed in a hybrid scheme (semigrand canonical ensemble), and thus energy conservation is not proven to hold. In terms of applications, ligand description with an implicit solvent model is a challenge, and the effects of explicit water and ions cannot be modeled.
To overcome the limitations of the hybrid-solvent scheme, much effort has been made in the development of all-atom constant pH methods over the past decade. The CHARMM program 21 contains the CPU implementations of the all-atom continuous constant pH method with generalized reaction field 22,23 or particle-mesh Ewald (PME) electrostatics for λ dynamics. 11 The multiple site λ dynamics (MSλD) 24 based constant pH method with a long-range electrostatic cutoff 10,25 was also implemented in CHARMM. 21 These methods have been validated using pK a calculations for a number of proteins 10,11,23 as well as RNAs. 25 The λ dynamics based constant pH method was also implemented in the GROMACS program, 26 although only the single-site titratable model was considered and performance for proteins remains to be demonstrated. 12 The NAMD program 27 contains an implementation of the all-atom constant pH method based on the non-equilibrium MD-MC approach, 17 which overcomes the issue of MC moves not accepted due to a large energy variation resulting from a sudden change in protonation states as in the aforementioned hybrid MC/MD constant pH approach.
The aforementioned all-atom continuous constant pH methods 10,11,22,25 are promising; however, the implementations were for CPU computing, which limits the simulation time scale and size of the systems that can be studied. In this work, we report the development and validation of the GPU all-atom continuous constant pH method in the pmemd.cuda engine of Amber program (version 2022). 28 Following the discussion of the model parameterization and validation, we present the data from the pK a calculations of benchmark proteins, including BBL, HEWL, SNase, RNase A, BACE 1, thioredoxin, and HMCK. In addition to comparison to experimental pK a values, we will discuss the coupled titration of catalytic residues, pH-dependent response of solvent exposure, titration of deeply buried sites. Finally, we will briefly discuss finite-size effects, and project future directions.

Methods and Implementation
The all-atom PME-CpHMD method. In contrast to conventional MD, the protonation state of a titratable site is treated as a dynamic variable (λ) of the system in the continuous constant pH MD method (hereafter referred to CpHMD). 8,29 When λ = 0, the proton is present and fully interacts with its environment, and when λ = 1, it is treated as a ghost particle without non-bonded interactions with its environment. The partial charges on the titratable residue are linearly scaled between those appropriate for a protonated state and those appropriate for a deprotonated state. 8,29 This differs from the multi-site λ dynamics 24 dynamics based MSλD CpHMD method, 10,25 which scales potential energies. Formally, the λ-dependent Lennard-Jones interaction energy between a titratable hydrogen i and another non-titratable atom j is given bỹ where λ i is the titration variable associated with the titratable hydrogen, and U LJ ij is the Lennard-Jones interaction energy between atoms i and j when the hydrogen is present.
Similarly, the Lennard-Jones interaction energy between two titratable hydrogens is given The charge of atom j in the titratable residue α is given by where q j (0) is the charge appropriate to the protonated form of the residue, and q j (1) is the charge appropriate to the deprotonated form.
The implementation presented in this paper uses fully explicit water molecules and treats the nonbonded electrostatic interaction energy between atoms with particle-mesh Ewald (PME) electrostatics. Because the λ values are treated as dynamic coordinates of the system, the derivatives of the energy with respect to the λ values are required. In the pmemd implementation 30 of PME electrostatics, this interaction energy is separated into several terms, where is the short-range component of the electrostatic energy, where n enumerates the copies of each atom from the neighboring periodic cells. This summation is performed only over those atom pairs i, j for which r ij falls within a small cutoff distance.
is the reciprocal space energy, where ν is the volume of the unit cell, m is reciprocal lattice vector, and S (m) is the structure factor, which can be approximated by where Q (k 1 , k 2 , k3) is a matrix constructed by interpolating the charge distribution in the simulation cell to a grid with the same dimensions k 1 , k 2 , k 3 , and F (Q) (m 1 , m 2 , m 3 ) is the fast Fourier transform of the Q matrix.
The V reciprocal can then be written as is a term that removes the self-interaction energies contained in V reciprocal , and where V is the volume of the unit cell is a term that counterbalances any net charge on the system. Implementation in the pmemd.cuda engine. As in our previous CPU implementation of the PME-CpHMD method in CHARMM, 11 the derivatives ofŨ LJ ij (λ i ) with respect to the titration variables can be written as Eqs. 1 and 2, which requires changes to be made to the Lennard-Jones forces between titratable atoms. In the present implementation, these modifications were made by making appropriate changes to the direct-force CUDA kernel in pmemd.cuda where the Lennard-Jones forces are computed. This kernel was also modified to compute the Lennard-Jones contributions to the forces on the λ titration variables. The electrostatic spatial forces on the atoms can be made to depend on the λ values by using the normal force calculations with the charges given in Eq. 3 according to the instantaneous values of the λ titration variables. Since V self and V plasma are independent of the spatial coordinates of the atoms, they are not computed during standard MD runs in pmemd.cuda. However, because these energies do depend on the λ titration coordinates, their derivatives with respect to the titration coordinates are required in CpHMD. The calculation of these derivatives was added to the kernel that interpolates the λ-dependent atomic partial charges, which was previously implemented by us for the generalized Born based CpHMD method. 31 This kernel otherwise required minimal changes for the present implementation. The derivatives of V direct with respect to λ are computed through appropriate changes to the direct-force kernel in pmemd.cuda. The derivatives of V reciprocal with respect to λ are computed with a new kernel that computes the derivatives given by Eq. 9 using the same method as outlined in our previous CPU implementation in CHARMM. 11 Modification of the force field parameters. The current constant pH methods are based on single topology, i.e., titration is represented by switching on/off the charge and Lenard-Jones interactions of the dummy hydrogen as well as by transforming between the protonated and deprotonated forms of sidechain partial charges. 7 The latter is straightforward to implement for the CHARMM force fields, 32 as the backbone partial charges are independent of the side chain. This is however not the case for the Amber force fields, 33 in which the backbone charges are dependent on the side chain protonation state. Due to the 1-4 interactions between the backbone and adjacent sidechain, this dependency presents a problem for constant pH simulations, as a single reference or model is used.
To circumvent this problem, we adopted a scheme used by Mongan et al. in the MC/MD based constant pH method 15 in which the backbone charges are fixed to the values of a one protonation state (in our case charged Asp/Glu and neutral His), and the residual change in the charge (ranging from 0.10 to 0.14 e for Asp, Glu and His) is absorbed onto the Cβ atom. In order to avoid introducing potential artifacts to conformational dynamics, such a scheme is only deployed for titration dynamics. For conformational dynamics, the partial charges are unmodified and the charge interpolation between the protonated and deprotonated states is made to both backbone and sidechain atoms.
Another compromise and approximation we made is in the bonded terms, which are not scaled between two protonation states as in the original CpHMD method. 8,29 For Asp and Glu, the bonded parameters for the protonated and deprotonated forms are different in all force fields. 32,33 The parameters of the deprotonated forms were used except for those related to the dummy hydrogens, which were taken from the protonated forms. For His, the bonded parameters for the protonated and deprotonated forms are the same in the Amber ff14SB 33 and CHARMM c22 32 force fields.
Finally, the use of two dummy hydrogens for Asp/Glu introduces a problem, namely, once an uncharged (ghost) dummy hydrogen rotates to the anti configuration, it loses the ability to titrate. Following our previous work, 29 the rotation barrier of the C-O bond in the carboxylate group of Asp/Glu was increased to 6 kcal/mol to keep the dummy hydrogens in the syn configuration. This is a limitation, as the anti configuration might become favorable in the protein although it is unfavorable in the peptide. 34 Potential of mean functions of titrating model peptides The linear response theory states that the charging free energy of an ion in polar solvent is quadratic in the charge perturbation. 35 Thus, the potential mean force (PMF) for protonation/deprotonation of a single-site titratable group (e.g., Cys and Lys) in explicit solvent can be approximated as a quadratic function in terms of λ. 11 Analogously, for a double-site titratable group, the model PMF in explicit solvent can be approximated as a quadratic function in terms of either λ or x, where x is a variable (x = sin 2 θ x ) that represents tautomeric interconversion. 11 For carboxylic acids Asp and Glu which have two equivalent protonation sites (carboxyl oxygens), the model PMF function can be written as 9 29,36 U mod where R 1 , ..., R 6 are parameters that can be determined by one-dimensional fitting of the corresponding mean forces (∂U/∂θ)| θx and (∂U/∂θ x )| θ calculated using thermodynamic integration (TI) at different combinations of θ and θ x values. The detailed derivation and protocol are given in Ref. 9,36 The model PMF function for His titration can be written as 9 The parameters in Eq. 14 are those in the one-dimensional PMF functions, where either λ or x is fixed at one of the end points (1 or 0). 9 Detailed protocols for obtaining the parameters are given in Ref. 9,36 3 Simulation Protocols Preparation of model peptide systems. Capped pentapeptides ACE-AAXAA-NH 2 (X = Asp, Glu, His, Cys, or Lys) were used to parameterize and validate the model PMF functions. First, each peptide structure was generated and placed in a cubic water box using CHARMM scripts (version c38b2). 21 The minimum distance between the heavy atoms of the peptide and the edges of the box was set as 10 Å. Next, to neutralize the system at pH 7.5, one Cl − counterion was added to the Lys pentapeptide system, and one Na + counterion was added to the Asp and Glu pentapeptide systems. The peptides were represented by the CHARMM c22, 32 Amber ff14SB, 33 or Amber ff19SB 37 force field.
Water was represented by the TIP3P water model. 38 Thermodynamic integration and titration simulations of the model peptides. We carried out an energy minimization in each pentapeptides system applying a force constant of 100 kcal mol −1 Å −2 to the peptide heavy atoms for 200 steps of SD followed by 300 steps of conjugate gradient method. Then, the system was heated from 100 to 300K using Langevin thermostat and a force constant of 5 kcal mol −1 Å −2 on the heavy atoms.
After heating, three stages of equilibration were performed with 250 ps each, whereby the ns. The TI simulations gave the mean forces, ∂U/∂θ| θx and ∂U/∂θ x | θ , which were used to obtain parameters in the PMF functions (Eq. 12, 13, and 14). The detailed protocols are given in a recent tutorial. 36 As validation of the model parameters, titration simulations were conducted for the model peptides at independent pH conditions, which were placed at 0.5-pH intervals in the range of 2-5.5 for Asp, 2.5-6 for Glu, 4.5-8 for His, 6.5-10 for Cys, and 8-11.5 for Lys model peptides. The equilibration and production runs of the peptide systems followed the same protocols as the protein simulations (see latter discussion). The production run at each pH lasted 20 ns and was repeated three times. With the CHARMM c22 force field, 32 we also performed pH replica-exchange simulations of 10 ns/replica for Asp, His, and Lys model peptides with the same pH conditions. Additional pH replica-exchange simulations were also performed with the hydrogen mass repartition scheme 39 (ribonucleas A or RNaseA), 44 1ERU (thioredoxin), 45 and 1I0E (human muscle creatine kinase or HMCK). 46 The coordinates were first processed using the convpdb.pl script from MMTSB Toolset 47 to remove hetero atoms, ions, water, ligands, and hydrogen atoms. The CHARMM c22 protein force field and CHARMM modified TIP3P water model were used to represent the protein and water, respectively. 32 The following steps were performed using the CHARMM package (c38b2). 21 The proteins were embedded in a pre-equilibrated cubic TIP3P water box with at least 10 Å cushion between the protein heavy atoms and the edges of the box. Sodium and chloride ions were added to neutralize the systems (assuming model pK a 's and pH 7.5) and to provide a physiological (0.15 M) or experimental salt concentration (0.1 M for SNase, 0.5 M for thioredoxin, and 0.06 M for RNase A). Using the HBUILD facility, missing hydrogens were added, and a custom CHARMM script is used to add two dummy hydrogens on the carboxylate oxygens. 9 The protein structures were energy minimized using 50 steps of steepest descent (SD) method with a harmonic force constant of 50 kcal/mol·Å 2 on the heavy atoms followed by 100 steps of adoptive basis Newton-Raphson (ABNR) method.
Equilibration of protein systems at independent pH conditions. The CHARMM22 topology and parameter files were converted to the Amber compatible format with the command chamber in ParmEd. 48 With the Amber input files prepared, a last round of minimization was performed in Amber22, 30 using 200 steps of SD followed by 300 steps of conjugate gradient method, whereby a force constant of 100 kcal mol −1 Å −2 was applied to the protein heavy atoms. Keeping the same restraint and with a time step of 1 fs, the system was then heated for 100 ps from the initial temperature of 100 K to 300 K using Langevin thermostat). Following heating, two stages of equilibration was performed. The

Production runs of protein systems with the pH replica-exchange protocol. For
CpHMD production runs, the asynchronous pH replica exchange algorithm 49 was employed to accelerate sampling convergence of conformational and protonation states and accelerate pK a calculations. 20 2 NVIDIA GTX 2080 Ti GPU cards were used. The pH range of the protein simulations was extended at least 1 pH unit below or above the lowest or highest experimental pK a values, and the pH spacing was 0.5 pH unit. Additional pH replica at 0.25 pH units were added in some cases to increase the probabilities of replica exchange. The exchanges between adjacent pH replicas were attempted every 2 ps (1000 MD steps). Each replica in the simulations of BBL, HEWL, SNase, thioredoxin, RNase A, and HMCK was run for 34, 40, 40, 50, 40, and 30 ns, respectively. The simulation length was sufficient to converge the pK a 's of all titratable sites (for HMCK we were only interested in Cys283). For SNase, additional simulations with larger box sizes were carried out. In these systems, the distance between the protein and edges of the where pK a and n are the fitting parameters, and S is defined where N unprot and N prot refer to the number of λ values representing the unprotonated and protonated states, respectively.
For two residues experiencing linked titration, the average number of protons bound to the two residues ( P ) are calculated at all simulation pH, and fit to the following coupled titration model to determine the macroscopic stepwise pK a 's, 22,50 P = 10 pK a2 −pH + 2 · 10 pK a1 +pK a2 −2pH 1 + 10 pK a2 −pH + 10 pK a1 +pK a2 −2pH (19) where pK a1 and pK a2 are the two stepwise pK a 's.
Finite-size corrections to the calculated pK a values for proteins. In our previous work, 11 we proposed a correction for pK a calculations for explicit-solvent constant pH simulations using lattice-sum methods under periodic boundary conditions. According to the analysis of Hünenberger and colleagues, the finite-size errors for the ligand charging free energies arise from four physical effects, among which the discrete solvent effect dominates when the protein's net charge is neutralized by counter-ions. 51 The discrete solvent effect arises from a homogeneous, constant potential that is applied to offset the potential generated by isotropically tumbling solvent molecules so that the average potential over the simulation box is zero. 51 This "offset" potential is positive for typical three-site water models, and needs to be corrected when calculating ligand charging free energies. 51 Hünenberger and colleagues developed an analytic correction to the ligand charging free energy 51 where k is the electrostatic constant, γ s is the quadrupole moment trace of the solvent model relative to a van der Waals interaction site. γ s is calculated as 0.764 e·Å 2 for TIP3P water model. 11 Q is the charge (-1 for charging to -1 e and +1 for charging to +1 e), N s is the number of solvent molecules, and V is the simulation box volume.
Now we consider the deprotonation reactions of protein titratable residues, which refers to the charging process of an acidic sidechain, e.g., aspartic acid, Asp −−→ Asp -, or the discharging process of a basic sidechain, e.g., histidine, His + −−→ His. Based on Eq. 20 (correction for the charging free energy), we obtain the correction for the deprotonation free energy of a titratable residue in a protein in reference to a model system where N s p /V p and N s m /V m refer to the solvent number density in the protein and model systems, respectively. Note, the minus sign in Eq. 20 and Q are absorbed due to the fact that Q is -1 for acidic residues, and for basic residues ∆G(deprot)=-∆G(charging). The corresponding pK a correction is where R is the ideal gas constant, T is the temperature. Since the solvent number density is higher in the model system than in the protein system, the pK a correction is negative for both acidic and basic sites.
Unless otherwise noted, the following pK a corrections were applied in the current data  Table S2.

Model parameterization and validation
Parameterization results. We obtained the parameters in the model PMF functions of Asp, Glu, and His (Eq. 13 and 14) for the CHARMM c22, 32 Amber ff14sb, 33 and ff19sb 37 force fields (Table S1). The model PMF parameters for Cys titration was also obtained for the CHARMM c22 32 and Amber ff14sb 33 force fields (Table S1). The fitting was generally very good (see an example fitting of His in Fig. 1), suggesting that the linear response theory holds, consistent with the results of PME-CpHMD in CHARMM. 11 We note, the parameters are more accurate when they are derived from fitting the mean forces ( ∂U/∂θ | θx and/or ∂U/∂θ x | θ ) instead of the potentials of mean forces (as in our previous work 11 ).
In the discussion that follows, we employed the CHARMM c22 force field. 52 A force field comparison study will be conducted in the near future.   Table 1). Trajectory analysis showed that the Asp and His sidechains can form hydrogen bonds (h-bonds) with the neighboring backbone group, resulting in meta-stable states. The carboxylate group of Asp is stabilized by h-bonding with the neighboring backbone amide, which contributes to the 0.3-unit underestimation of the target pK a value. This behavior was previously observed in both the GB and PME-CpHMD simulations in CHARMM. 9,11 To investigate if the proton-coupled conformational dynamics is adequately sampled for Asp and His in the independent pH simulations, we compared the results with those from three sets of pH replica-exchange simulations. The latter were conducted with the asynchronous pH replica-exchange scheme that was recently implemented for Amber simulations. 49 The previous work of us 11,20,31 and others 10,16 demonstrated that the pH replica-exchange protocol significant accelerates protonation state and conformational sampling. Indeed, the pK a convergence is significantly accelerated; the unprotonated  Table 2.
fractions generally plateau after about 5 ns, compared to more than 10 ns in the independent pH simulations ( Fig. S2 and S3). Interestingly, the resulting pK a 's (3.3 and 6.5) of Asp and His are very similar to those from the independent pH titration, which suggests that sampling is sufficient in the latter (Fig.S2 and S3). Note, to account for the difference between the calculated and target pK a of Asp pentapeptide, we changed the Asp reference pK a to 4.0 in the CpHMD parameter file (charmm22_pme.parm) for protein simulations.
In order to further accelerate simulations, we tested the sensitivities of pK a 's for 4-fs 10.0±0.01 10.4 a Independent pH simulations, whereby each simulation was conducted for 20 ns and repeated three times. b Three sets of pH replica-exchange simulations of 10 ns/replica. c Three sets of pH replica-exchange simulations of 10 ns/replica with the HMR scheme and 4-fs timestep. All pK a 's and errors were calculated from bootstrap. d Expt refers to the NMR derived pK a 's of the model pentapeptides from Thurlkill et al. 54 The His tautomer pK a 's are those estimated by Tanokura based on the NMR data of a model compound.  55 We note, evaluation of the 4-fs/HMR scheme for CpHMD simulations of proteins is not in the scope of the present work and will be conducted in the near future.

Titration simulations of proteins
Overall comparison of the calculated and experimental pK a values. To test the accuracy of the PME-CpHMD method for modeling protonation states of proteins, we calculated the pK a 's of Asp, Glu, His, and Cys residues in BBL, HEWL, SNase, RNase A, thioredoxin, and creatine kinase (HMCK) proteins, which have been previously used to benchmark CpHMD methods. 11,31,56,57 For a total of 67 residues, the root mean square error (RMSE) and the mean unsigned error (MUE) of the calculated pK a 's are 0.76 and 0.61, respectively, while the Pearson's correlation coefficient r is 0.85 (Figure 3). A more stringent test of the pK a prediction accuracy is to correlate the calculated and experimental pK a shifts (∆pK a ) with respect to model values, as the ∆pK a range is much smaller than the pK a range, exposing potentially problematic cases. Encouragingly, the r value for ∆pK a 's is 0.80, similar to the the r value for absolute pK a 's, suggesting that a good correlation with experimental is achieved and consistent for different residue types (see later discussion).

CHARMM.
The pK a 's of BBL, HEWL, and SNase have been previously calculated using the all-atom PME-CpHMD implementation in CHARMM (Table S4). 11 The r value of ∆pK a 's (from correlation with experiment) for these proteins from the present work is 0.80, which is nearly identical to that (0.78) using the CHARMM PME-CpHMD titration. 11 A comparison between the individual pK a 's shows that most pK a values agree within 0.2-0.3 units (Table S4); the agreement is especially remarkable for the pK a 's of the catalytic dyad in HEWL, which differ by 0.2 units for Glu35 and are identical for Asp52. 11 Note, the previous CHARMM PME-CpHMD simulations were run for 10 ns/replica, 11 whereas the present Amber PME-CpHMD simulations were run until full convergence for 30-40ns/replica. This analysis suggests that the pK a drift is small over time and the replica-exchange CpHMD simulations offer pK a calculations with good precision, consistent with our previous observations. 11 ,20 We further compared the calculated pK a 's of BBL and HEWL, which were previously reported with the MSλD method in CHARMM (Table   S4). Note, the MSλD uses a force-based cutoff for long-range electrostatics in λ dynamics and no finite-size correction for the pK a 's. 10 For BBL, the order of the two His pK a 's are in agreement with the experiment 58 and CpHMD results. Although the pK a 's from MSλD are 0.6-0.8 units higher, in better agreement with experiment, the simulation length was only 5 ns/replica and no finite-size correction was applied (which downshifts the pK a 's).
As to HEWL, the RMSE from MSλD 10 is nearly identical to the current work.  (c) Figure 3: Comparison between the calculated and experimental pK a 's and pK a shifts of the benchmark proteins. a) Calculated pK a 's vs. experimental pK a 's. b) Calculated vs. experimental pK a shifts with respect to the experimental model peptide pK a 's ( Table 1). The data for Asp, Glu, His, and Cys are shown in magenta, cyan, blue, and orange, respectively. Pearson's correlation coefficient (r) and RMSE are given. The solid black lines represent the linear regression. The shaded region indicates the calculated pK a 's within the overall RMSE (0.76 units) of the experimental values. To guide the eye, the dashed diagonal line (x=y) is shown. c) Histograms of the deviations between the calculated and experimental pK a 's for Asp (left), Glu (middle), and His (right) residues. .

Comparison of the calculated and experimental pK a 's of different residue types.
The Asp pK a 's vary the most in this dataset. Encompassing both down-and upshifted pK a 's, the experimental pK a range for Asp is 1.2 to 8.1, similar to the calculated range of 1.5 to 7.7 (Fig. 3a, magenta). The experimental pK a 's of Glu also display both the down-and upshifts, but the range is smaller than Asp, from 2.6 to 6.1, compared to the calculated range of 1.0 to 6.9 (Fig. 3a, cyan). The overall accuracy of the pK a calculation for Asp is slightly worse than Glu (RMSE of 0.76 and 0.69 respectively), but the r value for Asp (0.87) is somewhat larger than for Glu (0.72), which may be attributed to the larger pK a range. There are only 9 experimental pK a 's for His in the current dataset, which has a range of 5.2-6.5 and do not include upshifted values (Fig. 3a, blue). The calculated His pK a range is 4.2-6.6 ( Fig. 3a, blue), and there is a trend of systematic overestimation of pK a downshifts (Fig. 3c, right). In contrast, there is no clear trend for the pK a errors of Asp and Glu (Fig. 3c, left and  to the experimentally observed pK a downshift ( Table 2). The pK a downshift of His166 can be attributed to solvent exclusion and lack of hydrogen bonding (h-bonding) or electrostatic interactions (Fig. 4a). As pH decreases from 6 to 4, His166 undergoes a sigmoidal transition from the fully unprotonated neutral state to the fully protonated charged state (Fig. 4b). As expected, the fraction of the fully buried state also decreases (i.e., solvent exposure increases); however, the decrease does not appear to be sufficient, i.e., at low pH values the buried fraction doe not plateau (Fig. 4c). This analysis suggests that while the PME-CpHMD method is able to reproduce the experimental pK a downshift by capturing the pH-induced decrease in solvent exclusion or increase in solvent exposure of His166, sampling of the exposed state at low pH may be insufficient, which contributes to the pK a underestimation. Another potential contributor is an overestimation of the desolvation penalty by the CHARMM c22 force field. 32 Note, our previous CHARMM PME-CpHMD simulations gave a similarly underestimated pK a for His166 (by 1.3 units), 11  pK a calculation for HEWL: titration order of the catalytic dyad. HEWL is a small protein with 129 residues and 10 titratable sites; it is a popular test system for pK a prediction methods due to the abundance of experimental data. 60 The RMSE of the calculated pK a 's is 0.83 units. The two largest errors are for Asp52 and His15; the calculated pK a 's are 2 units over-and 1 unit underestimated, respectively ( Table 2). Despite the overestimation, the calculated pK a of Asp52 is 1.4 units lower than that of Glu35 (the second catalytic residue, Fig. 5), which indicates that Glu35 is a general acid and Asp52 is a general base in catalysis, in agreement with experiment (Table 2). Consistent with the CHARMM PME-CpHMD as well as the hybrid-solvent CpHMD simulations, 64 the titration events of Glu35 and Asp52 are uncoupled, as evident from the nearly identical stepwise pK a 's (7.0 and 5.5) from fitting to the two-proton titration model (Eq. 19; figure not shown). The lack of coupled titration is due to the relatively large distance between the carboxylate sidechains (>6.5 Å between the nearest carboxylate oxygens at any pH). Note, the calculated dyad pK a 's are nearly the same as the values from the CHARMM PME-CpHMD simulations. 11 The MSλD method in CHARMM gave a nearly identical pK a for Glu35, but a 1.1-unit lower pK a for Asp52. 10 To understand the pK a order of the catalytic dayd Glu35/Asp52 and the possible factors for the pK a overestimation of the latter, we compared the pH profiles of the sol-vent exposure, h-bonding and electrostatic interactions of the dyad residues with the pH-dependent titration curves (Fig. 5b-d). Both residues are partially buried. As Glu35 switches from being fully unprotonated to fully protonated in the pH range 7 to 10, the solvent exposure decreases from about 30% to just under 20% (Fig. 5b and c, blue). Asp52 has a similar behavior, except that the titration and change in solvent exposure are shifted to a lower pH range 7 to 4 ( Fig. 5b and c, red). Now we turn to h-bonding and electrostatic interactions that are also physical determinants of pK a shifts. 65 Glu35 does not form h-bonds below pH 7.5, and above pH 7.5, occasional salt-bridge interaction with Arg114 was observed, with an occupancy less than 20% (Fig. 5d, blue). In contrast, the carboxylate group of Asp52 can accept h-bonds from the sidechains of Asn46 and Asn59, and the occupancy increases to nearly 2 with increased deprotonation of Asp52 (Fig. 5d, red).
The analysis of solvent exposure and h-bond suggests that the latter is the major determinant for the lower pK a value of Asp52 relative to Glu35, in agreement with our previous work using the hybrid-solvent as well as the PME-CpHMD in CHARMM. 64 It is noteworthy that despite the correlation between charging of Asp52 and the increase in solvent exposure and h-bond formation, the pH profiles of solvent exposure and h-bond occupancy are more gradual than the titration curve, which may indicate that the pH-dependent conformational changes might be undersampled, contributing the overestimation of the pK a of Asp52.
pK a calculation for SNase: partially buried residues and coupled titration of Asp19 and Asp21. SNase has a large number of engineered mutants, which are popular model systems for testing pK a prediction methods. 65 Here we used the hyperstable, acidresistant form of SNase ∆+PHS (hereafter referred to as SNase) 43 which is only slightly (14 residues) larger than HEWL, but has 9 more titratable sites. The pK a 's of SNase are more challenging to predict than HEWL due to the fact that most titratable sites are partially buried. 20 The calculated pK a 's have a RMSE of 0.76 units, similar to the RMSE of 0.80 from the CHARMM PME-CpHMD simulations. 11 The largest error is for Glu129, for which the experimental pK a is 0.4 units down-and the calculated pK a is 0.8 units upshifted relative to the model value of 4.2. Curiously, simulation also fails to reproduce the direction of the experimental pK a shifts of Glu52, Glu57, and Glu101, although the magnitude of the errors is smaller (0.4, 0.7, and 0.4 respectively). Analysis showed that all these residues are partially buried, suggesting that desolvation penalty contributes to the pK a upshift. Based on the analysis of BBL's His166 and HEWL's Asp52, we hypothesized that simulation overestimates the desolvation penalty due to inadequate sampling of the solvent exposed state. To test this, we plotted the fractional SASA values vs. pH for Glu52, Glu57, Glu101, and Glu129 (Fig. S4). Deprotonation of glutamic acid is expected to induce larger solvent exposure. This is indeed the case for the more exposed residues Glu52 and Glu57 (fractional SASA about 60% at low pH), although the degree of increase is small. However, solvent exposure change with pH for the more buried residues Glu101 and Glu129, for which the fractional SASA values remain at about 40 and 20% throughout the entire pH range (Fig. S4). These data support the hypothesis that the solvent exposed state may be inadequately sampled, contributing to the desolvation related pK a upshift for Asp and Glu.
The NMR data 43 as well as our previous work 64 based on the hybrid-solvent and PME-CpHMD simulations in CHARMM suggest that the titration Asp19 and Asp21 is coupled.
The current simulations confirmed the strong coupling as a result of h-bond formation between the two residues (Fig. 6a). Fitting the titration data to a two-proton coupled equation (Eq. 19) gives the stepwise macroscopic pK a 's of 2.5 and 5.6 ( Fig. 6b), which are in good agreement with the experimental values of 2.2 and 6.5. 43 To assign the stepwise pK a 's to individual residues, we examine the pH-dependent probabilities of four microscopic states, doubly protonated (HH), singly protonated with proton on D19 (H−) or Asp21 (−H), and doubly deprotonated (−−) states (Fig. 6c). Above pH 7, Asp19/Asp21 are in the −− state (Fig. 6c, blue). As pH decreases from 7 to 5, the probability of −− decreases, while that of the −H or H− increases. Since the −H state (cyan) is more probable than the H− state (magenta) as protonation first occurs, Asp21 receives a proton first, which means the higher pK a should be assigned to Asp21. As pH further decreases from 5 to 2, both −H and H− states are possible; however, their combined probability decreases, while the probability of the HH state increases (Fig. 6c, red). Below pH 2, the latter state dominates. Analysis of h-bonding and electrostatic interactions (Fig. 6a) shows a network among Asp19, Asp21, Thr22, Thr41, and Arg35, consistent with the CHARMM hybrid-solvent and PME-CpHMD simulations. 31 pK a calculation for thioredoxin: the deeply buried Asp26 and Asp58. Thioredoxin has 105 residues with 18 titratable sites. The RMSE of the calculated pK a 's is 0.71. We first consider Asp26, which has one of the highest measured pK a 's of any carboxylic acids in proteins. Encouragingly, the calculated pK a of Asp26 is 7.7, in excellent agreement with the NMR-derived value of 8.1. The large pK a upshift of nearly 4 units can be attributed to the extremely low fraction of solvent exposure (below 7% at all pH, Fig. S5). Trajectory analysis showed that Asp26 does not form h-bonds with nearby residues. The only factor that may stabilize the deprotonated form is the salt-bridge formation with Lys39; however, the salt bridge is only formed above pH 8 and the solvent exposure is very low (<20% at pH 8, Fig. S5). A previous experimental study 61 suggested that the protonated Asp26 may be stabilized by donating a h-bond to the nearby Ser28; however, in the simulation the average distance from the hydroxyl oxygen of Ser28 to the nearest carboxylate oxygen of Asp26 is about 4.9 Å at pH 7, similar to the distance of 4.7 Å in the X-ray structure (PDB 1ERU). Thus, the dry environment, along with lack of polar interactions, results in the very large pK a upshift of Asp26.
The largest error in the calculated pK a 's of thioredoxin is for Asp58, whose direction of pK a shift is reproduced but the magnitude is 1.6 units too large (Table 2). Analysis showed that the Asp58 is also deeply buried, with ∼20% solvent exposure below pH 5, which explains the pK a upshift relative to the model. However, the solvent exposure only slightly increases to ∼30% at pH 8 before increasing steeply to over 50% at pH 10 ( Fig.   S6). H-bond analysis showed that the deprotonated Asp58 can accept h-bonds from the backbones of neighboring Asp60 and Asp61, which can stabilize the deprotonated state; however, the pH profile of the h-bond occupancy is irregular, showing a nearly 50% decreased occupancy in the pH range 4-8 (Fig. S6). The latter indicates a sampling issue, which explains the overestimation of the pK a of Asp58.
pK a calculation for RNase A: the deeply buried His12. RNase A has 124 residues with 14 titratable sidechains. The RMSE for the calculated pK a 's is 0.81, and the largest error is for His12 (Fig. 7a). The experimental pK a of His12 is 0.5 units downshifted relative to the model, and the simulation overestimated the downshift by 1.7 units (Table 2).
Analysis showed that His12 titrates over the pH range 3 to 6 (Fig. 7b), and the titration is correlated with two physical determinants, an increase in solvent exposure (decreased buried fraction) at lower pH (Fig. 7c) and an increase in h-bond formation at higher pH ( Fig. 7d). However, unlike in the previous GB-or hybrid-solvent CpHMD simulations, 7,31 the pH profiles of the buried fraction and the h-bond occupancy do not fully match the titration curve. Above pH 6, His12 is over 90% buried, and the decrease in the buried fraction at lower pH does not follow an expected sigmoidal curve. The buried fraction decreases by about 5% as pH decreases from 9 to 5 and remains constant between pH 4.5 and 2, before further decreasing to 80% at pH 0 (Fig. 7c). A major h-bond partner is the neighboring Asn11, which can donate a h-bond from its carboxamide group to the nitrogen of His12 to stabilize its deprotonated form (Fig. 7a). As pH increases from 3 to 6, the occupancy of the h-bond increases from zero to about 60%, and it further increases to nearly 100% at pH 8. Based on the above analysis, we suggest that the h-bond formation and solvent exposure of His12 may be insufficiently sampled in the pH range 3-6. The under-sampling of the solvent exposed state (buried fraction remains unchanged between pH 2 and 4.5) is particularly evident, which may be a major factor for the overestimation of pK a downshift of His12. . pK a calculation for HMCK: a buried active-site cysteine. To test the accuracy of Cys titration, we calculated the pK a of Cys283 in the active site of HMCK, which has a NMR measured pK a of 5.6, 63 one of the lowest in the literature. 66 Our simulations correctly reproduced the direction of the pK a shift relative to the model; however the downshift is 1.1 units underestimated compared to the experiment ( Table 2). Analysis showed that Cys283 is buried and its thiolate form can accept h-bonds from the sidechains and backbones of Ser285 and Asn286 ( Fig. 8a and b), consistent with the GB-based CpHMD titration simulation. 57 As Cys283 becomes deprotonated in the pH range 6 to 9, the total h-bond occupancy increases and plateaus at 1; however, the exposed fraction does not increase and instead remains at about 40% ( Fig. 8c and d). Since solvent exposure promotes the charged thiolate state and increases the pK a , we suggest that insufficient sampling of the solvent-exposed conformations may be a major contributor to the overestimated pK a of Cys283. Finite-size effect and corrections. Following the work of Hünenberger and colleagues, 51 we previously proposed an analytical pK a correction (Eq. 22) to correct for the effect of an offset potential introduced in PME simulations under periodic boundaries. 11 For the current simulations, the pK a corrections for Asp, Glu, His, and Cys are between -0.3 and -1.0 pH units (see Methods). To assess the effectiveness of the corrections and better understand the finite-size effect, we performed additional titration simulations of SNase with increased box sizes, i.e., adding more water to the simulation system. Table 3 summarizes the raw and corrected pK a 's using four different boxes, which have 10 (default), 12, 14, or 18 Å cushion space between the protein and edges of the box (minimum distance between the heavy atoms of protein and water oxygens on the box edges). The corresponding cubic box lengths are 68, 71, 76, and 84 Å, respectively.
We first examine the raw calculated pK a 's from simulations with different box sizes.
As expected, with increasing box size the raw pK a 's decrease for all but four residues (Fig. 9a). Increasing box size also leads to better agreement with the experimental pK a 's; the RMSEs of the raw pK a 's are 1.0, 1.0, 0.97 and 0.76 for boxes with 10, 12, 14, and 18 Å cushion space, respectively (Table 3 and Fig. 9a). The MUE also decreases from 0.86 (10 Å cushion) to 0.81 (12 Å cushion), 0.80 (14 Å cushion), and 0.62 (18 Å cushion) ( Table 3). Comparison of the raw pK a 's between the smallest (10 Å cushion) and largest (18 Å cushion) boxes shows that the pK a changes due to box size increase vary (Table 3, last column). Excluding the four residues (Asp19, Asp21, H121, and Glu135) that show very little pK a changes, the pK a 's mostly decrease by 0.3 to 0.7 units, as compared to the finite-size corrections of -0.70 to -0.8 units for the smallest box. The effect of box size is not clear for the coupled residues Asp19/Asp21, which have the raw calculated pK a 's of 3.2/6.3 with the smallest box; however, the pK a 's increase to 3.6/6.5 and 3.4/7.1 with the larger boxes (12 and 14 Å cushion space), and then decrease back to 3.1/6.4 with the largest box (18 Å cushion). Increasing box size has negligible effect on the downshifted pK a of His121. With the increasing box sizes, its raw pK a changes from 5.7 to 6.1, 5.7, and 5.8. Box size also shows little effect on the pK a of Glu135, which has the raw calculated pK a 's of 3.0, 3.4, 3.0, and 3.0 with the increasing box sizes. Now we examine the pK a corrections for the different simulation boxes. It is apparent that application of the finite-size correction removes the systematic overshift error (Fig. 9b). As the box size increases, the solvent number density increases and therefore the magnitude of the correction decreases (Eq. 21). The magnitude of the corrections decreases by about 0.4 pH units going from the smallest to the largest box. Interestingly, this difference is roughly the same as the average difference between the raw pK a 's (of all but the aforementioned four residues) calculated with the smallest and largest box (Table 3, last column), which suggests that the finite-size correction is valid. Another interesting observation is that the increasing box size does not significantly reduce the RMSE of the finite-size corrected pK a 's. The RMSE's are 0.76, 0.80, 0.80, and 0.70 with the increasing box sizes (Fig. 9b), which is another piece of evidence supporting the validity of the finite-size corrections.  Figure 9: Effect of box size on the calculated pK a 's. The errors of the raw (a) and finitesize corrected (b) pK a of SNase with different solvent cushion spaces, 10 (magenta), 12 (orange), 14 (green), and 18 Å (cyan). The corresponding RMSE values are shown next to the legends.  (10, 12, 14, and 18 Å) between the protein heavy atoms and edges of the water box, and the length of a cubic box (L box ) converted from the average system volume. The columns Raw and Corr refer to the pK a 's before and after the finite-size corrections (see Methods). The column ∆Box refers to the raw pK a difference between the largest and smallest boxes.

Concluding Discussion
We presented the first implementation, parameterization, and validation of the GPUaccelerated continuous constant pH particle-mesh Ewald molecular dynamics method in Amber22 (hereafter referred to as Amber PME-CpHMD). Titration parameters for three force fields (CHARMM c22, 32 Amber ff14SB, 33 and ff19SB 37 ) were derived and validated using model pentapeptides AAXAA, where X represents Asp, Glu, His, Cys, or Lys. To benchmark the performance and accuracy for constant pH simulations of proteins, we carried out titration simulations with the c22 force field for 6 proteins, including BBL, HEWL, SNase, RNase A, thioredoxin, and HMCK, which have NMR derived pK a values of Asp, Glu, His, and Cys residues. The asynchronous pH replica-exchange algorithm 49 was employed to enhance sampling of protonation and conformational states. The simulations were run for 30-50 ns per pH replica until all pK a 's were converged. The resulting RMSE and MUE with respect to the experimental pK a 's are 0.76 and 0.61, respectively, and the largest pK a deviation is 2 units. The Pearson's correlation coefficients for the calculated vs. experimental pK a 's and pK a shifts are 0.85 or 0.80, respectively. Importantly, the titration simulations quantitatively reproduced the experiment pK a orders of the catalytic dyad in HEWL and the coupled residues in SNase. Simulations also quantitatively captured one of the largest upshifted pK a 's of a deeply buried Asp in thioredoxin as well as the downshifted pK a of an active-site Cys in HMCK.
We compared the current validation data with those based on the CHARMM 21 CPU all-atom PME-CpHMD 11 and MSλD 10 simulations with the same c22 force field. 32 The Asp, Glu, and His pK a 's calculated from the CHARMM PME-CpHMD simulations of 10 ns per replica (much shorter than the present work) are in close agreement with the present work, suggesting that the pK a drifts over prolonged simulation time are small. Comparing to the calculated pK a 's of HEWL and the two His residues in BBL based on the MSλD simulations of 5-20 ns per pH replica (with a 12-Å electrostatic cutoff), 10 the overall RMSE is similar, and the pK a orders of the catalytic Glu35/Asp52 in HEWL and His142/His166 in BBL are consistent with the present simulations.
In agreement with the CHARMM PME-CpHMD simulations, 11 our data demonstrated that the finite-size effect needs to be taken into account for the accurate calculation of titration free energies with lattice sum methods under periodic boundary conditions. Following the work of Hünenberger and colleagues, 51 we previously proposed a pK a correction to account for a positive offset potential due to discrete solvent in periodic boxes, which is the dominant finite-size effect for all-atom PME constant pH simulations with a neutralized protein background. 11 Applying this correction to all calculated pK a 's, a systematic pK a upshift error was removed, and the overall agreement with experiment was improved. To further examine the finite-size effect and the validity of the correction, the pK a 's of SNase were calculated from simulations with four different box sizes. Consistent with the negative sign of the correction, increasing box size lowers the raw pK a 's of all but four residues that do not show significant changes. The RMSE of the raw pK a 's decreases from 1.0 with the smallest box to 0.76 with the largest box; the latter is identical to the RMSE (0.76) of the corrected pK a 's obtained from the simulation with the smallest box. The quantitative validity of the correction is also supported by a good agreement between the change in the finite-size correction and the average change of the raw pK a 's going from the smallest to the largest box size. As expected, the finite-size correction decreases with increasing box size, and consistently, the reduction in RMSE due to the correction also decreases.
Using the 18-Å cushion space, the correction is 0.3-0.4, and the RMSE (0.70) of the corrected pK a 's is only slightly smaller than the RMSE (0.76) of the raw pK a 's.
Although the overall box-size dependent trend of pK a 's is consistent with the positive offset potential being the dominant factor, 11 there are exceptions. The simulations of SNase showed that box size has negligible effect on the coupled pK a 's of Asp19 and Asp21 as well as the downshifted pK a 's of His121 and Glu135. We note that the effect of the offset potential and the corresponding pK a correction deal with an ideal situation in which a single residue titrates in a neutral background. Thus, it is possible that the correction is not valid for coupled pK a 's. However, with regards to the pK a 's of His121 and Glu135, the cause for the box size independence is difficult to speculate. Considering the above exceptions and the small magnitude of the correction and improvement of the RMSE using the 18 Å cushion space, we suggest that for simulations in large boxes (beyond 18 Å cushion space) a finite-size correction is not needed.
We analyzed the pH-dependent solvent exposure and formation of hydrogen bonds as well as electrostatic interactions of catalytic residues and those that exhibit larger pK a deviations from experiment. These analyses suggested while PME-CpHMD captures the proton-coupled conformational rearrangements, charging-induced increase of solvent exposure for buried residues is inadequate. This may be a major contributor to the pK a errors, including the overestimated pK a downshifts for buried His residues, e.g., His166 in BBL and His12 in RNaseA; the overestimated pK a upshifts for buried carboxyl residues, e.g., Glu57 in SNase and Asp58 in thioredoxin; and the overestimated pK a upshift for buried Cys, e.g., Cys283 in HMCK. Undersampling of the solvent-exposed state may also be related to the combination of c22/TIP3P force field, 32 which slightly biased solutesolute over solute-solvent interactions. 67 Overestimation of desolvation penalty may also be a source of error, which can be attributed to the low dielectric constant in the protein interior as a result of the lack of polarization in simulations with additive force fields. 68 Lack of polarization in the interior of protein may also lead to overly strong salt bridges, which may explain the overestimation of the pK a downshifts of Asp140 and Glu135 in SNase. While the use of polarizable force field for both protein and water is desirable, it may not be currently feasible due to speed. One intriguing idea worth exploring is to mix a polarizable water model such as OPC3-pol 69 with an additive force field to improve solute-solvent interactions. The present study did not examine the potential dependence on the additive force field. The force field related topics as well as the evaluation of PME-CpHMD for model proton-coupled conformational dynamics of catalytic residues in larger proteins (e.g., BACE1 56 ) will be explored in a future work.
By removing the reliance on the implicit-solvent model, the PME-CpHMD method can be applied to any system that has a force field representation. We anticipate the GPU accelerated PME-CpHMD to become a powerful tool for the investigation of a variety of proton-coupled dynamical phenomena that are poorly understood due to the current limitations in experimental and MD techniques, for example, secondary transport of ions/substrates across membrane transporter proteins and pH-dependent self-assembly of materials. Another important application of PME-CpHMD is to offer proper pH control, for example, by allowing protein and ligand to titrate while binding and unbinding, 5,6 or allowing His residues to fluctuate among the doubly protonated and two singly protonated tautomer states, which has been shown to affect the ligand binding mechanism and kinetics. 70,71 Supporting Information Available Supporting Information contains additional tables and figures.