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

From: Eric Lang <>
Date: Tue, 7 Feb 2017 13:22:16 +0000


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
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

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

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,


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

Eric Lang
BrisSynBio Postdoctoral Research Associate Modelling
Centre for Computational Chemistry
School of Chemistry - University of Bristol
Bristol BS8 1TS - United Kingdom
AMBER mailing list
Received on Tue Feb 07 2017 - 05:30:07 PST
Custom Search