[gmx-users] Re: Glutamate to Alanine Mutation
Sai Kumar Ramadugu
sramadugu at gmail.com
Wed Apr 10 16:57:24 CEST 2013
Prof Mobley,
Sorry to bother you, but do you have any suggestions for my previous post?
Thanks again,
Sai
On Sun, Mar 31, 2013 at 9:52 PM, Sai Kumar Ramadugu <sramadugu at gmail.com>wrote:
> 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