Volume 305 of the series Methods in Molecular Biology™ pp 493-515
Force Probe Molecular Dynamics Simulations
- Helmut GrubmüllerAffiliated withTheoretical and Computational Biophysics Department, Max-Planck-Institute for Biophysical Chemistry
10.1007/978-1-59259-912-7_23
Abstract
Many proteins are molecular nano-machines, which perform their biological function via well-coordinated structural transitions. Often, these motions occur on much slower time scales than those accessible to conventional molecular dynamics techniques, which are limited to submicrosecond time scales by current computer technology. This is also true for ligand binding and unbinding reactions. Force probe simulations (or steered molecular dynamics) provide a powerful means to overcome this limitation, and thus to get insight into the atomistic mechanisms that underlie biological functions such as ligand binding. This chapter provides a basic introduction into this method. It further sketches a simple nonequilibrium statistical mechanics treatment that shows how to relate the results of force probe simulations to atomic force microscopy (AFM) or optical tweezer experiments. As an example, enforced unbinding simulations of streptavidin/biotin complexes are detailed.
Key Words
Molecular dynamics simulation protein dynamics force probe simulation steered molecular dynamics ligand/receptor unbinding enforced protein unfolding rupture force calculation unfolding forces rescaling of loading rates nonequilibrium statistical mechanics23.1 Introduction
Chapter 22 provided an excellent introduction into the molecular dynamics (MD) simulation method on biomolecules, which is also the basis for the force probe simulation technique (1) (also known as steered molecular dynamics [2]) described herein. Therefore, the basics of MD simulations are not repeated here in detail, but only briefly sketched. We strongly recommend reading Chapter 22 in combination with this chapter.
A characteristic feature of protein dynamics is the occurrence of collective structural transitions, so-called conformational transitions (3), which are often crucial for ligand binding and for protein function in general. These differ from the thermal high-frequency vibrations, which fall into the picosecond range, in that they occur on much slower time scales from nanoseconds to seconds. The opening and closing of ion channel proteins in nerve cells is an example for such a functional conformational transition. It is characterized by millisecond rates and can be observed via an electric current through the channel by patchclamp measurements (4–6) on single proteins. Other examples are so-called allosteric effects (e.g., in hemoglobin), which are realized by conformational transitions and induced fit motions upon ligand binding, for example, of an antigen to an antibody. Last but no least, we discuss the elementary steps of motor proteins for muscle contraction, which also represent conformational transitions.
Functional processes in proteins occur on a large variety of different time scales, and so does protein dynamics and conformational dynamics that covers a wide hierarchy of dynamical processes: the high frequency part of the thermal fluctuations of single atoms around their average position is characterized by reciprocal frequencies of several ten to several hundred femtoseconds; collective motions of smaller groups of atoms range up to several tens of picoseconds. Larger structural changes occur on a hierarchy of time scales ranging from nanoseconds to hours. There is a convincing amount of experimental and theoretical evidence for hierarchical substates (3) for protein conformations, which are linked by conformational transitions on a wide range of time scales (7). These conformational dynamics apparently are a feature specific to proteins and other systems of comparable complexity (8).
The highly ordered, but heterogeneous structure of proteins in combination with the broad range of nonseparating time scales in protein dynamics renders the application of established concepts of many body physics extremely difficult; work in this direction, therefore, is a challenging field in theoretical physics.
Currently, the explicit, atomically resolved molecular dynamics simulation of the protein of interest, as described in Chapter 22, is still to be considered the only reliable theoretical description. The lack of more elegant treatments has a high price: such protein dynamics simulations are computationally extremely demanding and, therefore with presently available hardware restricted to relatively small simulation systems (several 100,000 atoms) and, most severely, to short time scales clearly below microseconds. Nevertheless, the method is now widely and successfully used for a large number of biological processes that fall into this category, and has yielded a large number of correct predictions and insights into protein function at the atomic level (9).
23.1.1 Protein Dynamics Simulations
Molecular dynamics simulations compute the motion of every single atom within the simulation system (e.g., a protein solvated in physiological solvent, Fig. 1 , left), which is determined by the interaction forces ( Fig. 1 , right) between all the n atoms of the system (10). Typically, these forces are described by a force field V(x 1,⋯, x n ) which serves to compute a trajectory {x i (t j )} i = 1⋯n;j=0 ⋯N via the (classical) Newtonian equations of motion. This trajectory specifies the position x i of each single atom i in the system within a (discretized) period of time, t 0 = 0, t 1=Δt,t 2 = 2Δt,⋯,t N = NΔt, for N so-called integration steps. From such a trajectory one can subsequently calculate the observables of interest, often via averages.
Left: A typical protein dynamics simulation system comprises a protein (dark gray atoms) solvated in water (light gray; only the O-H-bonds are shown). Salt ions (sparse light gray atoms) have been added to the aqueous solvent. Right: The inset shows part of the system, together with selected interatomic forces (for example, chemical binding forces, Coulomb forces between atoms that carry partial charges, Pauli repulsion, van der Waals attraction) that determine the molecular motion.
Interaction contributions to a typical force field. Bond stretch vibrations are described by a harmonic potential V B, the minimum of which is at the equilibrium distance b 0 between the two atoms connected by chemical bond i (see Eq. 1; the indices i, j, etc., are not shown in the figure). Bond angles and out-of-plane angles are also described by harmonic potential terms, V A and V E where Θ0 and κ0 denote the respective equilibrium angles. Dihedral twists are subjected to a periodic potential V D; the respective force coefficients are denoted by k’s with appropriate indices. Nonbonded forces are described by Coulomb interactions, V C, and Lennard-Jones potentials, V LJ = V P + V vdW , where the latter includes the Pauli repulsion, V P ~ r 12, and the van der Waals interaction, V vdW~ -r 6. (Adapted with permission from ref. 23.)
23.1.2 Force Fields for Protein Dynamics
The forces between the atoms of a solvated macromolecule ( Fig. 1 , left) are diverse ( Fig. 1 , inset, right). Chemical binding forces, symbolized by springs, enforce equilibrium distances or angles between chemically bound atoms (small arrows). Pauli repulsion forces (dark arrows) prevent atoms from moving through each other. Long-ranged interactions, particularly electrostatic forces (light arrows) between partially charged atoms (δ+, δ-) contribute substantially to the stability of protein structures and dominate the slow conformational dynamics. Hydrogen bonds, for example, are mainly of electrostatic origin and only to a lesser extent a quantum-mechanical effect, and contribute significantly to the stability of α-helices and β-sheets (22).
All these forces (and several more, which are not discussed here) determine the three-dimensional structure of a protein as well as the motion of each single atom, and therefore have to be considered in protein dynamics simulations. As mentioned before, this is achieved via the force field V(x 1 ,⋯, x n ), which describes the influence of the electronic dynamics on the motion of the nuclei. Most force fields are composed of empirically motivated interaction terms. Construction, parameterization, and testing of a force field is a challenging task, and substantial efforts will be required to further improve their accuracy.
Since the 1970s, several force fields for biomolecules and, specifically, for proteins and DNA/RNA have been developed. Well known and widely used are, for example, CHARMm, Gromos96, AMBER, and OPLS.
Figure 2
shows interaction terms from which these force fields are typically composed. They describe interactions that arise from covalent chemical bonds as well as from noncovalent interactions,
Covalent forces arise from changes of the length of chemical bonds (B) and of the angle between two bonds (A), from twisting chemical bonds (dihedral, D), and deviations of aromatic carbon atoms from their in-plane positions (extraplanar, E). As noncovalent interactions the Coulomb interaction (C), the Pauli repulsion (P) and the van der Waals interaction (vdW) are included; the latter two are conveniently described by a Lennard-Jones potential (LJ).
(A) Principle of single molecule AFM experiments; (B) electron microscopy picture of an AFM cantilever (picture kindly provided by Nanoscope); (C) force probe simulation of a single molecule AFM experiment; (D) typical force histogram obtained from a series of single molecule unbinding experiments and fit to theory (dashed).
Snapshots of enforced dissociation of a biotin molecule (white, thick contours) from the streptavidin binding pocket (only the few relevant residues of the binding pocket are shown as stick-models). Hydrogen bonds are shown as bold dashed lines, and few of the many water bridges with thin, dotted lines. The simulation length is one nanosecond; during this time the biotin molecule is moved by about one nanometer. (Reprinted with permission from ref. 1.)
Force exerted onto the biotin ligand during enforced dissociation (force profile, see Note 2 ). A large number of local force maxima can be seen, each of which can be attributed to the rupture of single noncovalent bonds like hydrogen bonds or water bridges. The bold dashed lines at the top denote the respective residues, which are also shown in Fig. 4 . The thin dashed lines denote rupture of water bridges. (Adapted with permission from ref. 1.).
Measured ([35] and Δ [36]) and calculated (•) unfolding forces of single titin molecules as a function of the loading rate kv. The solid line shows the fit of Eq. 14 to the computed unfolding forces and the spontaneous unfolding rate (◯). As can be seen, the experimental values are predicted very well. The six insets show snapshots of the protein during simulated unfolding. (Adapted with permission from ref. 37.)
Time dependent free energy in force probe experiments and simulations. A time-dependent potential V cant = (x,t), Eq. 2, describes the force exerted by the cantilever and modifies the unperturbed energy landscape G(x). As a result, the barrier that determines the spontaneous dissociation rate k 0 (black arrow) is reduced (left: dashed vertical arrow). For a more complex energy landscape (right), also the dissociation length Δx (dashed horizontal arrow), i.e., the distance between minimum and maximum of the free energy landscape, may change. (Adapted with permission from ref. 37.)
23.1.3 Force Probe Simulations
In recent years, we have seen dramatic improvements in single-molecule experiments, particularly atomic force microscopy (AFM) methods (28,29), described in Chapter 21, which motivated the force probe simulation technique described in this section.
In single-molecule force probe experiments (29,30), the cantilever of the AFM microscope is used as a piconewton force sensor ( Fig. 3B ). The cantilever can be positioned to subnanometer accuracy, and its deflection—measuring the applied force—is detected via a laser beam with equally high precision.
Figure 3A illustrates the experiment, in which a ligand (here, biotin, light gray) is forced by the cantilever to leave the specific binding pocket of the receptor (streptavidin, dark gray). A polymer linker connects the proteins (with empty binding pockets) with a surface (left), as well as several biotin molecules with the tip of the cantilever (right). As the cantilever approaches the surface several biotin/streptavidin complexes form, which dissociate one after the other upon subsequent retraction of the cantilever. Occasionally, one single complex remains intact until the very end of the experiment, in which case the force required to dissociate this last complex can be measured from the jump of the deflection of the cantilever to zero.
By repeating the experiment several hundred times one obtains a histogram of dissociation forces ( Fig. 3D ). Its maximum denotes the most probable dissociation force, in the shown example about 270 pN. This is the force that the noncovalent binding interaction between ligand and receptor can withstand at the time scale set by the experiment (typically milliseconds to seconds), and thus a measure for the binding strength (but see Note 6 ).
In a similar manner, molecular forces have been measured recently for a number of systems. One example is the force generated by single motor proteins (31) which is used to drive muscle contraction or to transport intracellular vesicles along filaments. Another one is the force generated by DNA polymerase upon transcription along the DNA primer (32). (For the application of torques rather than linear forces, see Note 4 ; for nondirected forces, see refs. 33,34.)
However, in those experiments, the underlying atomistic dynamics and interactions that generate the measured forces cannot be observed, which is the main motivation to simulate these experiments by means of atomistic protein dynamics simulations ( Fig. 3C ).
To that aim, and modeling the effect of the cantilever, an additional harmonic pulling potential,
is included within the molecular force field (1), such that it acts on that particular atom i of the ligand molecule, which in the real experiment is covalently connected to the cantilever via the polymer linker. This ensures that in the simulation the atom is subjected to the same force as in the experiment. In
Fig. 3
, this pulling potential is symbolized by a spring. The normalized vector n denotes the direction of pulling, and R
i
0
the position of atom i at the start of the simulation. As can be seen from Eq. 2, the minimum zcant = vt of the pulling potential (i.e., the equilibrium position of atom i) is subsequently moved with constant velocity v away from the binding pocket (arrow).
To avoid drift of the protein, the center of mass of the protein (consisting of n
p
atoms with positions R
i
, i = 1⋯n
p
) is kept in place at R
fix by a second harmonic (but stationary) potential,
with suitable spring coefficient k
fix. This ensures that the protein is free to undergo internal motions, such as induced-fit motions upon ligand unbinding, and that it can adopt to the pulling direction by rotations—as in the real experiment—but prevents translational motions.
As an example,
Fig. 4
shows a detailed model of streptavidin/biotin dissociation derived from force probe simulations (see
Note 5
for a caveat), highlighting the sequence and type of stretching and subsequent rupture of single noncovalent bonds between the ligand and the receptor (1,2). These complex sequence of localized rupture events give rise to a complex force profile (
Fig. 5
), that is the exerted force plotted during the simulation,
These force probe AFM simulations of enforced dissociation also allow one to determine dissociation forces from the maximum of the force profile, and to compare them with those measured in AFM experiments. In such comparisons, one has to carefully take into account the fact that the AFM measurements and the force probe simulations are carried out at quite different time scales, namely milliseconds vs nanoseconds. Therefore, in order to compare the respective dissociation forces, the computed force has to be rescaled to the measured ones, for example, for the simplest case using Eq. 14, as shown in Fig. 6 and described later.
This type of force probe simulations, also called steered molecular dynamics, is now widely used, for example, to elucidate at the atomic level the structural changes that are induced by mechanically unfolding proteins like titin (38–44) or by stretching various other polymers (45–47), and to explain the measured forces in terms of intramolecular interactions. In this context, it is worth pointing out that force profiles contain information on the underlying free energy landscape that governs unbinding or unfolding; indeed, information on this landscape can be obtained from force profiles (48–50). For more detailed information, we refer the reader to refs. 37, 51, 52.
23.1.4 Simple Nonequilibrium Statistical Mechanics of Enforced Dissociation
As we have seen in the previous example, both measured and computed dissociation or unfolding forces vary with the time scale at which the process is enforced to occur (1,53,54). Usually, the forces increase with shorter time scales, that is, increasing loading rates. For thermodynamic equilibrium, such time scale dependency should not appear, hence the underlying processes must be assumed to proceed far from equilibrium. This observation motivated the development of nonequilibrium theories of enforced unbinding reactions (2,53–57). One of these shall be sketched here.
We assume an effective free energy G(x) along a (one-dimensional) reaction coordinate x ( Fig. 7 , solid bold line), for example, the distance between the center of mass of the ligand and the one of the receptor, which governs the dissociation reaction.
As sketched in the two panels, the bound state (left minimum) is separated from the dissociated state (right) by a more or less structured energy barrier (global maximum). The height ΔG‡ of this barrier determines the rate coefficient k
0 for thermally activated spontaneous dissociation,
Here, β = 1/(k B T) denotes the reciprocal thermal energy and 0 the attempt frequency or Kramers’ prefactor of the bound state (58–60). For the above streptavidin/biotin dissociation, for example, the reciprocal rate coefficient is several months, whereas it can be as short as milliseconds for antibody-antigen dissociation. Note, that if unbinding is enforced at a time scale that is slower than the reciprocal rate coefficient no dissociation force is measured, because in this case the system has fallen apart already before any significant force has been applied.
As for the force probe simulations described previously, the potential exerted by the cantilever is assumed harmonic and moving with velocity v,
(dashed curve in
Fig. 7
), where k is the effective spring coefficient. Here, we restrict our discussion to the case of small spring coefficients, that is, soft springs, which covers most realistic situations. Therefore, V
cant(x,t) appears in
Fig. 7
as a nearly straight line, the slope of which increases linearly with time. Because the total potential G(x) + V
cant(x,t) now becomes time dependent, so does the barrier height ΔG‡ and hence, the dissociation coefficient k
0. This effect is demonstrated by the thin, solid lines. As can be seen, the barrier height decreases with time. Dissociation occurs when the barrier has become small enough to be overcome thermally activated at the time scale = ΔG‡/ (kvΔx) set by the loading rate kv.
For simple and localized barriers (left picture), we will neglect the slight shift of the position of the minimum that also occurs. In this case, a simple two-state model (53) is appropriate, which neglects the details of G(x) and assumes a linear decrease of the barrier height with time,
This implies a time dependent dissociation coefficient (see
Eq5),
which holds as long as ΔG‡- kvtΔx is larger than k
B
T and back reactions can be neglected. In analogy to the treatment of radioactive decay, the flux across the barrier decreases the probability P(t) to find the system in the bound state at time t,
with P(t = 0) =1. Solution of Eq. 9,
yields a distribution
of dissociation times t
D,
and dissociation forces F
D
= kvt
D
, respectively,
As the main result, the maximum (dP(F
D
)/dF
D
= 0) denoting the most probable dissociation force, reads
(where Eq. 5 has been used) and thus increases logarithmically with the loading rate kv.
This result links the rupture length Δx to the slope of the loading rate dependent dissociation force, and thus provides the basis for dynamic force spectroscopy (Fig. 8A). The larger Δx, the faster the energy barrier is decreased (see Fig. 7 ), and the stronger is the effect of the loading rate kv on the dissociation time t D and dissociation force F D .
Furthermore, this result serves to rescale dissociation forces from force probe simulations to the much slower time scales of the AFM experiments. For the very fast MD time scales, frictional forces can become relevant (which are negligible in the AFM experiments), and in the simplest approach (61) can be heuristically included into Eq. 13,
where γ is the effective friction coefficient. A more rigorous treatment is given, for example, in refs.
2,54.
For more complex energy landscapes as shown in Fig. 7 (right), the simple two-state model is not applicable. For the two barriers shown, there is a critical pulling force for which the left maximum becomes the global one at the time of dissociation (see the thin curve in Fig 7 ), and hence the rupture length Δx jumps to smaller values. Accordingly, the logarithmic slope of the dissociation force increases as shown in Fig. 8B . (For a more general treatment, the reader is referred to refs. 56,57.)
23.2 Materials
Most molecular dynamics packages can be used, as listed in Chapter 22. Good force probe implementations are in GROMACS (62), NAMD (63), and EGO (64). Small to medium-sized simulations can be run on Linux boxes; for large scale simulations we recommend Linux clusters (Beowulf) or special parallel machines such as the IBM SP Series.
Dynamic force spectra show the increase of the dissociation force F max with pulling velocity v of the cantilever or, equivalently, loading rate kv. (A) For the case of a well-localized energy barrier in G(x), Fig. 7 (left), the simple two-state model described in the text predicts a logarithmic increase. (B) For more complex energy landscapes such as the two barriers sketched in Fig. 7 (right), the dynamic force spectrum falls into two regimes. For small v, the larger barrier in G(x) dominates, and the resulting large Δx yields a slow increase of the spectrum with v. For large v, however, the (originally) smaller barrier to the left dominates at the time t D of dissociation, resulting in a steeper logarithmic increase of F max with v.
23.3 Methods
- 1.
Select protein-ligand complex and obtain pdb-structure from the Protein Data Bank (65), http://www.rcsb.org/pdb/, if available. If not, find a homologous structure and try homology modeling. This is a risky business, however, since the accuracy of homology models is generally questionable, and the atomic details that are crucial for the ligand-receptor interactions may be wrong. If the ligand is not present in the protein structure, try docking approaches; the same caveats apply here, however, as for homology modeling. Carefully study the Protein Data Bank entry: is the resolution sufficient (i.e., generally better than 2.5 Å)? Are there missing residues that need to be added by modeling? Have special compounds been added to facilitate crystallization that need to be removed? Is the ligand in the crystal the biologically relevant one, or is it an analog that needs to be changed?
- 2.
For the protein, the force field is typically implemented in the MD package (see Note 1 , however). For the ligand, often not. In this case, the chemical force field parameters can possibly be derived from the structure or from chemically related amino acid motifs. More challenging are the partial charges, which sensitively depend on the chemical environment. Often, it will be necessary to carry out quantum mechanical calculations to get sufficiently accurate charges (see Chapter 24 on density functional (DFT) calculations).
- 3.
Construct the native environment for the protein, typically a water box with solvated salt ions for soluble proteins or a solvated lipid bilayer for membrane proteins. There are several tools available for this task, typically included with MD packages, or stand-alone, for example, SOLVATE (66). The latter takes particular care of the thermodynamically correct placement of the ions according to a Debye-Hückel distribution. As the ion diffusion is rather slow at the MD time scale, correct placement can be crucial for the system to get close to thermal equilibrium. Pay particular attention to the selected system size: in contrast to conventional MD simulations, the solvation box may have to be chosen a bit larger to leave room for the ligand to move away from the protein without hitting the system boundaries or, via periodic boundary conditions, the protein again from the other side. This is particularly true if, rather than ligand unbinding, protein unfolding is studied with force probe simulations. In this case, step-wise increase of the box size as the protein unfolds, or, alternatively repeated truncation of already unfolded segments can help to keep the simulation system size reasonably small, and thus save computer time.
- 4.
Minimize and equilibrate the simulation system as in conventional MD simulations (see Chapter 22). To check whether the system is sufficiently equilibrated, just monitoring energies is typically insufficient. Therefore, always monitor equilibration via the root mean square (rms) deviation from the start structure (see Note 3 ).
- 5.Once a satisfactory simulation start system has been obtained, force probe simulations can be carried out, typically each starting from this equilibrated start structure. Crucial parameters are:
- a.
Pull selection: Which atom or atoms should be subjected to the pulling potential? If a single molecule AFM experiment is to be simulated, a natural choice is that particular atom of the ligand molecule, to which in the AFM experiment the cantilever is linked via a polymer (if known). If the chemistry in the experiment is not known, as a general rule, try various choices and check whether, and to what extent, these affect the results. Do not forget to avoid that the receptor protein is dragged together with the ligand, for example, by immobilization of the protein’s center of mass (seeEq. 3). Do not immobilize each protein atom separately, as this would suppress conformational motions that are often crucial for binding-unbinding of a ligand.
- b.
Pull direction: In the AFM experiment, the ligand-receptor complex has sufficient time to adapt, by rotation, to the pulling direction set by the AFM. Because of Stokes’ friction, such adaptation takes quite a while at the MD time scale, and thus should be avoided by carefully selecting the correct pulling direction right away, typically by visual inspection. If in doubt, do a prepulling phase and monitor rigid-body adaptation motions of the system. Start the production runs only after the adaptation motions have stopped.
- c.
Stiffness of the cantilever spring: This choice depends on the questions to be addressed. If an AFM experiment is to be simulated, choose the effective spring coefficient of the experimental set-up. Note, that just using the spring coefficient of the cantilever may be inadequate, as also other elements may contribute to the effective spring that acts onto the ligand-receptor complex. Typically, this will be the entropic elasticity of the linker polymer(s), which can, for example, be described to sufficient accuracy by the worm like chain model (67). If emphasis is on the unbinding pathway, a soft spring (k < 0.01 N/m) can be adequate, as it leaves the ligand the freedom to move as soon as the force has become large enough to overcome the next energy barrier. If, rather the force profile along the unbinding pathway is of interest, a stiff spring ( >1 N/m) is required to obtain sufficient spatial resolution,
- d.
Pulling velocity or loading rate: The simple rule is: as slow as possible, given the computational limitations. In fact, the pulling simulations will necessarily, almost always be carried out too fast, that is to say at too short time scales. This is particularly true if AFM or optical tweezers experiments are simulated: whereas those are carried out at a milliseconds to seconds time scale, the simulations fall into the nanoseconds scale, i.e., are many orders of magnitude too fast. To alleviate this discrepancy, one should carry out several force probe simulations with different pulling velocities, and study the variation of the calculated rupture force in order to estimate if the experimental value falls onto this curve (see the previous Theory Section). To change this sometimes a bit risky (unstable) extrapolation into an interpolation, a good trick is to also consider the independently measured spontaneous unbinding rate (koff), which corresponds to vanishing rupture force.
- a.
- 6.
Considerable effort has to be spent in analyzing the obtained force probe simulations in order to provide a causal picture of the unbinding process and to relate the observed forces to interatomic interactions (see Note 5 ). Monitoring interactions like individual hydrogen bonds, water bridges, hydrophobic contacts, salt bridges, electrostatics, or van der Waals contacts, as well as water accessible surface or larger structural rearrangements during unbinding are typical ways to gain microscopic insight. Be careful when relating forces to interaction energies ( Note 6 ). No general rules on what observables should be analyzed can be given, though, and one often has to resort to biophysical and biological intuition. Not sufficiently appreciated particularly in the physics community is the fact that from pure inspection of movies (generated from the trajectories), much insight can be gained. In particular, movies often are an invaluable tool to develop ideas for subsequent quantitative analyses.
23.4 Notes
- 1.
Histidines and S-bonds. One cannot overemphasize the necessity to carefully check protonation states—particularly of histidines (possibly, if in doubt, with Poisson-Boltzmann calculations)—and the correct placement of disulfide bonds. Errors here can hide frustratingly long, even until after publication!
- 2.
Smoothing of force profiles. Particularly when using stiff pulling potentials Vcant the raw force profiles as obtained from Eq. 4 contain both the unbinding forces of interest, but also high frequency contributions arising from the thermal fluctuation of the atom subjected to Vcant. Typically, these force fluctuations are much larger than the underlying unbinding forces, and therefore have to be filtered, for example, with Gaussian filters (i.e., by a convolution of the raw force profile with a Gaussian function exp(-t2/22) of width). Good results have been obtained for determined from the resonance frequency of Vcant, i.e., =(m/k)1/2, where m is the (total) mass of the atom or atoms subjected to Vcant. Note that the choice of affects the maximum of the force profile, and hence the calculated dissociation force Fmax. This slight arbitrariness entails an uncertainty for Fmax.
- 3.
Equilibration and rms deviation. Take care when checking if a simulation system is sufficiently equilibrated (68). Besides the fact that it is often not quite clearly explained what sufficiently means here (strictly, it means that the results should not change if the system is equilibrated even much longer, but this does not help much in practice), there is generally no watertight proof that a system is indeed sufficiently equilibrated (see also refs.69–73). Thus, do not uncritically interpret any short-time plateau in the root mean square (rms) deviation from the start structure as system is equilibrated. If in doubt, extend the equilibration period to see whether the rms really stays constant. Today, equilibration phases shorter than 1 ns are generally considered insufficient for medium-sized proteins. Note, also that larger rms deviations can originate both from systematic drift away from the start structure (i.e., system is not equilibrated) or from large fluctuations (in which case the system can, but need not, be in equilibrium). Careful analysis, for example, by excluding parts of the protein that show large fluctuations from the rms calculation, is mandatory here. Also helpful are residue-based rms-plots that help to trace where a large rms value comes from. Note finally that an rms of, for example, 2.5 Å may be small for a large protein but unacceptably large for a small, rigid protein.
- 4.
Exerting torque. In recent experiments, utilizing magnetic beads (74) torque has been applied to single DNA molecules rather than the linear forces discussed in the main text. This process can be simulated by subjecting n selected atoms (with positions x i ) onto which the torque shall act, to a rotating potential Vtorque (13),
The rotation matrix Ω(t) has to be specified such that it moves the minima of Vtorque along concentric circles around the desired rotation axis, and the spring coefficient k defines the size of the allowed fluctuations around the equilibrium positions of the atoms. The torques Θ
i
,(t) exerted on the individual atoms are readily computed from the deflection of the actual atomic positions x
i
(t) from the respective minima,
and the total torque exerted on the molecule, is just the sum of the individual torques. This approach has, for example, been used in MD simulations to mimic the torque exerted by the proton-motive force onto the central stalk (γ-subunit) of F1-ATP synthase, and thereby to study the mechano-chemical energy conversion involved in ATP synthesis (13). - 5.
Anecdotal events. The huge computational amount of computer time that is required to carry out just a single force probe simulation at sufficiently low loading rate leads into the temptation to extract information on unbinding reaction pathways from just one trajectory. However, in this case one cannot tell which features of the observed unbinding event are reproducible, and which are purely anecdotal. Therefore, with limited computational resources consider carrying out, instead of one very extended simulation, several at somewhat higher pulling speed.
- 6.
Forces differ from free energies mainly in two respects. First, the size of forces and that of free energies are, generally, uncorrelated. Without further support, large forces do not necessarily mean large free energy differences (e.g., between bound and unbound state) and vice versa. If distributed over large distances, as in nonlocal or multiple interactions along an extended reaction path (e.g., protein folding), small forces can build up large free energy differences. On the other hand, a very stiff bond does not need much energy to generate large rupture forces. Second, free energy differences ΔG are thermodynamic state functions, and therefore depend only on start and end state but not on the chosen reaction pathway. This is not true for forces, which do depend on the particular reaction pathway. Therefore, forces yield information on pathways whereas binding free energies do not. One consequence is that when analyzing the effect of point mutations, forces are not additive: whereas the binding free energy changes as a result of a double mutation is to a good approximation simply the sum of the binding free energy change of the two single-point mutations, this is not true for the respective dissociation forces (75).
