[gmx-users] Re: Glutamate to Alanine Mutation
Sai Kumar Ramadugu
sramadugu at gmail.com
Mon Apr 1 04:52:16 CEST 2013
Dear Prof Mobley,
I have some additional questions in understanding the free energy
calculations.
1) When doing a mutation from a larger aminoacid to a smaller one, one
necessarily ends up with dummy atoms and atoms that have different mass. In
my case when the glutamate is converted into a glycine or an alanine a
carbon will become a hydrogen (change in mass) and the rest of the excess
atoms will become dummies with the mass they had in Glu.
The change in mass has an effect in the kinetic energy part of the
Hamiltonian contribution to the free energy and this change is properly
computed by Gromacs as <dKe/dL>. However it is clear from this procedure
that because the dummies have mass, the mass of Ala or Gly will not be the
correct and therefore <dKe/dL> would be quite meaningless. I read papers
that suggest only integrating the potential energy component <dV/dL>. In
fact for example, the Amber11 manual indicates that Amber mutations only
change the potential but keep all masses the same as they were before the
alchemical change.
Clearly having dummies with no mass would fix this issue but it would be
problematic because F=m*a and with m=0 for a given finite F, the
acceleration on the dummies would be infinite.
In Monte Carlo all of this is irrelevant because one only computes the
potential contribution to the free energy difference but in MD the story is
different. What is the customary approach to deal with this? Just use
<dV/dL> and discard <dKe/dL>?
2) The g_bar analysis appears to assume that the interpolation between
states a and b is linear (g_bar issues a warning "using the derivative data
(dH/dlambda) to extrapolate delta H values. This will only work if the
Hamiltonian is linear in lambda." ). However when using softcore potentials
lambda also enters in the expression for "r". It would appear that the
interpolation is not really linear even though the LJ transformation also
looks like
Va *lambda+ (1- lambda)* Vb
Nonetheless I see g_bar being used for such transformations. Could you shed
light on this?
3) Also it is common to use small lambda spacings (close to lambda=0 and
lambda=1) but somewhat larger spacing between 0.1 and 0.9. How does g_bar
handle this? Is the integration and error analysis done correctly even when
the lambda spacing may be different at different lambda values?
Thank you for your time.
Regards
Sai
On Tue, Dec 4, 2012 at 5:27 PM, Sai Kumar Ramadugu <sramadugu at gmail.com>wrote:
> Dear Prof Mobley,
>
> Thank you very much for your time.
>
> While I was waiting for your opinion on my topology files, I did several
> tests to see what happens.
> The tests I did are:
>
> 1. Mutate only the amino acid (Glutamate -1 charge) and no mutation of K+
> ion, but have positional restraints on the ion.
> 2. Mutate only the amino acid (Glutamate -1 charge) and no mutation of K+
> ion and no positional restraints or constraints on the system
> 3. The actual simulation that caused me trouble, ie with position
> restraints, constraints and mutating the ion concurrently.
>
> Finally I decided to use the simulation data from no restraints or
> constraints and no mutation of K+ ion.
> With this approach now I am getting a positive deldelG value for the
> mutation going from Glutamate to Alanine.
>
> On Mon, Dec 3, 2012 at 4:07 PM, David Mobley <dmobley at gmail.com> wrote:
>
>> Hi,
>>
>> I had a look at your topologies and don't see anything obviously wrong.
>> You may however want to check that none of the bonded parameters involving
>> the transformed atoms are being unexpectedly perturbed. (In our setups we
>> typically repeat the A state parameters in the B state column for the
>> bonds, angles, and torsions, as some GROMACS versions would otherwise
>> assume these should be turned to zero).
>>
>
> I read in the gromacs manual that with opls_aa force field, the state A
> parameters are repeated if you dont have any parameters for state B. I saw
> those warnings issued by grompp and in order to avoid any problems I did
> repeat the state A parameters in state B columns as you do it.
>
>
>>
>> Sorry I can't be more help. I'm pretty swamped on this end.
>>
>>
>>
>>
> I once again, I thank you very much for the time you took for my problem.
>
> Regards
> Sai
>
>
>
>
>>
>> On Tue, Oct 30, 2012 at 10:35 AM, Sai Kumar Ramadugu <sramadugu at gmail.com
>> > wrote:
>>
>>> Prof Mobley,
>>>
>>> I just want to add that in the topology file named Ala_charge_off.top, I
>>> am *turning off* the charges on atoms of alanine 40 residue. The total
>>> integral of <dV/dl> from 0 to 1 gave me a value of 274.73 kJ/mol and I
>>> took the negative of this as this should be same as *turning on* the
>>> charges. To justify my approach I did add the partial charges on Ala40
>>> residue in a separate set of simulations (which finished today morning) and
>>> the total integral of <dV/dl> from 0 to 1 gave me a value of -274.705
>>> kJ/mol.
>>>
>>>
>>> Regards
>>> Sai
>>>
>>>
>>> On Sat, Oct 27, 2012 at 11:00 AM, Sai Kumar Ramadugu <
>>> sramadugu at gmail.com> wrote:
>>>
>>>> Dear Prof Mobley,
>>>>
>>>> I have attached the topologies for all three steps in presence of
>>>> ligand. The case for the absence of ligand is same.
>>>> The Glutamate happens to be residue 40 in the topology file (around
>>>> line 615).
>>>> In the case of LJ transformation, for the bonds with dummies, I am
>>>> using the same bond length, angles and dihedral angles as it is for their
>>>> respective original atoms.
>>>> I am using the two papers by Martin Karplus to justify this approach.
>>>> "The role of bonded terms in free energy simulations 1: Theoretical
>>>> Analysis", JPCA, 1999, 103, 103-118.
>>>>
>>>>
>>>> On Fri, Oct 26, 2012 at 1:32 PM, David Mobley <dmobley at gmail.com>wrote:
>>>>
>>>>> Hi,
>>>>>
>>>>> A few responses below. I'm afraid I don't have any brilliant insights
>>>>> right now but there are some things to look at...
>>>>>
>>>>> On Thu, Oct 25, 2012 at 11:05 AM, Sai Kumar Ramadugu <
>>>>> sramadugu at gmail.com> wrote:
>>>>>
>>>>>> Dear Prof Mobley,
>>>>>> I am sorry for the delay in answering your questions to my original
>>>>>> post.
>>>>>> The reason being, I was doing some calculations to answer some of
>>>>>> your questions.
>>>>>>
>>>>>> *My question*: I am getting negative values for vdw transformation
>>>>>> in gromacs when I am mutating sidechain of Glu to Ala. The overall free
>>>>>> energy difference should be positive but I am getting negative and this is
>>>>>> mostly effected due to the negative free energy values from LJ calculation.
>>>>>> Is my LJ transformation correct?
>>>>>>
>>>>>
>>>>> I don't know, honestly. If you're getting a sign that's different from
>>>>> what you expect, it certainly suggests a problem, though it's always
>>>>> possible that the force field disagrees with experiment. At this point,
>>>>> though, I'd be looking for problems with your topologies if I were you (I
>>>>> can't do this since you haven't sent them).
>>>>>
>>>>>
>>>>>> The values themselves:
>>>>>> delG1 (in presence of ligand) Step1 = 1041.31 kJ/mol; Step2: =
>>>>>> -27.08 kJ/mol; Step3 = -274.73 kJ/mol Total = 739.5
>>>>>> delG2 (in absence of ligand) Step1 = 1037.38 kJ/mol; Step2: =
>>>>>> -20.93 kJ/mol; Step3 = -271.86 kJ/mol Total = 744.59
>>>>>>
>>>>>> The relative free energy difference: delG1 - delG2 = 739.5 - 744.59 =
>>>>>> -5.09 kJ/mol. For the transformation of Glu to Ala this value has to be
>>>>>> positive.
>>>>>>
>>>>>> I'm skeptical that the differences in step 1 and step 3 are
>>>>> necessarily meaningful, though they do seem to cancel out. I think you're
>>>>> right to focus on step 2.
>>>>>
>>>>>>
>>>>>> My thermodynamic cycle:
>>>>>> delG1
>>>>>> Protein-Glutamate-Ligand--------------> Protein-Alanine_ligand
>>>>>> | |
>>>>>> |delG3 |delG4
>>>>>> | |
>>>>>> Protein-Glutamate ----------------> Protein-Ala
>>>>>> delG2
>>>>>> I am using delG1 - delG2 for the relative binding free energy of
>>>>>> ligand bound to protein with glutamate to protein with alanine.
>>>>>>
>>>>>> Initially I did a two step approach i.e., I have been doing the
>>>>>> uncharging and charging mutation of Glu--->Ala in one step.
>>>>>> Later I read a post by Prof Shirts in response to Fabian in this
>>>>>> users-list and followed a three step approach. (
>>>>>> http://gromacs.5086.n6.nabble.com/FEP-td4931063.html). These
>>>>>> calculations just started when you responded. My apologies!
>>>>>>
>>>>>> Step1: Uncharging of Glu
>>>>>> Step2: Vdw mutation of side chain of Glu to Ala
>>>>>> Step3: Charging of Ala
>>>>>>
>>>>>> If you want me to troubleshoot, it'd be nice to see topologies for
>>>>> these, especially step 2, but actually all of them. Most problems result
>>>>> from not setting up the right transformation.
>>>>>
>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>> The topology section of mutating the potassium ion K+ to K0 atom is
>>>>>>>> pasted below:
>>>>>>>>
>>>>>>>> [ moleculetype ]
>>>>>>>> ; Name nrexcl
>>>>>>>> KM 1
>>>>>>>>
>>>>>>>> [ atoms ]
>>>>>>>> ; nr type resnr residue atom cgnr charge mass
>>>>>>>> typeB chargeB massB
>>>>>>>> 1 opls-408 1 KM KM 1 1.0
>>>>>>>> 39.0983 DUM_408 0.0 39.0983
>>>>>>>>
>>>>>>>>
>>>>>>>> I am doing this calculation in presence and absence of the ligand.
>>>>>>>>
>>>>>>>> I take it your position restraints are imposed via posre.itp or
>>>>>>> something similar? Also, I hope you are holding the ion and protein apart
>>>>>>> from one another, rather than just holding the ion fixed and allowing the
>>>>>>> protein to drift? How exactly are you handlign restraints?
>>>>>>>
>>>>>>
>>>>>> For the position restraints, I am using the restraints below the
>>>>>> [molecule] and [atoms] directive with ifdef command and giving a value of
>>>>>> 1000 for fx, fy and fz. In the production mdp input file, I am using
>>>>>> -DPOSRES to implement the same. I am pasting the relevant part of the
>>>>>> topology for position restraints.
>>>>>>
>>>>>> #ifdef POSRES_ION
>>>>>> ; Position restraint for each Potassium ion
>>>>>> [ position_restraints ]
>>>>>> ; i funct fcx fcy fcz
>>>>>> 1 1 1000 1000 1000
>>>>>> #endif
>>>>>>
>>>>>>
>>>>>> Further I am not restraining the protein, but in all the cases, the
>>>>>> protein is not drifting during the simulation time I have which is 2ns. I
>>>>>> have checked this with distance between the atoms of mutated amino acid and
>>>>>> the mutating ion (and also other ions in the system).
>>>>>>
>>>>>
>>>>> Well, that's fine as far as it goes I suppose, but you've set up a
>>>>> case where the ion will drift with respect to the protein (as the protein
>>>>> moves) which will likely cause convergence problems in longer simulations.
>>>>> I don't think this is the problem, but it's certainly A (probably minor)
>>>>> problem.
>>>>>
>>>>
>>>> I am thinking to restart my simulations with restraints on the C-alpha
>>>> atoms of the protein.
>>>>
>>>>>
>>>>>
>>>>>>
>>>>>>>
>>>>>>>> Step 2: LJ transformation calculation
>>>>>>>>
>>>>>>>> The relevant part of the topology for LJ transformation is pasted
>>>>>>>> below:
>>>>>>>>
>>>>>>>> ; residue 40 GLU rtp GLU q -1.0
>>>>>>>> 552 opls_238 40 GLU N 191 0.0
>>>>>>>> 14.0067 ;
>>>>>>>> 553 opls_241 40 GLU H 191 0.0
>>>>>>>> 1.008 ;
>>>>>>>> 554 opls_224B 40 GLU CA 191 0.0 12.011
>>>>>>>> ;
>>>>>>>> 555 opls_140 40 GLU HA 191 0.0 1.008
>>>>>>>> ;
>>>>>>>> 556 opls_136 40 GLU CB 192 0.0 12.011
>>>>>>>> ;
>>>>>>>> 557 opls_140 40 GLU HB1 192 0.0 1.008
>>>>>>>> ;
>>>>>>>> 558 opls_140 40 GLU HB2 192 0.0 1.008
>>>>>>>> ;
>>>>>>>> 559 opls_274 40 GLU CG 193 0.0 12.011
>>>>>>>> opls_140 0.0 1.008 ; qtot 0.78
>>>>>>>> 560 opls_140 40 GLU HG1 193 0.0 1.008
>>>>>>>> DUM_140 0.0 1.008 ; qtot 0.84
>>>>>>>> 561 opls_140 40 GLU HG2 193 0.0 1.008
>>>>>>>> DUM_140 0.0 1.008 ; qtot 0.9
>>>>>>>> 562 opls_271 40 GLU CD 194 0.0 12.011
>>>>>>>> DUM_271 0.0 12.011 ; qtot 1.6
>>>>>>>> 563 opls_272 40 GLU OE1 194 0.0 15.9994
>>>>>>>> DUM_272 0.0 15.9994 ; qtot 0.8
>>>>>>>> 564 opls_272 40 GLU OE2 194 0.0 15.9994
>>>>>>>> DUM_272 0.0 15.9994 ; qtot 0
>>>>>>>> 565 opls_235 40 GLU C 195 0.0
>>>>>>>> 12.011 ;
>>>>>>>> 566 opls_236 40 GLU O 195 0.0 15.9994
>>>>>>>> ;
>>>>>>>>
>>>>>>>> Since I am only changing the side chain of Glu to side chain of
>>>>>>>> Ala, I have mutated only the part of the side chain that is different
>>>>>>>> between the two amino acids.
>>>>>>>>
>>>>>>>
>>>>>>> Have you now turned off the LJ interactions on your dummy atoms? Or,
>>>>>>> when do you finish turning off your atoms?
>>>>>>>
>>>>>>
>>>>>> I have turned off the sigma and epsilon at this stage in the
>>>>>> ffnonbonded.itp file. For all the dummy atoms in my system, sigma and
>>>>>> epsilon are zero.
>>>>>>
>>>>>
>>>>> Didn't you say in your LJ transformation you are turning Glu into Ala?
>>>>> Just to confirm, you are turning a carbon into a hydrogen at this stage as
>>>>> well as changing atoms into dummies, right?
>>>>>
>>>>
>>>> Yes, I am mutating the C-gamma to Hydrogen and the rest of the side
>>>> chain atoms of Glutamate to dummies.
>>>>
>>>>>
>>>>>>
>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>> At the same time I am also mutating the K0 atom to a dummy atom.
>>>>>>>> The topology section of mutating the potassium atom to a dummy atom is
>>>>>>>> pasted below:
>>>>>>>>
>>>>>>>> [ moleculetype ]
>>>>>>>> ; Name nrexcl
>>>>>>>> KM 1
>>>>>>>>
>>>>>>>> [ atoms ]
>>>>>>>> ; nr type resnr residue atom cgnr charge mass
>>>>>>>> typeB chargeB massB
>>>>>>>> 1 opls-408 1 KM KM 1 0.0
>>>>>>>> 39.0983 DUM_408 0.0 39.0983
>>>>>>>>
>>>>>>>>
>>>>>>>> I am doing this transformation in presence and absence of the
>>>>>>>> ligand.
>>>>>>>>
>>>>>>>> After these two steps:
>>>>>>>> For the analysis I am just using the values of dV/dl printed for
>>>>>>>> every 10 steps from the simulation from 0 to 1 in lambda and integrating
>>>>>>>> the dV/dl w.r.t. lambda.
>>>>>>>>
>>>>>>>> Step 1 charge in presence of Ligand = 790.109 kJ/mol
>>>>>>>> Step 2 vdw in presence of Ligand = -29.324 kJ/mol
>>>>>>>>
>>>>>>>> The sum of two steps in presence of ligand = 760.785 kJ/mol
>>>>>>>>
>>>>>>>> Step1 charge in absence of Ligand = 787.33 kJ/mol
>>>>>>>> Step 2 vdw in absence of Ligand = -21.8127 kJ/mol
>>>>>>>>
>>>>>>>> The sum of two steps in absence of ligand = 765.517 kJ/mol
>>>>>>>>
>>>>>>>> I'm confused. It looks to me like you should be taking three steps
>>>>>>> Glu -> Ala: Turn off charges on atoms in Glu being deleted while turning
>>>>>>> off charges on ion; second, turn off LJ on any atoms in Glu being deleted;
>>>>>>> third, change any charges on atoms retained in Ala from what they were in
>>>>>>> Glu to what they should be in Ala. I think you are missing some of that
>>>>>>> unless there's something you haven't explained.
>>>>>>>
>>>>>>
>>>>>> Here is the part where I did only two steps and later read a post by
>>>>>> Prof Shirts and started doing the three step approach as you wrote. The
>>>>>> delay for the reply is because of doing new set of calculations. For step 1
>>>>>> I have topology and coordinates corresponding to Glu residue and for step 3
>>>>>> I have topology and coordinates for Ala residue and in both cases I am
>>>>>> changing only the partial charges on the atoms. Did I misunderstand
>>>>>> anything here? For step 2, I have topology corresponding to Glu and I am
>>>>>> deleting the additional atoms in glutamate to make it alanine.
>>>>>>
>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>> The relative free energy of the mutation Glu-->Ala = -4.732 kJ/mol
>>>>>>>>
>>>>>>>> My main concern is with respect to the LJ transformation. Is my
>>>>>>>> approach correct in terms of modifying the side chain of glutamate to
>>>>>>>> alanine? The doubt arises because the relative free energy difference is
>>>>>>>> negative where as the experimental value is close to 20 kJ/mol. I am way
>>>>>>>> under-predicting the value and with a negative sign.
>>>>>>>>
>>>>>>>
>>>>>>> You haven't discussed how you get the total free energy of this. I
>>>>>>> assume you're doing some solvent calculation as well?
>>>>>>>
>>>>>>
>>>>>> I am using the <dV/dl> vs lambda for lambda 0 to 1 at various lambda
>>>>>> values ( I have 38 lambdas for vdw calculations and 21 lambda values for
>>>>>> charging/uncharging calculations) and integrating the curve from 0 to 1. I
>>>>>> have tried the gbar approach too. But my integrated values and gbar values
>>>>>> for free energy were almost same. But I can use gbar if needed.
>>>>>>
>>>>>> I am using the thermodynamic cycle as shown above (similar to gromacs
>>>>>> manual figure 3.10 A). What additional solvent calculation do I need to do
>>>>>> when I am taking the difference delG1 - delG2?
>>>>>>
>>>>>
>>>>> Sorry, now that I see your thermodynamic cycle I see you don't need a
>>>>> "solvent" calculation (or rather you have a "solvent" calculation where
>>>>> it's just the protein in solution with no ligand).
>>>>>
>>>>>>
>>>>>>
>>>>>>> Thanks.
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>> When I did these same calculations using AMBER 11, I am getting
>>>>>>>> 6.945 kcal/mol which is still less but it does not have a negative sign.
>>>>>>>> The other thing I observed is the coulomb and the vdw dV/dl vs
>>>>>>>> lambda curves for OPLS-AA (gromacs) and AMBER 99SB (AMBER11) have a very
>>>>>>>> similar trend only shifted in the values of dV/dl in the y-axis.
>>>>>>>>
>>>>>>>> I can attach the graphs and include more details, if needed.
>>>>>>>> Let me know.
>>>>>>>>
>>>>>>>>
>>>>>>>> Thanks for your time,
>>>>>>>>
>>>>>>>> Regards
>>>>>>>> Sai
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> --
>>>>>>> David Mobley
>>>>>>> dmobley at gmail.com
>>>>>>> 949-385-2436
>>>>>>>
>>>>>>>
>>>>>>
>>>>>
>>>>>
>>>>
>>>> Thank you for your time.
>>>>
>>>> Regards
>>>> Sai
>>>>
>>>>>
>>>>> --
>>>>> David Mobley
>>>>> dmobley at gmail.com
>>>>> 949-385-2436
>>>>>
>>>>>
>>>>
>>>
>>
>>
>> --
>> David Mobley
>> dmobley at gmail.com
>> 949-385-2436
>>
>>
>
More information about the gromacs.org_gmx-users
mailing list