Non equilibrium method and AFM simulation w gromacs

Atomic Force Microscopy

Atomic Force Microscopy (AFM) encompasses a variety of techniques related to the use of an extremely small mechanical probe - the probe tip can have a radius as small as 2 nm - that interacts physically with the sample of interest. In addition, the probe tip can be modified with biomolecules of interest, such as single-stranded DNA, that will interact with other biomolecules on a surface. One measurement an AFM can make is the force required to pull the tip away from a surface - in theory, this can be used to obtain the binding energy between two molecules.

This picture is an example of an AFM tip that is attached to the end of a cantilever.

This figure shows an average of 20 withdrawal traces between two complementary DNA strands. One strand is attached to an AFM tip, and another strand is attached to an immobile substrate.

---

Comparison of AFM experiment to Steered Molecular Dynamics

Steered Molecular Dynamics (SMD) refers to simulations in which an external force is applied to a group of atoms or molecules, closely resembling the real physical system of AFM experiments. However, due to the time constraints of MD, SMD will proceed 6 to 9 orders of magnitude faster than AFM experiments. To review basic thermodynamics, when work is being done on a system, the ensemble average of the work is related to the change in the Helmholtz free energy by

${\overline {W}}\geq \Delta F$

In this case ${\overline {W}}$ denotes an average over an ensemble of measurements of W.

When the work is done infinitely slowly, the system is allowed to equilibrate at each step leading to the equality

$\ W=\Delta F$

This is a somewhat reasonable model for AFM experiments in which pulling velocities are on the order of 1 micron/s, but time constraints prevent this in SMD, in which pulling velocities are currently $\geq$ 0.1 m/s. In order to compare SMD with AFM experiments, methods must be developed in order to differentiate between the equilibrium free energy change and energy lost due to operating in the non-equilibrium regime. Friction is one major source of energy loss in SMD that is less important in AFM. It should be noted that AFM experiments must be carefully designed so as not to violate this quasi-equilibrium assumption, or the measured force as a function of separation distance will also be too large.

In fact, the average AFM trace shown above was not acquired in the equilibrium regime. Another group later went back and applied the Jarzynski equality (described below) to recover the real equilibrium force data:

AFM pulling in GROMACS

The implementation details are given in the GROMACS manual version 3.3 in section 6.1.

Briefly, GROMACS requires four additional files: two input files and two output files

-pi pull.ppa

If this file is specified the pull code will be used. It contains the parameters that control what type of calculation is done. A full explanation of all the options is given below.

-pn index.ndx

This file defines the different groups for use in all pull calculations. The groups are referred to by name, so the index file can contain other groups that are not used as well.

-po pullout.ppa

A formatted copy of the input parameter file with the parameters that were actually used in the run.

-pdo pull.pdo

The data file with the calculated forces (AFM pulling, constrained) or positions (umbrella sampling).

AFM pulling attaches a spring to the center of mass of the pulled group. This spring moves in relation to the center of mass of the reference group. It the reference group moves, then the spring does as well. If the reference group is left blank, then pulling will be in absolute coordinates, with the reference group position taken to be [0 0 0].

afm rate1 =
afm rate2 =

The rate at which the spring moves in nm/ps for each group.

afm k1 =
afm k2 =

The force constant of the spring for each pulled group in kJ mol−1 nm−2.

afm dir1 =
afm dir2 =

Unit vector describing the direction of pulling.

afm init1 =
afm init2 =

Vector describing the initial position of the spring relative to the reference group. To start a simulation with zero initial force on the pulled group, the initial position should be set to the position of the pulled group relative to the reference group.

AFM

10.5 3.5 3.63 3.6

The output file consists of the time, the position of the reference group, and the position of the pulled group, and the position of the spring. The position of the pulled group and spring are in absolute coordinates. Subtract the reference coordinate to get the position relative to the pulled group. Note that only the dimensions for which pulldim is Y are output.

Jarzynski Equality

One goal of ligand-protein AFM experiments and SMD simulations is to obtain the potential of mean force (PMF), which gives the free energy change as a function of separation distance between two entities of interest. Sample PMFs of the interaction between glucose and benzene, in vacuum and water, are shown below.

However, as indicated previously, SMD will introduce energy dissipation during the pulling simulation, and so the work done by the virtual spring will not directly relate to the free energy of unbinding (or whatever process is being studied). The Jarzynsky equality is one method of obtaining the free energy from the simulation data:

This states that the equilibrium free energy difference $\Delta F$ between states A and B can be obtained by an exponentially weighted average over all possible realizations of the external process that moves the system from the equilibrium state A to the same external conditions as that of the equilibrium state B.

This is an important equality for both experimental and computational results. The equality was mathematically supported and extended by Hummer and Szabo to deal with experimental AFM results:

While the technique appears to be well-suited for use with AFM results in which the work done at any time is on the order of $k_{B}T$, it is not as clear that it can be easily extended to simulation data. The limitations of this equality have been stated in several ways:

• In the initial derivation, Jarzynski states that the equality only holds for cases in which the system of interest is only weakly coupled to the reservoir into which excess energy is dissipated.
• Others have stated that the equality is limited to microscopic applications where the magnitude of thermal fluctuations are significant with respect to the pulling force.
• Another paper indicates that the variation in the work function must be $\leq k_{B}T$ in order to avoid statistical error, and further states that this is difficult to accomplish via simulation.

Langevin Modeling

An alternative method to recover the PMF relies on a form of the Langevin equation. This is described in several papers; for consistency the following equations are pulled from a single paper which uses Langevin modeling to recover the PMF for ions being pulled rapidly through a channel.

This work provides a definition for the PMF of a single ion as:

${PMF}\left(z\right)=-k_{B}Tln\left(P\left(z\right)/P_{0}\right)$

Here $P(z)$ is the equilibrium probability to find an ion at position z, and $P_{0}$ sets the zero point of the PMF, in this case when the ion is in the bulk phase.

The Langevin equation is given by

$\ F_{ext}$ is any externally applied force and $\ \gamma$ is the friction coefficient. $\ \sigma \xi (t)$ is a Gaussian white noise term that is defined as $\left\langle \xi (t)\xi (t')\right\rangle =\delta (t-t')$ and the root mean square amplitude $\ \sigma$. $\ \sigma$ and $\ \gamma$ are related via $\ \sigma ^{2}=2\gamma k_{B}T$

$\ \gamma$ is estimated using a velocity autocorrelation method:

For realistic values of $\ \gamma$, the first term in [6] will dominate by a factor of around 1000. The autocorrelation function then approximates as follows:

As mentioned, this method was used to simulate ions K+, Cl-, and Na+ being pulled through the membrane pore Vpu. Sample autocorrelation curves are shown below. The authors chose an evaluation time of 0.03ps, during which all the trajectories are fairly similar. There is no rigorous support for this evaluation time.

The major results are a set of snapshots of ions being pulled through the pore, and the PMF associated with each ion for the pore in two different configurations with two different spring constants.

The calculated PMFs are shown below.

The authors suggest that these supports weakly support the experimental data which suggest that the pore is weakly cation selective.

A comparison of methods for calculating the PMF

In order to determine which method of reconstructing the PMF is ideal, Gullingrsud et al. performed a comparison of WHAM (umbrella sampling), Gaussian Drift based on van Kampen's $\ \Omega$-expansion, and a least squares minimization of the Onsager-Machlup action with respect to the choice of PMG. The second two methods assume a Langevin description as above in order to account for the nonequilibrium character of the SMD.

These methods were applied to several different theoretical systems with PMFs corresponding to either a sinusoid, two Gaussians, or a step function. The PMFs were estimated for each method and then plotted along with the true potential U(x).