Re: FEP/binding free energy difference WT vs Mut

From: Chris Chipot (
Date: Sun Feb 15 2009 - 20:03:27 CST

Dear Luis,

I would add one comment to Jerome's message. However appealing the
access to free-energy component is, one should be cognizant that
such a decomposition is not always feasible. In particular, let us
write the potential-energy function as a sum, U(λa,λb)=U0+Ua(λa)+
Ub(λb). When turning to the thermodynamic-integration formalism,
the relevant free-energy changes look like ΔAa=∫<Ua>dλa with λb=1 and
ΔAb=∫<Ub>dλb with λa=1, which you can see as, say the electrostatic
and the van der Waals components of the potential energy function.
Now, let us think differently: λ=λa=λb, and ΔAa'=∫<Ua>dλ and ΔAb'=
∫<Ub>dλ. Here, the total free-energy change is ΔA=ΔAa'+ΔAb'. If we
go back to the first definition, we have:

ΔAa=∫{∫Ua exp[-β(U0+λaUa+Ub)]dx}/{exp[-β(U0+λaUa+Ub)]dx}dλ and
ΔAb=∫{∫Ub exp[-β(U0+Ua+λbUb)]dx}/{exp[-β(U0+Ua+λbUb)]dx}dλ

which is at variance with:

ΔAa'=∫{∫Ua exp[-β(U0+λUa+λUb)]dx}/{exp[-β(U0+λUa+λUb)]dx}dλ and
ΔAb'=∫{∫Vb exp[-β(U0+λUa+λUb)]dx}/{exp[-β(U0+λUa+λUb)]dx}dλ

We have actually discussed this point from the perspective of
2nd-order perturbation theory:


Yet, with the alternate definition, we have:


When you sum up the two contributions to obtain the total free-energy
change, you can easily see that things don't quite match at the level
of the correlation terms, (<UaUb>0-<Ua>0<Ub>0), which ought to remind
us that if free energy is a state function, its components are not.
The final result is, therefore, inherently path dependent.


Jerome Henin a écrit :
> Dear Luis,
> Yes, the thermodynamic cycle you designed is correct. Just in case, I
> will add that in the dissociated state (computing deltaG2), there is
> no need to include the unbound DNA in the actual simulation, as it
> should not be interacting with the protein anyway.
> Now let me try to answer your specific questions:
> 1) If you can afford it, yes, it is better to use PBC with PME.
> However, in some cases, using a solvent bubble may be an acceptable
> simulation environment. This is not my field of expertise, so I leave
> it to others to comment on this if they want to.
> 2) If you don't use PME, there is no need to keep the system globally
> neutral (see point 1). It you do, what I would recomment is having a
> chloride counterion vanish (or a sodium appear) while the arginine
> disappears. Ideally this counterion should have a positional restraint
> ("constraint" in NAMDese) to keep it somewhere in bulk solution.
> 3) You can obtain individual free energy contributions if you use the
> TI implementation currently available in the development version of
> NAMD (CVS). I would recommend that version in any case, as in contains
> many improvements to the code in general and free energy calculations
> in particular.
> Best,
> Jerome
> On Fri, Feb 13, 2009 at 5:43 PM, Luis Cunha <> wrote:
>> Dear NAMD users,
>> I need your advice. I use NAMD to study ligand binding to my protein, but so
>> far I've only been interested in conformation changes, not energies. I need
>> advice with the following problem:
>> I have an xray structure of a protein (actually a domain of the protein)/DNA
>> complex. In this complex, an arginine makes a Hbond with a specific DNA
>> base. We've found patients with a specific genetic disease who have a
>> mutation (Arg to Cys) of this very arginine. Obviously this mutation has
>> profound effects in the way the protein interacts with the DNA, hence it
>> causes a disease. I'm hypothesizing that the Hbond removal in the mutant
>> severely affects the binding capacity. An assay is being developed by
>> someone else to test this in vitro. Meanwhile, I would like to estimate the
>> difference in the binding energy between the mutant and wild type proteins.
>> Please let me know if I got this correctly: What I need is a thermodynamic
>> cycle:
>> deltaG1
>> WT-DNA ---------> WT + DNA
>> | deltaG4 | deltaG2
>> | |
>> MUT-DNA --------> MUT + DNA
>> DeltaG3
>> And the deltadeltaG (the difference in free energy of binding between the
>> mutant and the wild type) that I'm interested in would be given by
>> subtracting deltaG4 from DeltaG2. These 2 deltaGs would be calculated in
>> two alchemical FEP/MD experiments.
>> 1) If this concept is correct, should I do the simulation in solvated
>> systems, with PBC, PME?
>> 2) since the positive charge of the arginine disappears, I guess I need a
>> ion to appear in lambda=1?
>> 3) is it possible to evaluate the individual contributions to the free
>> energy difference? (electrostatic, VdW)
>> thank you in advance to any help provided,
>> Luis

Chris Chipot, Ph.D.
on leave from Nancy Université, CNRS
Theoretical and Computational Biophysics Group
Beckman Institute
University of Illinois at Urbana-Champaign
405 North Mathews                                 Phone: (217) 244-5711
Urbana, Illinois 61801                            Fax:   (217) 244-6078
       To sin by silence when we should protest makes cowards out of men
                                                     Ella Wheeler Wilcox

This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:52:22 CST