Hi Jason,
Thanks a lot for this precisions. I definitely explains a lot.
Many thanks also for the information on how to refine the reference energy,
I though it would be much more complicated with the need to run TI. I am
going to try this.
Many thanks!
Eric
On 9 February 2017 at 03:19, Jason Swails <jason.swails.gmail.com> wrote:
> Hi Eric,
>
> I don't check the Amber list nearly as often as I used to and only just
> saw this question.
>
> On Tue, Feb 7, 2017 at 8:22 AM, Eric Lang <eric.lang.bristol.ac.uk> 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.
>>
>
> ​The first thing I would have suggested is precisely this :). In explicit
> solvent especially, pH-REMD is absolutely critical to breaking out of local
> minima and sampling phase space. This is discussed to some extent in the
> original paper. Basically, the only way to get the ion to "let go" of the
> titratable site is to let the configuration diffuse to a lower pH to force
> protonation and then diffuse back up after the Na+ diffuses away (or vice
> versa for Cl-). Another way of thinking about this is that large barriers
> in configuration space sampling can be circumvented by traversing state
> space (the thermodynamic state here being solution pH) in which the
> barriers are much lower.
>
> So this is consistent with what I've observed.
> ​
>
>> 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.
>>
>
> ​As Vinicius said, the only affect that the explicit ions could *possibly*
> have on the calculated pKa is by modifying the configuration ensemble.
> When you run short simulations that get stuck in a local minimum, these
> changes can be pronounced for the model compounds (which is why you saw
> large differences in your first simulations). When you sample more
> completely, these differences largely disappear.
>
> 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.
>>
>
> A very good level of accuracy for CpHMD is to be within 0.8-0.9 pK units
> (although for the model compounds it should be better). The glutamate
> result is a bit worrisome, but no exceptionally so. The reference energy
> may need some level of correcting (the second set of simulations you ran
> are more thorough than the set I used to tweak the reference energies for
> explicit solvent), but the correction should not be that large.
>
> Histidine is another story. The 6.5 value stored in cpinutil.py was one
> of the two sites. The other site has a pKa of 7.1. So the pKa should be
> something slightly less than 6.5 (which it is), since when the site with a
> higher pKa is protonated, the other site can still be deprotonated with a
> small probability.
>
> 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.
>>
>
> ​This is good to know.
> ​
>
>> 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?
>>
>
> ​Maybe.
> ​​
>
>> 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?
>>
>
> ​The reference energy may benefit from some adjustment.​
>
> 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.
>>
>
> ​It sounds like fine-tuning the parameters would be helpful. It can
> actually be done in a single step. Knowing that the total reference energy
> is
>
> Ref Energy = TI - R*T*ln(10)*pKa
>
> If we want to solve for the "correct" TI energy, you can isolate that by
> calculating
>
> TI = Ref Energy + R*T*ln(10)*pKa
>
> So in your case, take the Reference energy from the cpin file, add
> 0.00199*300*log(10)*4.82, and that value is what the TI energy *should*
> have been. You can then insert that value back into titratable_residues.py
> where refene2.solvent_energies() is set for the GL4 residue.
>
> Recompile, regenerate your cpin file, and rerun your calculations and you
> should see that the pKa comes out more or less perfect.
>
> Hope this helped,
> Jason
>
> --
> Jason M. Swails
>
--
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
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Fri Feb 10 2017 - 03:00:02 PST
