Re: [AMBER] CpHMD – Effect of explicit ions when calculating the pKa of model compounds

From: Adrian Roitberg <>
Date: Tue, 7 Feb 2017 08:28:11 -0500

Dear Eric

I am sorry for the delay answering your query.

My student, Vinicius Cruzeiro, gave a draft of a response to you and I
sat on it too long. I copied it below, and point that that it was all
his work, while I agree with his conclusions.



Dear Eric Lang,

Thanks for bringing this interesting question up. In order to better
address your question, I have ran the explicit CpHMD simulations again
but longer,for 50 ns. In my simulations I place the ACE-CYS-NME system
inside a 1680 water molecules box, of volume 59052.89 A^3. I have done
the simulations randomly adding no ions, 1 Na+ Cl- pair, or 4 Na+ Cl-
pairs. From the explicit solvent point of view, this correspond to the
following ionic strengths: 0 M,0.028 M, and 0.112 M. I’ve also ran the
simulations for two different GB salt concentrations: 0.1 M or 0.4 M.
These are the results:

GB Salt Concentration


No ions


1 Na+ Cl- pair


4 Na+ Cl- pairs

0.1 M










0.4 M










There are two points I would like to discuss here:

1) Although the MD is done with Explicit Solvent, the MC steps are done
based on GB calculations. Although the explicit ions are not directly
included in the GB calculations, they might affect the conformational
distribution of the ACE-CYS-NME system, mainly if there a closer
interaction. What the results show is that this effect is minimum on the
pKa., as you can see the results are basically the same by adding
explicit ions or not.

2) We see the GB Salt Concentration change from 0.1 to 0.4 M affects the
pKa a little bit. This effect is expected by the construction of GB (you
would be able to predict this behavior from the GB equations).
Experimentally speaking, there is also a generally small ionic strength
dependence of the pKa, see:

The reference Free Energy necessary for AMBER’s CpHMD methodology needs
to be computed reproducing as close as possible the environment in which
the experimental pKa reference value was obtained. Once this value is
obtained, it should be good to be used to make predictions even on other
GB salt concentrations. The table above shows that, for the system
studied, adding ions explicitly makes no significant difference in the
pKa results (thus no significant difference too on the reference Free
Energy calculation). I believe this conclusion most hold true for other
systems unless there is a fix and close interaction between the explicit
ion ions and the system.

Best regards,


On 2/7/17 8:22 AM, Eric Lang wrote:
> Hello,
> This post is partly a reply to my previous post but also extends on it.
> I decided to test more carefully the effect of having explicit Na+ and Cl-
> ions at a concentration of 0.1 M with the model compounds for Lys, His, Glu
> and Cys (shamelessly the titrable residues I am currently interested in).
> To do this, I prepared the systems ACE-X-NNE (where X is the titrable
> residue) and solvated them in a 20 A TIP3P solvent buffer and added 7 or 6
> explicit Na+ and Cl- ions (depending on the number of water molecules added
> to solvate the residues of interest). After minimisation, heating and
> equilibration, I ran a 10 ns CpHMD simulation at pH=pKa and then pH-REMD at
> pH = pKa ± 0.1 pKa ± 0.2 pKa ± 1.2 (as in the Jason’s paper on explicit
> solvent pH-REMD) for 20 ns. Results between CpHMD and pH-REMD are
> consistent and here I will just show the pH-REMD results.
> Basically, I do not see significant differences in the calculated pKa
> whether I add explicit ions or not (_salt here means in the presence of 0.1
> M of explicit NaCl):
> Lys Lys_salt Cys Cys_salt Glu Glu_salt His
> His_salt
> pKa 10.35 10.38 8.35 8.34 4.82 4.82 6.31 6.31
> n 1.01 1.00 1.02 1.01 0.96 0.97 1.01
> 1.01
> I do find it a bit surprising as, after all, I would have expected that the
> presence of ions interacting with the titrable residues would make some
> differences as my early results suggested, nonetheless these much more
> thorough simulations suggest the effect is negligible.
> Please note that I did not include error bars by either dividing the
> pH-REMD trajectories into small chunks or by running multiple simulations.
> However, the good agreement between the predicted pKa values obtained from
> pH-REMD and from CpHMD (the values agree at a maximum of ± 0.03 pK units)
> suggests that dividing the pH-REMD into chunks would lead to similar
> variations.
> However, I am now a bit concerned by the shift between calculated and
> theoretical pKa, especially in the case of Glu (4.8 calculated instead of
> the 4.4 it has been parametrised for), and to a lesser extend His (6.3 vs.
> 6.5). Are those kinds of deviations within the accuracy of the method? I
> would have expected deviations of no more than ± 0.1 pK units.
> I thought this could be due to the force field, as I am using ff14SB (see
> this post for more details:,
> so I repeated the simulations with the original ff10 force field
> Lys_ff10 Cys_ff10 Glu_ff10 His_ff10
> pKa 10.39 8.52 4.90 6.39
> n 1.01 1.01 0.97 1.02
> In the case of Lys, Glu and His, the results are similar between the two
> force fileds, suggesting that ff14SB can be used for constant pH MD.
> However, the pKa of Cys is slightly higher with ff10 than with ff14SB,
> could it be somehow because the unprotonated Cys is not defined in the
> standard ff14SB?
> Nonetheless the pKa of Glu is still quite high compared with the
> theoretical one. Could it be a problem with the parametrisation of GL4 in
> amber 16?
> I would really appreciate if you could share your opinion on either or not
> these differences should be considered problematics and if some fine tuning
> of the parameters / reparametrisation would be useful.
> I am happy to share my input/output files if needed.
> Many thanks in advance,
> Eric
> On 30 January 2017 at 12:33, Eric Lang <> wrote:
>> Hello,
>> I have noticed in the original paper on CpHMD in explicit solvent (
>> that the pKa of the model
>> compounds were calculated with the model compounds solvated in explicit
>> water but in the absence of explicit ions. However, when switching to the
>> GB potential in the explicit solvent CpHMD workflow, an ionic concentration
>> of 0.1 M is used. Therefore, the reference energies are calculated assuming
>> a 0.1 M ionic concentration but the dynamics of the model compounds are
>> performed without ions.
>> In contrast, in the hen egg white lysozyme example of the paper, Na+ and
>> Cl- ions are added to neutralise the system and have a 0.1 M concentration
>> of NaCl. This is also done in Jason’s tutorial on Explicit solvent MD.
>> I was therefore wondering if there was a reason for not adding ions to the
>> water box when calculating the pKa of the model compounds as it would make
>> sense to do it to remain consistent with the GB calculations.
>> I ran a quick 1 ns CpHMD test of the Cys model compound (ACE-CYS-NME)
>> following the same protocol as reported in the paper.
>> If I do not add any explicit ions, as per the paper, I end up after at pH
>> =8.5 with the expected pKa/protonated fraction:
>> Solvent pH is 8.500
>> CYS 2 : Offset 0.007 Pred 8.507 Frac Prot 0.504 Transitions 3403
>> Average total molecular protonation: 0.504
>> However, if I add 1 Na+ and 1 Cl- to approximate a 0.1 M NaCl
>> concentration and run the same CpHMD, then I obtain the following:
>> Solvent pH is 8.500
>> CYS 2 : Offset -0.187 Pred 8.313 Frac Prot 0.394 Transitions 3923
>> Average total molecular protonation: 0.394
>> If I increase the simulation time to 3 ns this result is confirmed:
>> Solvent pH is 8.500
>> CYS 2 : Offset -0.171 Pred 8.329 Frac Prot 0.403 Transitions 11383
>> Average total molecular protonation: 0.403
>> When I look at the trajectory, the Na+ ion interacts with the S- of Cys,
>> so it is not due (or only due) to a change in dynamics due to the presence
>> of the ions but because of electrostatic interactions between the
>> deprotonated Cys state and the sodium ion.
>> I understand that this is a simple test and more work would be required to
>> indeed acuratly capture the effect of adding explicit ions to the model
>> compounds simulations, having for example a larger solvent box would enable
>> to better approximate the 0.1 M NaCl concentration, and running pH-replica
>> exchange would give more accurate results. In addition, I didn’t try any
>> other titrable residues at this stage either.
>> Nonetheless, I am a bit worried about these results and the effect of not
>> having explicit ions when the reference energies where calculated. Overall
>> is this not going to lead to a small error in the calculated pKa due to the
>> fact that the model compounds were simulated in the absence of explicit
>> ions? When running the simulations of a protein of interest with explicit
>> ions, the same ions-charged residues interaction will occur. But perhaps
>> I am missing something important there, so I would be grateful to have your
>> opinion on this matter.
>> Many thanks in advance,
>> Eric
>> --
>> Eric Lang
>> BrisSynBio Postdoctoral Research Associate Modelling
>> Centre for Computational Chemistry
>> School of Chemistry - University of Bristol
>> Bristol BS8 1TS - United Kingdom

Dr. Adrian E. Roitberg
University of Florida Research Foundation Professor.
Department of Chemistry
University of Florida
AMBER mailing list
Received on Tue Feb 07 2017 - 05:30:09 PST
Custom Search