Archive for the ‘Biophysics’ Category

Thermodynamics integration

Sunday, January 16th, 2011

In free energy calculation, thermodynamics integration, which is alternative to free energy perturbation has generated some steam recently. This post is a brief introduction to this method. Free energy is a state function, if there are two states, state 0 and state 1, the free energy of transferring from state 0 to state 1 is written as:

 \Delta G(0 \rightarrow 1) = -k_{B}T\ln \frac{Q_{1}}{Q_{0}}=-k_{B}T\ln\frac{\int{e^{-\beta H_{1}(x,p)}dxdp}}{\int{e^{-\beta H_{0}(x,p)}dxdp}}

Since free energy is independent of actual path, we introduce an alchemical order parameter, \lambda, the Hamiltonian of the system is then H_{\lambda}(x,p), and we need make sure that

H_{\lambda}(x,p)=H_{0}(x,p) when \lambda=0, and H_{1}(x,p) when \lambda=1.

So we have free energy dependent on \lambda,

\Delta G(0 \rightarrow \lambda)=-k_{B}T\ln \frac{Q_{\lambda}}{Q_{0}}. Its derivative is

 \frac{d\Delta G}{d\lambda} = \frac{\int {\frac{\partial H_{\lambda}}{\partial \lambda}e^{-\beta H_{\lambda}(x,p)}dxdp}}{\int e^{-\beta H_{\lambda}dxdp}} =\ll \frac{\partial H_{\lambda}}{\partial \lambda} \gg_{\lambda}

The free energy between 0 and 1 is

\Delta G(0 \rightarrow 1)=\int_{0}^{1}\frac{d\delta G(\lambda)}{d\lambda}d\lambda = \int_{0}^{1}\ll \frac{\partial H_{\lambda}}{\partial \lambda} \gg_{\lambda}d\lambda \approx \sum_{k} \ll \frac{\partial H_{\lambda}}{\partial \lambda} \gg_{\lambda}\Delta \lambda_{k}.

In general, there is a thermodynamics cycle to carry out the actual free energy difference. For example, in solvation free energy calculation, there is such a cycle:

[Solute]_{vacuum} \longrightarrow [solute]_{solvent} : \Delta G_{solvation}

\Downarrow \Delta G^{0}_{annihilation}                 \Downarrow \Delta G^{1}_{annihilation}

[Nothing]_{vacuum} \longrightarrow [Nothing]_{solvent} :\Delta G=0

From the above cycle, we have \Delta G_{solvation}= \Delta G^{0}_{annihilation}-\Delta G^{1}_{annihilation}.

Generalized acceptance ratio

Wednesday, December 8th, 2010

The generalized acceptance ratio scheme has been used intensively in MD to get free energy difference between two states, in particular, in thermodynamic integration or free energy perturbation simulations. It is post simulation analysis tool and uses a combination of forward and backward sampling to minimize the free energy error. The idea was introduced by Charles H. Bennett in 1976 and the original version is called Bennett Acceptance Ratio (BAR), here we focus on the generalized version, though very subtle difference.

Suppose we have a system at two states, 0 and 1, the transferring energy for from state 0 to state 1 is

 \Delta A = \beta^{-1} ln \frac{Q_{0}}{Q_{1}}

where Q_{0} and Q_{1} are the partition functions of state 0 and 1, respectively.

Let’s do a little trick here

\frac{Q_{0}}{Q_{1}}=\frac{Q_{0}}{Q_{1}}\frac{\int d\vec{r}^{N} w(\vec{r}^{N})e^{-\beta (U_{0}+U_{1})}}{\int  d\vec{r}^{N} w(\vec{r}^{N})e^{-\beta (U_{0}+U_{1})}}

=\frac{Q_{0}}{Q_{1}} \frac{\int d\vec{r}^{N} w(\vec{r}^{N})e^{-\beta U_{0}}e^{ -\beta U_{1}}}{\int d\vec{r}^{N} w(\vec{r}^{N})e^{-\beta U_{0}}e^{ -\beta U_{1}}}

=\frac{\ll w e^{-\beta U_{1}}\gg_{0}}{\ll w e^{-\beta U_{0}}\gg_{1}}

=e^{\beta \Delta A} [Eq.1]

So we introduce an arbitrary function w into the free energy difference calculation, and so far nothing gain and nothing loss. Before we continue to work on this new and trivial formula, we need fresh up our knowledge of error estimation using delta function.  Suppose we have a function of measurement f(X_{N}) of a sample X_{N} consisting of N independent configurations, then the expectation value of f(X_{N}) can be written in the following expansion manner:

E[f(X_{N})] = f(\overline{X}) + \sum_{j=1}^{n}\frac{f^{j}(\overline{X})}{j!}E(X_{N}-\overline{X})^{j}+O(N^{-(n+1)/2}) where f is bounded, and so be the first n+1 derivatives, \overline{X} is the true value of the average X and f^{j}(\overline{X}) denotes the jth derivative of f at \overline{X}. If we take only the first order, we have

E[f(X_{N})] = f(\overline{X}) + \frac{\partial f}{\partial X}\vert _{X=\overline{X}} E(X_{N}-\overline{X})

the error for f(X_{N}) is then written as:

\sigma^{2}_{f(X)}=[\frac{\partial f}{\partial X}]^{2}\sigma^{2}(X) [Eq.2].

Now let’s use Eq.2 to analyze the error in Eq.1 for the free energy difference. Very simple substitution and derivation leads to:

\sigma^{2}_{\beta \Delta A} = \frac{\ll [w e^{-\beta U_{1}}]^{2}\gg_{0}-\ll w e^{-\beta U_{1}}\gg_{0}^{2}}{N_{0}\ll w e^{-\beta U_{1}}\gg_{0}^{2}} + \frac{\ll [w e^{-\beta U_{0}}]\gg_{1}-\ll w e^{-\beta U_{0}}\gg_{1}^{2}}{N_{1}\ll w e^{-\beta U_{0}}\gg_{1}^{2}} [Eq.3]

where N_{0} and N_{1} are the numbers of independent configurations in state 0 and 1, respectively. Notice that the two denominators in the right side of the equation are from [\frac{\partial f}{\partial X}]^2 as shown in Eq.2.

Let us fill out the integration in Eq.3:

\sigma^{2}_{\beta \Delta A}= \int d \vec{r}^{N} (\frac{Q_{0}}{N_{0}}e^{-\beta U_{1}} + \frac{Q_{1}}{N_{1}}e^{-\beta U_{0}} ) w^{2} e^{-\beta (U_{0}+U_{1})} \frac{1}{(\int d \vec{r}^{N} w e^{-\beta (U_{0}+U_{1})})^{2}}-\frac{1}{N_{0}}-\frac{1}{N_{1}}

One observation is that if we multiple w with a constant, the variance doesn’t change at all and since we said before w is an arbitrary function, we define the following constant:

\int d \vec{r}^{N} w e^{-\beta (U_{0}+U_{1})}=constant [Eq.5].

Then minimizing the variance is formulated as minimizing the variance with the Eq.5 equality constraint. We invoke Lagrange multiplier to solve such a minimization problem.

\Lambda(w,\lambda)=\sigma^{2}_{\beta \Delta A}+\lambda (\int d \vec{r}^{N} w e^{-\beta (U_{0}+U_{1})}-C_{1})

\sigma^{2}_{\beta \Delta A}= \frac{1}{C_{1}}\int d \vec{r}^{N}  (\frac{Q_{0}}{N_{0}}e^{-\beta U_{1}} + \frac{Q_{1}}{N_{1}}e^{-\beta  U_{0}} ) w^{2} e^{-\beta (U_{0}+U_{1})} -\frac{1}{N_{0}}-\frac{1}{N_{1}}

The essence of Lagrange multiplier tells us that

 \nabla_{w,\lambda} \Lambda(w,\lambda)=0, this gives us the solution for w:

 w = \frac{C}{\frac{Q_{0}}{N_{0}}e^{-\beta U_{1}} + \frac{Q_{1}}{N_{1}}e^{-\beta  U_{0}}} [Eq.6]

Now we substitute w in Eq.1, we have

 \frac{Q_{0}}{Q_{1}}=\frac{\ll [1+e^{\beta (U_{0}-U_{1}+C)}]^{-1}\gg_{1}}{\ll [1+e^{\beta (U_{1}-U_{0}-C)} ]^{-1}\gg_{0}} e^{\beta C} [Eq.7]

where \beta C=\ln \frac{Q_{0} N_{1}}{Q_{1} N_{0}} and \beta \Delta A = \ln\frac{Q_{0}}{Q_{1}}.

Let us reorganize Eq.7 and make it looking nicer:

 e^{\beta \Delta A} = \frac{\ll [1+e^{\beta (U_{0}-U_{1}+C)}]^{-1}\gg_{1}}{\ll [1+e^{\beta  (U_{1}-U_{0}-C)} ]^{-1}\gg_{0}} e^{\beta C}

C= \Delta A+ \beta^{-1} \ln \frac {N_{1}}{ N_{0}}

Eq.7 is the so called generalized acceptance ratio. It is both forward and backward sampling because the sign difference of \Delta U = U_{1}-U_{0} in the numerator and denominator of right hand side of Eq.7.

The free energy difference can be get using the self-consistent way, that is, since \Delta A shows up in both sides of the nonlinear equation Eq.6, we use iteration to get the actual value for the free energy difference.


1. Free energy calculations: theory and applications in chemistry and biology, Christophe Chipot and Andrew Pohorille, Springer

2.Understanding molecular simulations: from algorithms to applications, D. Frenkel and B. Smit Academic

3.Good practices in free-energy calculations, Andrew Pohorille, Christopher Jarzynski, and Christophe Chipot, J. Phys. Chem. B 2010, 114, 10235-10253

The three references are exceptionally good in their own ways and I absolutely recommend them!

Protein-ligand binding equlibria

Friday, October 1st, 2010

Reversible binding of ligands to proteins is similar to reversible folding of proteins and involves reaching equilibrium. The binding reaction is

Protein \cdot Ligand \leftrightharpoons Protein + Ligand

The dissociation constant which describes the equlibria is

 K_{D} = \frac{[Protein][Ligand]}{[Protein \cdot Ligand]}=\frac{[P][L]}{[P\cdot L]c^{0}}, \ with\ c^{0}=1 mol\ L^{-1}

So when it is reported K_{D}=10^{-9}mol\ L^{-1}, it should be just K_{D}=10^{-9}.

And Gibbs free energy change of the binding is related to the dissociation constant:

\Delta G^{binding} = \Delta H - T \Delta S = -RT\ lnK_{D}

A good reference of this post is from “Introduction to Protein Science: Artchitecture, Function, and Genomics by Arthur M. Lesk”.

Two more interesting papers

Tuesday, August 3rd, 2010

To read:

A three-dimensional model of the yeast genome(link)

Sequence space and the ongoing expansion of the protein universe (link)

Briefs of two experiments: Differential Scanning Calorimetry and Pressure Perturbation Calorimetry

Tuesday, March 2nd, 2010

I happened to read a paper on the measurement of  \frac{\Delta V_{pr}} {V_{pr}} . And the heat capacity measurement was used, so I dig into it and found how it is done.

Two pans are connected to a single heater to make sure the heating rate is exactly same. One pan with sample molecule plus solvent, the other pan is reference pan, so just with same amount of the solvent in the sample pan. Then heat up both pans simultaneously and at exact same speed. Suppose we have T1 and T2 for increased temperatures in the sample and reference pans, respectively. And the total amount of heat for each one is \Delta Q. Then we have the heat capacity of the molecule determined by:

C_{p}=\frac{\Delta Q}{T_{1}}(1-\frac{T_{1}}{T_{2}})

This is essentially what “Differential Scanning Calorimetry” experiment is doing. A good reference that I used is here.

Regarding the Pressure Perturbation Calorimetry experiment. We need start with the thermodynamics laws.

We have \Delta Q_{reversible} = T \Delta S , therefore (\frac{dQ_{rev}}{dP})_{T} = T(\frac{dS}{dP})_{T} .

Maxwell relation tells us that

(\frac{dS}{dP})_{T}=-(\frac{dV}{dT})_{P} (1).

So we have:

(\frac{dQ_{rev}}{dP})_{T} = -T (\frac{dV}{dT})_{P} or (\frac{dQ_{rev}}{dP})_{T} = - TV\alpha (2), where,

 \alpha = \frac{1}{V}(\frac{V}{T})_{P}.

Integrate Eq. 2, we have:

Q_{rev} = - T V \alpha \Delta P.

The Pressure Perturbation Calorimetry measures the thermal expansion coefficient \alpha. Same as in DSC, there two pans, one is for the sample and the other is for the reference. Suppose the heat flow into two pans are Q_{sample}, Q_{reference}. The total volume of the sample pan is V_{total} = g_{0}V_{0} + g_{s} V_{s} , and g_{0} is the mass of solvent and V_{0} is its volume, g_{s} is the mass of the sample. So the PPC measures Q_{sample} = -T(g_{0}V_{0} \alpha_{0} + g_{s} V_{s} \alpha_{s})\Delta P , and Q_{reference}=-T(g_{0}V_{0} \alpha_{0} + g_{s} V_{s} \alpha_{0})\Delta P. Simplify both equations we have:

\alpha_{s}=\alpha_{0}-\frac{\Delta Q}{Tg_{s}V_{s}}, \Delta Q = Q_{sample}-Q_{reference}.

That is how PPC measures the thermal expansion coefficient. A good reference that I used is here.