Hello,
I am actually implementing the changes in titratable_residues.py only now,
and I think I am not 100% sure about what I should do.
Jason, you said that :
*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.*
So in my cpin file, the Reference energy is STATENE right? So in my case I
have STATENE=0.000000,21.248624,21.248624,21.248624,21.248624,
So it means that the reference energy I should use is 21.248624?
If we do *Ref Energy + R*T*ln(10)*pKa =*
21.248624+0.00199*300*log(10)*4.82= 27.8744
so the TI energy should have been 27.8744 to give me the correct pKa, is
that correct?
Here is an extract of the glutamate entry in titrable_residues.py:
# Glutamate
refene1 = _ReferenceEnergy(igb1=0, igb2=0, igb5=0, igb7=0, igb8=0)
refene1.solvent_energies()
refene1.dielc2_energies(igb1=0, igb2=0, igb5=0, igb7=0, igb8=0)
refene1.dielc2.solvent_energies()
refene2 = _ReferenceEnergy(igb1=3.89691326, igb2=8.4057785, igb5=8.0855764,
igb7=5.305949, igb8=8.3591335)
refene2.solvent_energies(igb2=15.20019319)
refene2.dielc2_energies(igb2=3.455596, igb5=3.957270)
refene2.dielc2.solvent_energies()
refene2.set_pKa(4.4, deprotonated=False)
So in order to "correct" the TI energy in titrable_residues.py, I simply
need to replace refene2.solvent_energies(igb2=15.20019319) with
refene2.solvent_energies(igb2=27.8744) and that should give me the correct
pKa when I re-run the simulation. Is that it? I feel like I am missing
something... Especially because it is a quite large difference.
If I want to do the same for Lys were I got a pKa of 10.38 and
STATENE=-0.845507,0.000000,
I should do TI = -0.845507 + 0.00199*300*log(10)*10.38= 13.4233 and then I
would need to replace refene1.solvent_energies(igb2=-15.1417977) with
refene1.solvent_energies(igb2=13.4233)? That doesn't look correct.
Shouldn't it be actually a subtraction in both cases, i.e. *Ref Energy -
R*T*ln(10)*pKa*, in which case it would give an energy of 14.6228
and -15.1143 for Glu and Lys respectively, which would be much more in line
with the original refene2.solvent_energies of15.20019319 and -15.1417977
respectively?
Sorry I am a bit lost with this equation and the different values, I would
really appreciate if you could clarify this a bit.
Many thanks in advance,
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 Tue Mar 14 2017 - 08:00:03 PDT
