solution ph often has 
 a dramatic impact on biomolecular systems . 
 by modulating the protonation state equilibria of various , titratable functional groups present in the system , small 
 changes in solution ph can affect the charge distribution within the 
 biomolecule . this charge distribution , in turn , often has a profound 
 impact on the fundamental structure and function of biomolecules . this effect can be so pronounced that some proteins native 
 states are stable only in a narrow ph range , even denaturing completely 
 in extreme ph environments . because biomolecular 
 behavior can depend very strongly on the ph - dependent 
 protonation states of various titratable residues , accurate computational 
 models designed to treat such systems must somehow account for ph 
 effects . while the traditional approach of assigning a fixed protonation 
 state for each titratable residue at the beginning of the simulation 
 is still the most common approach , numerous methods have been developed 
 in an attempt to treat ph effects in biomolecules more quantitatively . of particular interest in this study is 
 the constant ph molecular 
 dynamics ( cphmd ) technique , of which there are several variants . cphmd is a method that leverages the ability of classical molecular 
 dynamics to sample conformational space while simultaneously sampling 
 from the available protonation states according to the semigrand canonical 
 ensemble . by adopting this approach , cphmd 
 simulations can overcome the limitations imposed by constant - protonation 
 md simulations by constructing an ensemble whose protonation state 
 distributions are properly weighted via the thermodynamic constraint 
 imposed by a constant chemical potential of hydronium ions . within the myriad of available cphmd methods , there are two fundamentally 
 different approaches  continuous protonation states and discrete protonation states . in the former , a continuous titration particle 
 is introduced at each protonable site that is propagated as part of 
 the standard md according to a ph - dependent force acting on this particle . discrete protonation state models , on the other hand , employ occasional 
 metropolis monte carlo ( mc ) exchange attempts between different protonation 
 states throughout the course of the md simulation . for 
 the purposes of this study , we will focus on cphmd models using discrete 
 protonation states  specifically as implemented in the amber 
 software suite . in the original 
 amber implementation of cphmd , the molecular dynamics is propagated treating solvent effects implicitly via the generalized 
 born ( gb ) method . periodically throughout 
 the dynamics , a trial move changing the protonation state of either 
 one or two closely interacting titratable residues is evaluated based 
 on the difference in electrostatic and solvation energies calculated 
 via gb , after which the coordinates are propagated according to this 
 same potential . recently , however , machuqueiro and baptista 
 raised concerns about 
 pka predictions inheriting problems related 
 to the model compound definition and inaccuracies in the underlying 
 force field . in particular , force field 
 deficiencies have been shown to result in incorrect  even unphysical so far , reported applications of mongan et al.s method have 
 shown good results because the simulations were too short to reveal 
 the full extent of the weaknesses in the generalized born model being 
 used . when we implemented mongan et al.s method on graphics processing 
 units ( gpus ) in order to run long simulations , the unphysical states 
 sampled over long time scales degraded the pka predictions of the hen egg - white lysozyme ( hewl ) . this erroneous 
 sampling may be addressed to some extent by using an explicit representation 
 of the solvent to propagate the dynamics . while most of the 
 physics - based methods designed to describe a 
 biomolecular system at constant ph use an implicit solvent representation 
 of the solvent , several cphmd methods have been extended to sample , 
 at least conformations , in explicit solvent with both the discrete and continuous protonation models . the methods proposed by baptista et al . and wallace and shen use an implicit 
 solvent potential to sample protonation states while the methods developed 
 by goh et al . and donnini et al . perform -dynamics on the titration a more recent approach by wallace and 
 shen uses a -dynamics approach in pure explicit solvent and 
 adds a counterion whose charge is changed simultaneously with a titratable 
 residue in order to maintain charge neutrality in the unit cell . discrete protonation methods use molecular 
 dynamics to propagate 
 the spatial coordinates , while occasionally interrupting the dynamics 
 to attempt changes to the protonation states of the titratable residues 
 using a metropolis monte carlo criteria . the cphmd method implemented 
 in amber ( and later implemented in charmm ) performs md in gb solvent , periodically attempting 
 to change the protonation state of one or two interacting residues 
 roughly every 10 fs . in the stochastic 
 titration method described by baptista et al . , dynamics is run in 
 explicit solvent for 2 ps , after which 
 a cycle of protonation state change attempts are evaluated using the 
 poisson  boltzmann ( pb ) equation to treat solvation effects 
 for every titratable residue and interacting titratable residue pair . about 40  000 full cycles are attempted each time protonation 
 state changes are attempted . afterward , 
 the solute is held fixed while md is propagated on the solvent to 
 reorganize the solvent distribution to the new set of protonation 
 states . implicit solvent models  in this case average 
 over all solvent degrees of freedom , thereby instantly incorporating 
 the effects of solvent relaxation around discrete protonation state 
 changes . therefore , mc moves in which a protonation state change is 
 attempted have a reasonable probability of succeeding when the solution 
 ph is set close to the pka of the titratable 
 group . when explicit solvent molecules are present , however , the solvent 
 orientation around any solvent - exposed , titratable residue will oppose 
 any protonation state changes . on average , the solvent distribution 
 tends to resist protonation state changes by imposing a barrier on 
 the order of 100 kcal / mol as estimated by measurements in our lab 
 and in others , making titration 
 with discrete protonation states difficult directly in explicit solvent . in this study , we present a new method of performing cphmd simulations 
 in explicit solvent using discrete protonation states . this method 
 is similar in some regards to that used by baptista et al . , and we evaluate its performance on the model 
 compounds , a pentapeptide , and the hewl protein . to enhance the sampling 
 capabilities of this new cphmd method , we use replica exchange in 
 the ph - dimension ( ph - remd ) , whose theory and performance were discussed 
 previously in the context of implicit solvent calculations . this paper is organized as follows : we will first describe the method 
 and its implementation in the theory and methods section , followed by a description of the calculations we performed 
 in the calculation details section . afterward , 
 we will evaluate its performance as well as sensitivity to the method s 
 tunable parameters in the results and discussion section . in this section , 
 we will discuss the details of our proposed method 
 and highlight how it differs from the approach used by baptista et 
 al . the theoretical foundation of our 
 cphmd method is described in detail , as well as the ph - remd method 
 we used in our simulations . in cphmd , structures are sampled from the 
 semigrand canonical ensemble , 
 whose probability distribution function is given by1where  = 1/kbt ,  * is the chemical potential of hydronium 
 ( directly related to the solution ph ) , q is the generalized 
 coordinates of the system particles , p is the conjugate 
 momenta , and n is the total number of titratable 
 protons present in that state . when bold , n refers to 
 the protonation state vector , specifying not only the total number 
 of protons present but on which titratable sites those protons are 
 located . the denominator in eq 1 is the partition function of the semigrand canonical ensemble . to sample from the probability function  in eq 1 , discrete protonation state methods use md with a fixed set 
 of protonation states to sample coordinates and momenta coupled with 
 a mc - based protonation state sampling at fixed conformations throughout 
 the trajectory . baptista et al . showed that standard md , which samples 
 (q , p|n ) , used in conjunction 
 with mc moves on protonation states , which samples from (n|q , p ) , properly samples from the 
 desired probability distribution function (q , p , n ) . in this notation , 
 (q , p|n ) is the conditional 
 probability function of the positions and momenta with fixed protonation 
 states whereas (n|q , p ) is the conditional probability function of the protonation state 
 vector at a fixed protein conformation . in explicit solvent , 
 (n|q , p ) is difficult 
 to sample directly , since the solvent orientation 
 is relaxed with respect to the current protonation state vector . following 
 the arguments of baptista et al . , the system coordinates ( and momenta ) 
 can be separated into solute and solvent degrees of freedom . the protonation state sampling is then performed 
 according to the conditional probability2where qsolvent and psolvent are relaxed solvent distributions 
 of positions and momenta around the protonation state vector , n. the energy differences resulting 
 from the relaxed solvent distributions around the different protonation 
 states are quantities that implicit solvent models strive to reproduce . 
 therefore , the distribution function  in eq 2 can be approximated using continuum models , such 
 as the pb or gb equations , thereby avoiding the otherwise costly solvent 
 relaxation calculation associated with each attempted protonation 
 state change . performing relaxation md on the solvent degrees of freedom should be done after protonation state changes to generate the uncorrelated , 
 relaxed solvent distribution required by eq 2 . manousiouthakis 
 and deem have shown that simply satisfying the weaker balance condition 
 is valid . contrary to the stochastic 
 titration method that calculated solvation 
 free energies using the pb equation to evaluate protonation state 
 changes , we chose to use the gb implicit 
 solvent model for three main reasons . first , amber has numerous gb 
 models readily available , allowing us to use the existing code to evaluate protonation state 
 change attempts . second , results from the original gb - based cphmd 
 implementation by mongan et al . , and from a number of previous studies 
 using the method , have been promising . furthermore , gb was 
 shown to be effective when used in a hybrid solvent method with continuous 
 protonation states , is computationally 
 cheaper than pb , and gb is more easily parallelizable , allowing longer 
 simulations to be performed in the same amount of time . the 
 process of the cphmd method presented here can be divided into three 
 repeating steps , summarized in the workflow diagram in figure 1 . this workflow is very similar to the one presented 
 in ref ( 5 ) ( figure 
 2 ) , although the nature of the mc protonation state move is different . in the proposed method standard md in explicit solvent is carried 
 out using a constant set of protonation states ( an initial set must 
 be provided at the start of the simulation ) . at some point the md is stopped , the solvent ( including any nonstructural ions ) are stripped , 
 the potential is switched to an available gb model , and a set of n protonation state changes are attempted where n is the number of titratable residues . while in principle 
 the md can be stopped randomly with a predetermined probability at 
 any step , in this iteration of our proposed method we run md for a 
 set time interval , md , similar to the stochastic 
 titration method . after the md is 
 halted and the solvent stripped , protonation state 
 changes are proposed for each titratable residue once , in random order , 
 choosing from the available protonation states of that residue excluding 
 the currently occupied state . the electrostatic energy difference 
 between the proposed and current protonation states , as well as the 
 mc decision regarding whether or not to accept the proposed state , 
 are calculated the same way as in the original gb implementation . if the protonation state change is accepted , 
 the  current state is appropriately updated , and the 
 next residue , chosen at random without replacement , is titrated with 
 this new state . for each residue that is titrated , there is 
 a 25% chance that a 
 so - called multisite titration will occur with a neighboring residue ; 
 that is , the proposed change will involve changes to the protonation 
 state of both the residue and its neighbor . . if either residue has more 
 than one titrating proton , the two residues are neighbors if the minimum 
 distance between any pair of titrating hydrogens meets the cutoff . including multisite protonation state jumps is important for systems 
 that have closely interacting titratable residues . without these multisite 
 moves , proton transfers between adjacent titratable residues involved 
 in a hydrogen bond would never occur due to the high penalty of disrupting 
 the interaction by adding another proton or removing the proton involved 
 in the hydrogen bond . this feature was actually present in the initial 
 gb cphmd implementation , and while no mention of it was made in the 
 original paper , a small note was made in the amber users manual . if any of the protonation state change 
 attempts were accepted , 
 the solute is frozen while md is performed on the solvent ( and any 
 ions ) to relax the solvent distribution around the new protonation 
 states . the length of this relaxation is a tunable parameter of the 
 method , which we will call rlx . after the relaxation is complete , the velocities of the solute atoms 
 are restored to their values prior to the relaxation and the standard 
 dynamics is continued . it is worth noting that as the protonation 
 states change during 
 the course of the cphmd simulations , so too does the net charge on 
 the system . because we are using periodic boundary conditions with 
 a lattice - sum method to compute electrostatics , the investigations focusing on charge - dependent finite - size effects 
 have identified artifacts affecting primarily computed free energies 
 and pressure and are larger for smaller unit cells . in the proposed method , the protonation state changes are 
 attempted in implicit solvent with a gb potential , which is entirely 
 unaffected by these finite - size effects . furthermore , since replica 
 exchange simulations in amber require the use of constant volume , 
 no pressure corrections are required , either . workflow of the proposed 
 discrete protonation cphmd method in explicit 
 solvent . following the standard md , the solvent , including all nonstructural 
 ions ( as determined by user - input ) , after that , the solvent 
 and the original settings are restored for the remaining steps . the underlying 
 theory behind replica exchange in ph - space with md run in explicit 
 solvent is unchanged from the version we implemented in implicit solvent 
 earlier . replicas are ordered by their solution ph 
 parameter , and adjacent replicas attempt to exchange their ph periodically 
 throughout the md simulations . the probability of accepting 
 these replica exchange attempts , given by eq 3 , depends only on the difference in the number of titrating protons 
 present in each replica and their respective difference in ph.3where ni is the total 
 number of titratable protons present 
 in state i. because the probability of accepting 
 replica exchange attempts 
 depends only on the number of titratable protons that are present 
 in the system , the number of replicas necessary to obtain efficient 
 mixing in ph - space does not increase as explicit solvent is added . this , coupled with the more accurate sampling found with explicit 
 solvent simulations , makes ph - remd an 
 effective tool for explicit solvent cphmd . to evaluate the performance of the proposed method , 
 we applied 
 it to the amino acid model compounds , a small pentapeptide ( acfca ) , 
 and a protein commonly used in pka calculation 
 studies - the hen egg - white lysozyme ( hewl ) . in order 
 to allow simulations to be run for hundreds of nanoseconds , we implemented 
 mongan et al.s cphmd method with ph - remd on gpus in gb implicit solvent . we used the hewl structure solved in pdb code 1aki(38 ) as our starting structure . the carboxylate residues were 
 initially set in the deprotonated state , and histidine 15 was started 
 in the double - protonated state according to amber defaults . the structure was minimized using 10 steps of the steepest descent 
 algorithm followed by 990 steps of conjugate gradient with 10 kcal / mol . the minimized structure was then heated for 1 ns , varying 
 the target temperature linearly from 10 to 300 k over 667 ps . positional 
 restraints ( 5 kcal / mol ) were placed on the 
 backbone atoms . the heated 
 structure was then equilibrated at 300 k for 2 ns using 
 langevin dynamics with a friction coefficient of 10 ps with 0.1 kcal / mol positional restraints 
 on the backbone . we then ran 500 ns of ph - remd , attempting to change 
 protonation states every 10 fs and attempting replica exchanges every 
 20 fs . we used 14 replicas spanning the ph range 17.5 with 
 a 0.5 ph - unit spacing between adjacent replicas . absolute pka values 
 are very difficult to calculate in 
 solution ; they are impossible using classical force fields . as a result , 
 every physics - based cphmd method uses the idea of a model 
 compound whose experimental pka is easy to measure with a high level of accuracy . an empirical parameter  the reference energy is then added so that cphmd reproduces 
 the experimental pka values of these compounds . the reference energy must be set for the solvation model that is used 
 during the simulations , which was the same gb model 
 that mongan et al . the model compounds have the residue sequence ace - x - nme , where ace is a neutral acetyl capping residue , x is the titratable residue , and nme is 
 a neutral methyl amine capping residue . the available titratable residues in amber are aspartate ( as4 ) , 
 glutamate ( gl4 ) , histidine ( hip ) , lysine ( lys ) , tyrosine ( tyr ) , and 
 cysteine ( cys ) , which are all defined as described by mongan et al . a 10  tip3p solvent buffer was added in a truncated octahedron around each model 
 compound that we simulated . the aspartate model compound was also 
 simulated with larger box sizes15  and 20  buffers  to 
 determine if it had any effect on the calculated pka . after the system topologies were generated , each 
 system was minimized 
 using 100 steps of steepest - descent minimization followed by 900 steps 
 of conjugate gradient minimization . they were then heated at constant 
 pressure , varying the target temperature linearly from 50 to 300 k 
 over 200 ps . the solvated model compounds were then run for 2 ns at 
 constant temperature and pressure . each model compound system 
 was simulated at constant ph and volume 
 for 2 ns , setting rlx = 200fs . each simulation employed ph - remd using six replicas 
 with the solution ph set to pka  
 0.1 , pka  0.2 , and pka  1.2 . to evaluate the effect of the solvent relaxation 
 time , the cysteine model compound was run with rlx set to 10 fs , 40 fs , 100 fs , 200 fs , and 2 ps . a pentapeptide with the sequence ala - cys - phe - cys - ala ( acfca ) was solvated with a 15  
 buffer of tip3p molecules around the solute in a truncated octahedron . the system was minimized using 100 steps of steepest descent minimization 
 followed by 900 steps of conjugate gradient . the minimized structure 
 was heated by varying the target temperature linearly from 50 to 300 
 k over 200 ps at constant pressure . the resulting structure was then 
 simulated at 300 k at constant temperature and pressure to stabilize 
 the system density and equilibrate the solvent distribution around 
 the small peptide . the resulting structure was then used in 
 simulations at six different 
 solution ph values7.1 , 8.1 , 8.3 , 8.7 , 8.9 , and 9.9 . these 
 ph values were chosen because the pka of 
 the cysteine model compound is 8.5 , so the two cysteines of acfca 
 were expected to titrate in this ph range . to demonstrate the effect 
 that ph - remd had on the titration of acfca , two sets of simulations 
 were run cphmd with no exchanges and ph - remd  with each 
 replica being run for 2 ns . we used the 
 structure solved in pdb code 3lzt as the starting structure for our simulations . all eight aspartate residues were renamed as4 , 
 both glutamate residues were renamed gl4 , and histidine 15 was renamed 
 hip in preparation for the ph - remd simulations . all disulfide 
 bonds were added manually in tleap , and the system 
 was solvated with a 10  tip3p water buffer surrounding the protein 
 in a truncated octahedron . we added 26 sodium ions and 11 chloride 
 ions in random locations in the simulation cell to provide an ionic 
 atmosphere and neutralize the system in its initial , default protonation 
 states . the system was minimized using 1000 steps of steepest 
 descent minimization 
 followed by 4000 steps of conjugate gradient , with 10 kcal / mol positional restraints applied to the backbone . the structures 
 were then heated at constant volume , varying the target temperature 
 linearly from 10 to 300 k over 400 ps using langevin dynamics with 
 a 5 ps friction coefficient . the heated structures 
 were then equilibrated for 4 ns at constant temperature and pressure 
 using langevin dynamics with a 2 ps friction coefficient 
 to maintain the temperature and 2 ps coupling 
 constant for the berendsen barostat . the system was then subjected 
 to 112 ns of equilibration md at constant volume and temperature using 
 langevin dynamics with a 2 ps friction coefficient 
 to maintain a constant temperature of 300 k. following the setup 
 stages of the simulations , two sets of ph - remd 
 simulations were carried out for 122 and 150 ns with 12 replicas spanning 
 the ph range from 3 to 8 at 1 ph - unit intervals to characterize 
 the acidic - range titration behavior of hewl . all systems were 
 parametrized using the amber ff10 force field , which is equivalent 
 to the amber ff99sb force field for proteins . the tleap program of the ambertools 12 program 
 suite was used to build the model compound and acfca molecules , to 
 add hydrogen atoms hewl and to solvate each system that was run in 
 explicit solvent . all simulations were performed using either 
 the sander or pmemd module of a 
 development version of amber 12 . langevin 
 dynamics was used in every simulation to maintain constant temperature 
 with collision frequencies varying from 1 ps to 
 5 ps , and the random seed was set from the computer 
 clock to avoid synchronization artifacts . the berendsen barostat was used to maintain constant pressure for 
 the equilibration dynamics with a coupling constant of 1 ps . all molecular dynamics , including the solvent relaxation 
 dynamics , 
 are run with a 2 fs time step , constraining bonds containing hydrogen 
 using shake . replica exchange attempts between 
 adjacent replicas were made every 200 fs for all ph - remd simulations . long - range electrostatic interactions were treated 
 with the particle - mesh 
 ewald method using a direct - space and van 
 der waals cutoff of 8  . defaults were used for the remaining 
 ewald parameters . the gb model proposed by onufriev et al . , indicated 
 by the parameter igb = 2 in sander , was used to evaluate 
 the protonation state change attempts to be consistent with the original 
 implementation in implicit solvent . the 
 intrinsic solvent radius of the carboxylate oxygen atoms was reduced 
 from the standard value of 1.5  to 1.3  to compensate 
 for the effect of having 2 dummy protons present on each oxygen in 
 the syn- and anti - positions in the explicit solvent simulations . the deprotonation fraction ( fd ) and 
 ph for each simulation  and each window of the running averages was 
 fitted to the hill equation ( eq 4 ) using the 
 levenberg marquardt nonlinear least - squares algorithm implemented 
 in scipy to compute the pka and hill coefficient ( n ) . all pka values reported for titratable residues in this paper 
 correspond to the value computed by fitting fd from the simulations at every ph to eq 4 over the specified time interval.4 first , we will analyze the long 
 simulations in implicit solvent 
 that demonstrate the weaknesses of implicit solvent models . next , 
 we will analyze our proposed cphmd and ph - remd methods as well as 
 ways to optimize its overall performance . we will start by discussing 
 the behavior of the model compounds when the size of the unit cell 
 and the length of the relaxation dynamics ( rlx ) is varied . we will follow this discussion with a similar 
 analysis on a slightly larger system  acfca  before discussing 
 the application of our proposed method to hewl . long simulations 
 are required to observe rare events or phenomena that occur on long 
 time scales , but often expose deficiencies in a computational model 
 that are frequently missed in shorter simulations . therefore , we ran cphmd gb simulations of hewl for 600 
 ns  more than an order of magnitude longer than what has previously 
 been published  to determine if these calculations were stable 
 and could be used to study the dynamical behavior of proteins over 
 long time scales . the backbone rmsd with respect to the starting 
 crystal structure , shown in figure 2 for the 
 replicas at ph 1 , 4 , and 7 , reaches as much as 10  at lower 
 ph values , suggesting that the native state is unstable in gb implicit 
 solvent at long time scales . as the lysozyme is active at the lysosomal 
 ph around 4.5 , cphmd simulations near this ph should reflect this 
 stability . as the backbone rmsd from the crystal structure increases , 
 the predicted pka value worsens . we computed 
 the root mean squared error ( rmse ) of the computed acidic residue 
 pka values compared to experiment . the 
 predicted pka was computed for each residue 
 by taking a running average with a 10 ns window of the deprotonation 
 fraction for each residue at each ph in the ph - remd simulation . the 
 rmse of the calculated pka values from 
 experiment averaged over all titratable residues is shown in figure 3 . rmsd compared to the 1aki crystal structure for the ensembles 
 at 
 ph 1 , 4 , and 7 . the time series is shown on the left with the normalized 
 histograms shown on the right . rmse of all acidic titratable residue pka values compared to experiment during 
 the course of the simulation . a 10 ns window was used for the running 
 average of the computed pka . to study the effect 
 that the unit cell size has on titrations in our proposed method , 
 we prepared three simulations of the aspartate model compound with 
 different tip3p solvent buffers surrounding it . we prepared systems 
 with a 10  , 15  , and 20  tip3p solvent buffer around 
 the model aspartate . because protonation state sampling takes 
 place in gb solvent without periodic boundary conditions , any effect 
 of the box size on calculated pka values 
 can only arise due to alterations of the structural ensembles induced 
 by artifacts from the box size . the calculated pka values of the three systems were 4.02  0.07 , 4.05 
  0.08 , and 4.12  0.07 for the 10  , 15  , and 
 20  solvent buffer systems , respectively . to estimate the uncertainties 
 we divided each simulation into 100 ps chunks and took the standard 
 deviation of the set of 20 pka values 
 calculated from those segments . to further demonstrate the insensitivity 
 of box size to ph - remd 
 titrations , we plotted the solvent radial distribution functions ( rdfs ) 
 around the center of mass of the carboxylate functional group in figure 4 . the insensitivity of the pka and solvent structure with respect to the model compound 
 provides strong evidence that no undue care is necessary when choosing 
 the size of the solvent buffer for these types of simulations . radial distribution 
 functions ( rdfs ) of solvent oxygen atoms ( o ) 
 and hydrogen atoms ( h ) with different unit cell sizes . the shown measurements10 , 
 15 , and 20   represent the size of the solvent buffer 
 surrounding the solute . rdf plots for three different phs are shown , 
 highlighting the ph dependence of the solvent structure around the 
 carboxylate of the aspartate model compound and its invariance to 
 box size . an important approximation in 
 the proposed method is that the protonation state sampling  
 from eq 2 can be replaced using an implicit 
 solvent model followed by relaxation md to generate the relaxed solvent 
 positions and momenta . the question then becomes how long this relaxation 
 dynamics should be run . to address this , 2 ns of constant protonation 
 molecular dynamics simulations were run on the model cysteine compound 
 in both protonation states  protonated and deprotonated after 
 the same minimization and heating protocols were used as for the other 
 model compound simulations . the protonation state was then swapped 
 for the final structures of both simulations , and md was performed 
 while constraining the solute position for 20 ns , equivalent to the 
 relaxation dynamics protocol in our proposed method . the optimum 
 value for rlx is 
 the time after which the energy of the relaxation trajectory stabilizes 
 and the simulation loses all memory of its initial configuration . to be truly equivalent to having been chosen at random , the final , 
 relaxed solvent distribution must be completely uncorrelated from 
 the initial distribution at the time the protonation state was changed . to probe the necessary time scales for these relaxation dynamics , 
 the energy of each snapshot in the relaxation trajectory is plotted 
 alongside the autocorrelation function of that energy in figure 5 to clearly demonstrate the first , the model compounds are fully solvent - exposed due to their 
 small size , which results in a worst - case scenario in terms of the 
 number of water molecules that must be reorganized during the relaxation 
 dynamics . the optimum rlx value 
 for model compounds is expected to be an upper - bound on the values 
 required for larger systems . second , cysteine is the smallest and 
 simplest of the titratable amino acids , eliminating potential complications 
 from tautomeric states compared to aspartate , glutamate , and histidine . the relaxation energies plotted in figure 5 begin to stabilize after 4 to 6 ps of relaxation dynamics and are 
 largely uncorrelated within that same time scale . however , because 
 4 ps of md  corresponding to 2000 steps of dynamics with a 2 
 fs time step  adds dramatically to the cost of cphmd simulations 
 in explicit solvent , we will analyze the approximation of using a 
 significantly smaller value for rlx . both the relaxation energies and autocorrelations drop very 
 sharply 
 at the start of the relaxation dynamics , so the majority of the benefit 
 gained by relaxing the solvent is realized within the first few steps . for both simulations , the first 200 fs of relaxation dynamics resulted 
 in 70% of the total relaxation energy in calculations . the correlation 
 of the relaxation energy decays similarly , so the assumption that 
 the relaxed solvent distribution is uncorrelated from its starting 
 point is a reasonable approximation . the relaxation of the protonated state , 
 starting from the protonated 
 trajectory , is shown in blue with its autocorrelation function shown 
 in purple . the relaxation of the deprotonated state from an equilibrated 
 snapshot from the protonated ensemble is shown in red with its autocorrelation 
 function shown in green . here , pr and dr stand for protonated - relaxation and deprotonated - relaxation , respectively . to validate the use of a shorter 
 rlx , we titrated the aspartate 
 model compound using ph - remd with 
 five different values for rlx10 
 fs , 40 fs , 100 fs , 200 fs , and 2 ps . the calculated pka for each simulation was 4.10 , 4.08 , 4.07 , 4.10 , and 
 4.05 , respectively , for the listed relaxation times . furthermore , 
 comparing the solvent radial distribution functions of the different 
 solvent relaxation times ( figure 6 ) shows little 
 dependence of the solvent distribution on the value of rlx . rdfs of water oxygen atoms ( o ) and hydrogen 
 atoms ( h ) around the 
 center - of - mass of the carboxylate group of the model aspartate molecule 
 at different solution phs . the small 
 peptide chain acfca , described in section 3.3 , was chosen as a test due to its small size and predictable titration 
 behavior . the simplicity of the system makes it an ideal test ; its 
 small size mitigates the conformational sampling problem , and the 
 simple titrating behavior of cysteine further simplifies protonation 
 state sampling . unlike aspartate and glutamate , which have the four 
 defined tautomeric states defined by anti- and syn - protonation on 
 each of two carboxylate oxygens , and histidine which has two tautomeric 
 states on the imidazole , cysteine has only one protonated and one 
 deprotonated state , presenting fewer degrees of freedom that must 
 be exhaustively sampled . each cysteine is in a slightly different 
 microenvironment due to the different charges of the n- and c - termini . 
 due to its proximity to the n - terminus , cys-2 is expected to experience 
 a negative pka shift with respect to the 
 model compound due to the electrostatic interactions with the positively 
 charged terminus . cys-4 , on the other hand , is expected to experience 
 a pka shift in the opposite direction 
 due to the electrostatic interactions with the negatively charged 
 c - terminus . we ran simulations at ph 7.1 , 8.1 , 8.3 , 8.7 , 8.9 , 
 and 9.9 to sufficiently 
 characterize the titration behavior of both cysteine residues around 
 their pka values . one set of replicas 
 was run with ph - remd while the other set was run using cphmd ( i.e. , 
 without attempting exchanges between the replicas ) . the titration 
 curves for both sets of simulations , shown in figure 7 , demonstrate the importance of using ph - remd in constant 
 ph simulations in explicit solvent . as expected , the ph - remd simulations 
 revealed pka shifts of 0.2 pk units for cys 2 and + 0.9 pk units for 
 cys 4 with respect to the model cys compound . even for a simple 
 system such as acfca , using ph - remd on top of 
 standard cphmd simulations results in a drastic improvement in titration 
 curve fit  a result of improved protonation state sampling . the residual sum of squares ( rss ) , a quantity that measures how well 
 an equation fits a data set , shows drastic improvement using ph - remd . the rss for cys-2 and cys-4 using cphmd was 9  10 and 7  10 , respectively . for the ph - remd 
 simulations , on the other hand , the rss was reduced by several orders 
 of magnitude to 7  10 and 9  10 for cys-2 and cys-4 , respectively . . results from 
 cphmd ( no replica exchange attempts ) 
 and ph - remd are shown in the plots on the left and right , respectively . hewl is 
 a common benchmark for pka calculations 
 because it has been studied extensively both experimentally and theoretically , and it has a large number of titratable residues  some 
 with a marked pka shift compared to the 
 isolated model compound . the simulations run in explicit solvent 
 revealed far more stable trajectories over all ph values than their 
 analogs run in implicit solvent over the 150 ns time scale of the 
 explicit solvent simulations . in figure 8 , 
 the plots of backbone rmsds are bounded below 3.5  for the duration 
 of the ca . the predicted pka values for the titratable 
 residues , summarized in table 1 for both sets 
 of simulations , show good agreement to experiment . the agreement is 
 significantly better than the predictions from the implicit solvent 
 calculations over a similar time scale . with the exception of aspartate 
 119 ( and aspartate 52 in the 150 ns simulation ) , the predicted pka values of all residues were within 1 pk unit from the experimental values given in ref ( 49 ) . backbone rmsd distributions 
 for hewl simulations in explicit solvent 
 with solution ph set to 1 , 4 , and 7 . time series are shown on the 
 left and histograms are shown on the right . furthermore , the large fluctuations in the rmse throughout 
 the 
 course of the simulation  between 0.70 and 1.25 seen in figure 9suggest that differences in reported rmse 
 around 0.1 pk units are statistically insignificant 
 for short cphmd and ph - remd simulations on the order of 10 ns . the 
 standard deviation of the pka rmse plotted 
 in figure 9 is 0.12 pk units 
 with a correlation time of roughly 25 ns . the standard error of the 
 mean , given by ( /n ) where  
 is the variance and n is the number of independent 
 samples , is 0.05 pk units over 150 ns . therefore , 
 there is no statistically significant difference between methods with 
 reported pka rmses within 0.10 pk units of each other ( 0.05 pk units 
 for each simulation ) over 150 ns . ( results are shown for the 
 final 500 ns of the 600 ns implicit solvent simulation ) . unlike the implicit solvent calculations , the quality of the predicted 
 pka values in the proposed method does 
 not appear to degrade over time in long simulations . the deviation 
 of each residue compared to experiment and the total rmse is shown 
 in figure 9 . the pka at each time step is computed from a running average of deprotonation 
 fraction for replicas at each ph using a window size of 5 ns and fitting 
 to eq 4 the same way as described for figure 3 . the smaller backbone rmsd and drastically 
 improved pka predictions are strong evidence 
 that the current method 
 improves conformational sampling by employing an explicit representation 
 of the solvent . because the protonation state sampling is performed 
 with the same gb potential in both the original and proposed methods , 
 any differences in the predicted pka values 
 must be caused by differences in conformational sampling . given the 
 rigorous nature of the experimental measurements , this provides strong evidence that the proposed method 
 improves significantly upon the original by improving conformational 
 sampling with an explicit solvent representation . ( left ) plot showing the 
 rmse of all acidic titratable residue pka values compared to experiment during 
 the course of the simulation . a 5 ns window was 
 used for the running average of the computed pka . ( right ) plot showing the autocorrelation function of the 
 rmse time series , showing that the correlation time between statistically 
 independent samples is roughly 25 ns . we have 
 extended the constant ph molecular dynamics method developed 
 by mongan et al . so that the dynamics 
 can be run in explicit solvent . we tested a wide range of parameters 
 in our proposed method for their effect on the conformational and 
 protonation state sampling of small test systems . because these test 
 systems are small and their titratable sites are completely solvent - exposed , 
 they likely represent the highest level of sensitivity to these various 
 parameters . in particular , we found that the box size of the 
 unit cell had 
 no discernible effect on the titration behavior of the aspartate model 
 compound , given cell sizes that ranged from 20  in diameter one 
 of the smallest sizes permissible when using the minimum image convention 
 with an 8  cutoff  to 40  in diameter . another 
 key aspect of the current method is the necessity to relax 
 the solvent around any new protonation state selected by the mc moves 
 carried out in gb . by analyzing the decay of the potential energy 
 in the solvent relaxation dynamics , we determined that 4 ps of md 
 was sufficient to stabilize the energy of the solvent distributions 
 and generate relaxed solvent conformations whose energies are uncorrelated 
 from the initial arrangements . however , given the expense of such 
 a long relaxation period , we investigated using fewer relaxation steps 
 to increase the simulation efficiency and found shorter times  down 
 to 0.2 ps  had no measurable effect on the calculated pka and very little effect on the solvent distribution 
 around the model cysteine compound . further tests on a small 
 pentapeptide test system with two titratable 
 sites ( acfca ) showed the importance of using ph - remd over conventional 
 cphmd with the proposed method . while we showed that the enhanced 
 protonation state sampling of ph - remd results in smoother titration 
 curves for complex proteins in implicit solvent , even the simplest systems in explicit solvent require ph - remd 
 to obtain a smooth titration curve . we tested the proposed method 
 on hewl , a very common pka benchmark system . we found that the proposed method 
 of using gb implicit solvent to evaluate protonation state changes 
 and explicit solvent to propagate dynamics yielded stable trajectories 
 whose predicted pka values agreed well 
 with experiment . our results show that using the proposed method leads 
 to a significant improvement in how systems are modeled at constant 
 ph compared to the original method that used gb to propagate system 
 dynamics . often , the most interesting titratable residues in 
 biological systems 
 have a large pka shift compared to the 
 model compound . these highly perturbed residues have environments 
 drastically different than the one provided by bulk solvent , and the 
 conformational sampling must be both accurate and extensive to yield 
 accurate pka predictions . in future work we will explore the use of enhanced sampling techniques in conjunction 
 with ph - remd in an attempt to improve the efficiency of the conformational 
 sampling in explicit solvent , such as accelerated md and temperature - based 
 remd .
by 
 utilizing graphics processing units , we show that constant ph 
 molecular dynamics simulations ( cphmd ) run in generalized born ( gb ) 
 implicit solvent for long time scales can yield poor pka predictions as a result of sampling unrealistic conformations . 
 to address this shortcoming , 
 we present a method for performing constant 
 ph molecular dynamics simulations ( cphmd ) in explicit solvent using 
 a 
 discrete protonation state model . 
 the method involves standard molecular 
 dynamics ( md ) being propagated in explicit solvent followed by protonation 
 state changes being attempted in gb implicit solvent at fixed intervals . 
 replica exchange along the ph - dimension ( ph - remd ) helps to obtain 
 acceptable titration behavior with the proposed method . 
 we analyzed 
 the effects of various parameters and settings on the titration behavior 
 of cphmd and ph - remd in explicit solvent , including the size of the 
 simulation unit cell and the length of the relaxation dynamics following 
 protonation state changes . 
 we tested the method with the amino acid 
 model compounds , a small pentapeptide with two titratable sites , and 
 hen egg white lysozyme ( hewl ) . 
 the proposed method yields superior 
 predicted pka values for hewl over hundreds 
 of nanoseconds of simulation relative to corresponding predicted values 
 from simulations run in implicit solvent .