Hello,
I was wondering, if I want to correct the reference energy for Histidine
residue can I use the same formula TI = Ref Energy + R*T*ln(10)*pKa ? In
that case how do I deal with the two sites that have a different pKa?
Can I just run the pH-REMD (at pH = pKa ± 0.1 pKa ± 0.2 pKa ± 1.2) using
pKa=6.5 and then do the same for pKa=7.1 and then correct TI for each pKa?
It sounds too simple to do that in the case of His and I fell like there is
more hidden complexity because of the two sites?
Many thanks in advance for your help,
Eric
On 10 February 2017 at 10:34, Eric Lang <eric.lang.bristol.ac.uk> wrote:
> 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
>
--
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 Mar 03 2017 - 07:00:02 PST
