Single-molecule pulling: phenomenology and interpretation
Single-molecule pulling techniques have emerged as versatile tools for probing the noncovalent forces holding together the secondary and tertiary structure of macromolecules. They also constitute a way to study at the single-molecule level processes that are familiar from our macroscopic thermodynamic experience. In this Chapter, we summarize the essential phenomenology that is typically observed during single-molecule pulling, provide a general statistical mechanical framework for the interpretation of the equilibrium force spectroscopy and illustrate how to simulate single-molecule pulling experiments using molecular dynamics.
Published in Nano and Cell Mechanics: Fundamentals and Frontiers, edited by H.D. Espinosa and G. Bao (Wiley, Microsystem and Nanotechnology Series, 2013), Chap. 14, pages 359–388.
- I Introduction
- II Force-extension behavior of single molecules
- III Single-molecule thermodynamics
- IV Modeling single-molecule pulling using molecular dynamics
- V Interpretation of pulling phenomenology
- VI Summary
In recent years a number of techniques have been developed that allow us to access the properties of single molecules (see, e.g., Refs. Ashkin (1997); Rief et al. (1997); Weiss (1999); Moerner and Orrit (1999); Strick et al. (2000); Nitzan and Ratner (2003); Ritort (2006); Camden et al. (2008)). These windows into the single-molecule world are teaching us with unprecedented detail how molecules behave, including how they move, their mechanical behavior, their conductivity, reactivity, optical properties, and other properties. While measurements made on bulk matter are averages over a whole ensemble of molecules, single-molecule measurements highlight the contributions of the constituent parts to the ensemble. A key feature of the properties of single molecules is that fluctuations in the observables are comparable with the average values. This contrasts with macroscopic behavior (where such fluctuations are usually negligible) and indeed fluctuation-induced emergent phenomena can arise that, simply put, cannot be observed in the macroscale.
In this Chapter we discuss basic aspects of single-molecule pulling experiments in which laser optical tweezers or atomic force microscopes (AFMs) are employed to exert mechanical forces on single molecules Rief et al. (1997); Liphardt et al. (2001); Oberhauser et al. (2001); Evans (2001); Bustamante et al. (2005), and the resulting stress/strain response is used as a probe of the molecular potential energy surface. The focus will be on how to simulate single-molecule pulling experiments using molecular dynamics, on what class of information can be extracted from the resulting force-extension isotherms and on how to interpret the basic phenomena typically observed upon pulling.
Figure 1 shows a schematic of a single-molecule pulling experiment using AFMs (an equivalent setup exists for experiments that employ laser optical tweezers). In it, one end of a molecule is attached to a surface, while the other end is attached to an AFM tip coupled to a cantilever. The distance between the surface and the cantilever is controlled while the molecular end-to-end distance is allowed to fluctuate. The pulling is performed by varying in some prescribed fashion. The force exerted during the process is determined by measuring the deflection of the cantilever from its equilibrium position , where is the cantilever stiffness. For a discussion of experimental aspects of single-molecule force spectroscopy see Ref. Neuman and Nagy (2008).
The experiments exert mechanical control over the molecular conformation while simultaneously making thermodynamic measurements of any mechanically-induced unfolding events. The stress maxima observed during pulling have been linked to the breaking Evans (2001); Freund (2009) of noncovalent interactions, thus providing insight into the secondary and tertiary structures of macromolecules. From the force-extension data it is possible to reconstruct the molecular free energy profile along the extension coordinate, quantifying in this way thermodynamic changes undergone by the molecule during folding.
The structure of this Chapter is as follows: Section II introduces the basic phenomenology observed during single molecule pulling through a discussion of representative experiments. Section III establishes the statistical mechanical theory behind single-molecule pulling experiments and describes a method for reconstructing the molecular free energy profile along the extension coordinate from the force-extension data. Section IV, in turn, introduces direct and indirect methods that can be used to model single-molecule pulling experiments using molecular dynamics. Section V focuses on the interpretation of single-molecule pulling experiments through an analysis of the basic structure of the molecular free energy profile. Section VI summarizes the main topics discussed in this Chapter.
Ii Force-extension behavior of single molecules
In this section, we consider a series of examples of single-molecule pulling experiments that illustrate the basic phenomenology that can be observed. Consider first the mechanically induced unfolding of a RNA hairpin Liphardt et al. (2001). The hairpin is attached to polystyrene beads via RNA/DNA hybrid handles and the pulling is performed in an optical tweezers arrangement. Figure 2 shows the experimentally determined force-extension behavior and the details of the hairpin structure. For fast pulling speeds hysteresis in the force-extension traces is evidenced by the fact that the traces observed upon pulling and subsequent contraction do not coincide. However, upon reduction of the pulling speed the traces obtained during extension and contraction do coincide, indicating that for all practical purposes the process occurs under reversible conditions. During pulling the force initially increases monotonically, then at pN it shows a drop, and then it increases again. In the process, the RNA hairpin undergoes a conformational transition from the hairpin structure to an extended coil. The drop in the force corresponds to the molecular unfolding event where the hydrogen bonding network holding the RNA hairpin together is broken. The local maximum in the force-extension isotherm is interpreted as the force required to break the secondary structure of the molecule.
The force-extension behavior of the RNA hairpin illustrates two basic phenomena that are commonly observed during single-molecule pulling. During pulling there is a region of mechanical instability where the average force decreases with the extension . This region occurs when the pulled molecule undergoes a conformational transition from a stable folded structure to an extended system. Note that around the unfolding region the molecule hops between the folded and the unfolded state. That is, around the region of mechanical instability the molecule plus cantilever system exhibits dynamical bistability where the molecule unfolds and refolds, leading to blinking in the force measurement between a stronger force when the molecule is folded and a weaker force regime when extended. In Sec. V we discuss the origin of these two phenomena.
Figure 3 shows the force-extension characteristics of a related system, a 6838 base-pair long DNA hairpin, also obtained in a optical tweezers arrangement Huguet et al. (2010). Even for this complex system it is possible to find experimentally accessible pulling speeds (10 nm/s) where reversible behavior is recovered in the mechanically induced melting process. Melting in this case refers to dehybridization of the base pairs, and we see that unlike Fig. 2, there is an extended unfolding region where the average force is in the 15-19 pN range while the length increases. This corresponds to sequential dehybridization (believed to occur a few base pairs at a time) rather than a concerted breaking of all the base pairs at the same time. The inset in the figure shows the raw data obtained in the experiment while the smooth curves shown are after further temporal averaging. The averaged data evidences the changes in the pulling force required to mechanically induce melting. It quantifies the strength of the type of hydrogen bonds encountered as the elongation proceeds.
The raw data evidences yet another fundamental aspect of this class of experiments: the substantial thermal fluctuations in the force measurements that are typically observed during pulling indicating that the system is not in the thermodynamic limit. This observation has the important consequence of leading to a non-equivalence between different statistical ensembles. In particular, it leads to a non-equivalence between the isometric ensemble (where is controlled and fluctuates) and the isotensional ensemble (where is controlled and fluctuates). Throughout we will focus on the isometric ensemble since is the most common way to perform these experiments. For recent discussions of the nonequivalence between ensembles see Kreuzer and Payne (2001); Süzen et al. (2009). See Kirmizialtin et al. (2005) for an interpretation of the pulling phenomenology in the isotensional ensemble.
As a third example, consider the reversible force-extension behavior of a synthetic polymer pulled in poor and good solvent conditions using an AFM Gunari et al. (2007). Figure 4 shows the force-extension isotherm of an individual polystyrene (PS) chain pulled in toluene (good solvent) and water (poor solvent). In toluene (Fig. 4a), the deformation of PS exhibits a worm-like chain Rubinstein and Colby (2003) behavior with a persistence length of (0.25 0.05) nm. In water (Fig. 4b) the behavior is qualitatively different. In this case, the force initially rises linearly up to pN, then exhibits a plateau force at pN, followed by a third regime where the force rises again. This behavior is well predicted by a model Halperin and Zhulina (1991) in which the initial linear restoring force is due to deformation of a collapsed spherical globule to an ellipsoid during which the polymer’s surface energy increases. In this model, the plateau is associated with a conformational transition involving the coexistence of the globule and a stretched chain of polymer or polymer globules. For large stretching, where coexistence is no longer possible, the force extension behavior of a Gaussian chain Rubinstein and Colby (2003) is recovered. Note that large-scale fluctuations in the force measurements are also evident in this system.
Consider last the forced unfolding of a recombinant polyprotein composed of eight repeats of the Ig27 domain of human titin obtained in an AFM setup Imparato et al. (2008). Figure 5 shows the force-extension profile when pulling at a constant speed under nonequilibrium conditions. As the polyprotein is extended the force exerted on the molecule increases until one of the Ig27 domains suddenly unfolds. Unfolding of a domain reduces the overall force, but with increasing extension another local subunit unfolds. This occurs until all the domains are unfolded and the protein-tip interaction is broken. This results in the sawtooth pattern shown in Fig. 5 where each peak in the force-extension profile marks the unfolding of one subunit in the polyprotein and the last peak signals the detachment of the protein from the tip.
Iii Single-molecule thermodynamics
In what follows we specify what we mean about thermodynamics in the context of single-molecule pulling experiments Bustamante et al. (2005) and the class of thermodynamic information that is accessible from the force-extension isotherms. The thermodynamic system consists of the molecule plus cantilever coupled to a thermal bath, typically the surrounding solvent. As such, the equilibrium state of the system is well described by the canonical ensemble with configurational partition function
where denotes the coordinates of the atoms of the molecule, is the inverse temperature, and
is the potential energy function of the molecule plus restraints. Here is the molecular potential energy plus potential restraints not varied during the experiment, and is the potential due to the cantilever at extension . To a good approximation, the quantity is given by a harmonic function
of stiffness . The instantaneous force exerted by the cantilever on the molecule is determined by measuring the deflection of the cantilever from its equilibrium position and it is given by
where the gradient is with respect to the coordinate. In writing Eqs. (3) and (4) the vector nature of , and has been obviated since these quantities are typically collinear in the single-molecule pulling setup. For notational simplicity, in Eq. (1) and throughout this manuscript we ignore the constant multiplicative factor that makes partition functions dimensionless since they are irrelevant for determining relative free energy changes.
iii.1 Free energy profile of the molecule plus cantilever
Because the state of the molecule plus cantilever under equilibrium conditions is well described by a canonical ensemble, the force exerted during the pulling yields the associated change in the Helmholtz free energy
for the composite molecule plus cantilever system. To see this, note that in the canonical ensemble the average force at a given extension is given by
Therefore, if the pulling is done under reversible conditions, such that the state of the system is well described by a canonical ensemble at each step of the pulling, the change in when pulling from to is determined by the reversible work exerted in the process
Remarkably, even when the pulling is performed under nonequilibrium conditions it is still possible to determine the free energy changes during pulling by means of the Jarzynski non-equilibrium work fluctuation relation Jarzynski (1997, 2008, 2011)
where the average in this expression is over nonequilibrium realizations of a given pulling protocol starting from a system initially prepared in the canonical ensemble. Eq. (8) relates the nonequilibrium work with the equilibrium free energy changes of the system and holds arbitrarily far away from equilibrium. Note that by using Jensen’s inequality (), Eq. (8) implies that the average nonequilibrium work is greater than or equal to the free energy change
which is the usual statement of the second law of thermodynamics
The above analysis associates a free energy change to the pulling process although the system is not in the thermodynamic limit. This association relies on the assumption that the equilibrium state of the system is well described by a classical canonical ensemble. This feature has been experimentally tested: In Liphardt et al. (2002) the authors experimentally demonstrated the validity of the Jarzynski equality [Eq. (8)]. The Jarzynski equality supposes that the initial equilibrium state of the system is given by a classical canonical partition function and hence the test of the Jarzynski equality also tests whether this is a good description of the equilibrium state of the system.
iii.2 Extracting the molecular potential of mean force
The properties that are measured during the pulling are those of the molecule plus cantilever. Thus, direct integration of the force-extension isotherms yields the free energy profile of the composite system [see Eq. (7)]. However, what one is really interested in are the properties of the molecule itself. Specifically, one is interested in the molecular Helmholtz free energy profile (also known as the potential of mean force PMF Kirkwood (1935)) along the extension coordinate . The PMF succinctly captures thermodynamic changes undergone by the molecule during folding and determines the basic phenomenology that is observed upon pulling Kirmizialtin et al. (2005); Franco et al. (2009).
A useful method to extract the molecular PMF from the force versus extension measurements is the Weighted Histogram Analysis Method (WHAM) Ferrenberg and Swendsen (1989); Kumar et al. (1992); Roux (1995); Frenkel and Smit (2002); Newman and Barkena (2007). The WHAM analysis combines all the force-extension data collected during the pulling to properly remove the bias due to the cantilever and extract the PMF. In fact, one can view the single-molecule pulling process as an experimental realization of the WHAM methodology using harmonic potential restraints. We now describe how to extract the PMF from equilibrium force-extension data using the WHAM. A schematic of the process is shown in Fig. 6.
The PMF is defined by Roux (1995)
where and are arbitrary constants and
is the configurational partition function of the molecule, and
In the context of the pulling experiments is the probability density that the molecular end-to-end distance function adopts the value in the unbiased (cantilever-free) ensemble. The quantity determines the PMF up to a constant and can be estimated from the force measurements, as we now describe.
At this point it is convenient to discretize the pulling process and suppose that extensions are visited during pulling. The potential of the system plus cantilever at extension is , where is the potential bias due to the cantilever. In the -th biased measurement knowledge of the force and length gives the molecular end-to-end distance. The probability density of observing the value at this given extension is given by
where is the configurational partition function for the system plus cantilever at the -th extension, and
where we have exploited the fact that . In practice this approach will not work because the range of values where and differ significantly from zero need to overlap Frenkel and Smit (2002). Nevertheless, the experiments do provide measurements at a wealth of values of that can be combined to estimate :
where the are some normalized () set of weights. In the limit of perfect sampling any set of weights should yield the same . In practice, for finite sampling it is convenient to employ a set that minimizes the variance in the estimate from the series of independent estimates of biased distributions. Such a set of weights is precisely provided by the WHAM prescription Ferrenberg and Swendsen (1989); Kumar et al. (1992); Frenkel and Smit (2002)
where is the statistical inefficiency Chodera et al. (2007) and the integrated autocorrelation time Ferrenberg and Swendsen (1989); Chodera et al. (2007). The difficulty in estimating the for each measurement generally leads to further supposing that the ’s are approximately constant and factor out of Eq. (18), so that
Neglecting the ’s from Eq. (19) does not imply that the resulting estimate of is incorrect; but simply that the weights selected do not precisely minimize the variance in the estimate. Experience with this method indicates that if the ’s do not differ by more than an order of magnitude their effect on is small Kumar et al. (1992); Hub et al. (2010).
The computation of the PMF then proceeds as follows. Suppose that force measurements are done for each extension . Knowledge of the force () and of gives the molecular end-to-end distance in each of these measurements. The numerator in Eq. (19) is then estimated by constructing a histogram with all the available data,
where is the bin size, and if and zero otherwise. Estimating the denominator in Eq. (19) requires knowledge of the ’s, the configurational partition functions of the system plus cantilever at all extensions considered. There are two ways to obtain these quantities. If reversible force-extension data is available, the most direct one is to employ the Helmholtz free energies for the system plus cantilever obtained through the thermodynamic integration in Eq. (7), as they determine the ’s up to a constant multiplicative factor. Alternatively, it is also possible to determine the ratio of the configurational partition functions between the -th biased system and its unbiased counterpart using :
Equations (19) and (21) can be solved iteratively. Starting from a guess for the , is estimated using Eq. (19) and normalized. The resulting is then used to obtain a new set of through Eq. (21), and the process is repeated until self-consistency. This latter approach is the usual procedure to solve the WHAM equations. For computational implementations of the WHAM methodology see Refs. Felberg et al. (2010); Grossfield (); Minh (); Hub et al. (2010).
iii.3 Estimating force-extension behavior from
The PMF is central to the interpretation of the phenomena observed upon pulling and in the estimation of the force-extension behavior when pulling with cantilevers of arbitrary stiffness Franco et al. (2009). Its utility relies on the fact that the configurational partition function of the molecule plus cantilever at extension [Eq. (1)] can be expressed as
where we have neglected the constant and -independent multiplicative factor that arises in the transformation since it is irrelevant for the present purposes. The above relation implies that the extension process can be viewed as thermal motion along a one-dimensional effective potential determined by the PMF and the bias due to the cantilever :
Further, by means of the PMF it is possible to estimate the force vs. extension characteristics for the composite system for any value of the force constant . This is because the average force exerted on the system at extension can be expressed as:
where, for convenience, we have introduced in the logarithm just to make the argument dimensionless. Since is a property of the isolated molecule it is independent of . Hence, Eq. (24) can be employed to estimate the - isotherms for arbitrary . Another quantity of interest is the probability density in the force measurements as a function of the extension , given by
as this quantity highlights the fluctuations in the force measurements during pulling.
Iv Modeling single-molecule pulling using molecular dynamics
iv.1 Basic computational setup
Computationally, the forced unfolding of single molecules can be studied using “pulling” or “steered” molecular dynamics (MD) simulations Izrailev et al. (1998); Isralewitz et al. (2001); Heymann and Grubmüller (1999). The general setup of this type of simulations is schematically described in Fig. 7. The stretching computation begins by attaching one end of the molecule to a rigid isotropic harmonic potential that mimics the molecular attachment to the surface. Simultaneously, the opposite molecular end is connected to a dummy atom via a virtual harmonic spring. The position of the dummy atom is the simulation analogue of the cantilever position , and is controlled throughout. In turn, the varying deflection of the virtual harmonic spring measures the force exerted during the pulling. The stretching is caused by moving the dummy atom away from the molecule at a constant speed. The pulling direction is defined by the vector connecting the two terminal atoms of the complex. Since cantilever potentials are typically stiff in the direction perpendicular to the pulling, in the simulations the terminal atom that is being pulled is typically forced to move along the pulling direction by introducing appropriate additional restraining potentials.
In the simulations, the molecule plus any surrounding solvent are assumed to be weakly coupled to a heat bath and the effect of the heat bath is modeled through a thermostat. Thermostats that are useful in simulating the force spectroscopy include the Nosé-Hoover chain Nosé (1984); Hoover (1985); Martyna et al. (1992); Frenkel and Smit (2002) and the Langevin thermostat Thijssen (2007). In turn, velocity rescaling thermostats, like the Berendsen thermostat Berendsen et al. (1984), can be problematic because they can lead to a spurious violation of the second law during the pulling because the ensemble generated by them does not satisfy the equipartition theorem Harvey et al. (1998). The thermal dynamics of the molecule, cantilever and any surrounding solvent are then followed in a NVT or NPT ensemble using MD. Freely available MD codes such as NAMD Phillips et al. (2005) and GROMACS Hess et al. (2008) have pulling routines implemented in them. We have, in addition, developed Franco et al. (2009) a pulling routine for TINKER Ponder (2004) and created a graphical interface for it called MOLpull Felberg et al. (2010). MOLpull is accessible through the NanoHub (http://nanohub.org/) and can be used to perform the pulling and reconstruct the PMF along the extension coordinate of arbitrary molecules.
iv.2 Modeling strategies
One of the challenges in simulating pulling experiments using MD is to bridge the large disparity between the pulling speeds that are employed experimentally and those that can be accessed computationally. While typical experiments often use m/s, current computational capabilities require pulling speeds that are several (5 to 8) orders of magnitude faster.
One possible strategy that can be used to overcome this difficulty is to strive for pulling speeds that are slow enough that reversible behavior is recovered. Under such conditions the results become independent of , and the simulations comparable to experimental findings. The way to determine if the pulling is done under quasistatic conditions is by making sure that the force-extension curves during extension and subsequent contraction essentially coincide, such that the net work exerted during a pulling/retraction cycle is negligible. The advantage of this approach is that there is a systematic way to test if quasistatic behavior is recovered by estimating the amount of work dissipated in a pulling/retraction cycle. Its limitation is that the pulling speeds that are required to obtain reversible behavior may be computationally unfeasible for pulling large macromolecules. In fact, the examples that exist in the literature where reversible behavior has been computationally recovered are for relatively small molecular systems Franco et al. (2009, 2011a); Schäfer et al. (2007) using pulling speeds between nm/ps. When pulling large macromolecules, present computational limitations require using pulling speeds of nm/ps Lee et al. (2010) which is far too fast for recovering reversible behavior even in simple systems.
A second possible strategy is to pull under irreversible conditions. A single nonequilibrium trajectory can recover some qualitative structural features encountered during the pulling Lee et al. (2010); Heymann and Grubmüller (1999), but it is not sufficient to reconstruct the molecular potential of mean force. In order to reconstruct the PMF it is necessary to pull many times and use a nonequilibrium work fluctuation theorem Jarzynski (1997); Crooks (1999); Jarzynski (2008, 2011) to reconstruct the free energy profile Hummer and Szabo (2001); Park and Schulten (2004); Hummer and Szabo (2005); Harris et al. (2007); Imparato et al. (2008); Minh and Adib (2008); Pohorille et al. (2010); Hummer and Szabo (2010). The advantage of this strategy is that it can be used when equilibrium data is difficult to obtain, for example during the forced unfolding of a complex protein like the one in Fig. 5. Note, however, that this strategy requires adequately estimating the average of an exponential quantity (Eq. (8)) that often has poor convergence properties Jarzynski (2006); Pohorille et al. (2010). For a computational implementation of some nonequilibrium methods to reconstruct the PMF see Minh ().
As a third indirect alternative it is possible to reconstruct the molecular potential of mean force without pulling continuously but rather by sampling the fluctuating force at selected points along the extension coordinate. For each the system is allowed to thermally equilibrate and the force is measured for a given amount of simulation time. The data obtained is used to reconstruct the PMF using WHAM as described in Sec. III, and then the PMF is employed to estimate the force-extension behavior of the molecule plus cantilever. While the free energy reconstruction using continuous equilibrium force-extension data uses a data set with a large (essentially continuous) set of values with very few data points for the force at each , this strategy works in a different limit where one has few extensions but extensive sampling of the force at each extension. The advantage of this procedure is that it can be used to reconstruct the PMF of molecules that are computationally challenging to pull under equilibrium or near equilibrium conditions. The limitation is that there is no entirely satisfactory way to judge when the system has achieved a state of thermal equilibrium in the MD trajectory. This contrasts with the direct strategy involving explicit pulling in which it is easy to judge when quasistatic behavior has been recovered by measuring the amount of dissipative work in extension/contraction cycles. As in the direct equilibrium pulling, the sampling for each has to be sufficient to capture all relevant molecular events that are consistent with the extension.
Another important aspect that has to be taken into consideration when modeling single-molecule pulling is the cantilever stiffness employed. Experimentally accessible AFM cantilevers stiffnesses are in the pN/Å range Neuman and Nagy (2008). While experiments tend to favor soft cantilever potentials because they provide high resolution in the force measurements, in the simulations stiff cantilevers are preferable. This is because when employing stiffer cantilever potentials: (i) a shorter range of values is required to explore all the molecular extension coordinates thus reducing the simulation time; (ii) the possible values of that are accessible for the molecule for a given are further restricted making it faster to sample all accessible conformational space. As a consequence of (ii), the amount of dissipative work in a pulling/retraction cycle decreases with increasing cantilever stiffness for fixed pulling speed Marsili and Procacci (2010). In our experience with molecular pulling, a rule of thumb that we frequently use is that a cantilever stiffness of pN/Å can be regarded as soft and one with stiffness of pN/Å as stiff. These numbers, however, will clearly depend on the particular system that is being pulled.
Direct modeling: Oligorotaxanes
As an example of an explicit simulation of a single-molecule pulling experiment consider the forced unfolding of a rotaxane Franco et al. (2009, 2011b); Basu et al. (2011). Figure 8 shows the force-extension curves during the extension (blue) and subsequent contraction (red) for two different pulling speeds and a schematic of the structure of the molecule that is being pulled. The oligorotaxane
The force-extension behavior predicted by the MD simulations for the oligorotaxane is qualitatively similar to that experimentally observed for the RNA hairpin (recall Fig. 2). Specifically, there is a region of mechanical instability where the force drops with increasing extension (). Further, if one fixes around the region of mechanical instability ( Å in this case) and allows the system to evolve, the molecule undergoes conformational transitions between a folded globular state and an extended coil, see Fig. 9. The right panel in Fig. 9 shows the probability density of the distribution of values obtained from a 20 ns trajectory. The system exhibits a clear dynamical bistability along the coordinate. As in the RNA case, the bistability along the end-to-end distance leads to blinking in the force measurements from a high force to a low force regime during the pulling.
The molecular PMF of the rotaxane reconstructed from the force-extension data is shown in Fig. 10. The thermodynamic native state of the oligorotaxane, i.e. the minimum in , corresponds to a globular folded structure with an average end-to-end distance of 7 Å. The PMF consists of a convex region (a region of positive curvature) for Å when the molecule is folded, another convex region for Å when the molecule is extended and a region of concavity (where ) for Å where the conformational transition occurs. The region of concavity along the PMF arises because the molecule is inherently bistable along some natural unfolding pathway. As discussed in Sec. V, this characteristic structure of the PMF is responsible for many of the interesting features observed during pulling.
Indirect modeling: DNA dimer
As an example of indirect modeling of the force spectroscopy, we now consider the forced extension of the DNA dimer shown in Fig. 11 along the O5’-O3’ coordinate in explicit water McCullagh et al. (2011). The dimer is composed of two guanine-cytosine base pairs capped by a trans-azobenzene linker. Computationally, it is described by the CHARMM27 force field and the total system is composed of the 178 DNA hairpin atoms, 4 sodium ions, and 1731 TIP3P water molecules in a ÅÅÅ box in the NVT ensemble at 300K (see Ref. McCullagh et al. (2011) for details). Even for this relatively simple system it is challenging to computationally obtain equilibrium force-extension data by pulling directly because the pulling has to be slower than any characteristic timescale of the system. An indirect approach is thus more convenient.
The sampling along the extension coordinate required to reconstruct the potential of mean force proceeded as follows: The distance, , between the surface and the cantilever was fixed at several different values along the extension coordinate. The system was first allowed to equilibrate for 4 ns. Subsequently, the dynamics was followed for 8 ns and the end-to-end distance recorded every 1 ps. The cantilever was taken to have a stiffness of pN/Å. In total, 155 extensions were simulated. Values of ranged from Å to Å while ranged from Å to Å. Total analyzed simulation time was s.
Figure 12 shows some of the sampled [defined in Eq. (14)] that are used as input for the PMF reconstruction. The resulting Helmholtz free energy profile is shown in Fig. 13. It consists of a convex region for Å representing the folded state, a mostly concave region for 17 Å 40 Å where the molecule unfolds, followed by another region of convexity that corresponds to the unfolded state. The unfolding region exhibits some convex intervals signaling marginally stable intermediates encountered during the unfolding.
Using the PMF one can then estimate the force-extension characteristics of the system when pulling with cantilevers of arbitrary stiffness as described in Sec. III.3. Figure 14 shows the predicted average and probability density in the force measurements when the DNA dimer is pulled using a cantilever stiffness of 1.1 pN/Å. Note the degree of detail of the pulling process revealed by the probability density.
V Interpretation of pulling phenomenology
A striking feature of single-molecule pulling is that the elastic behaviors of vastly different molecules often have qualitatively similar features. In Secs. II and IV we discussed two of these features: (i) mechanically unstable regions where and (ii) regions of dynamical bistability where blinking in the force measurements for fixed is observed. In this section, we discuss the origin of both of these effects. Emphasis will be placed on the basic requirements on the molecule and the pulling device for the emergence of the phenomena.
v.1 Basic structure of the molecular potential of mean force
Recall the basic structure of the PMF for the oligorotaxanes (Fig. 10): It consists of a convex region () where the molecule is folded, another region of convexity when the molecule is extended and a region of concavity () where the molecule unfolds. This basic structure in the PMF is not unique to the rotaxanes but it is actually a common feature of the extension behavior of single molecules of sufficient structural complexity. Regions of convexity mark mechanically stable conformations while regions of concavity mark molecular unfolding events. We had already encountered this very same basic structure when studying the elastic properties of the DNA dimer shown in Fig. 13. Perhaps a more spectacular example of this structure is provided by the polyprotein composed of eight repeats of the Ig27 domain of human titin that we introduced in Fig. 5. The PMF reconstructed from experimental force-extension data is shown in Fig. 15. In this case also, mechanically stable conformations are described by regions of convexity in the PMF, while regions of concavity mark unfolding events. The reason why in this case there are 8 regions of concavity along the extension pathway is that there are 8 unfolding events, each marking the unfolding of one of the Ig27 domains in the polyprotein. Below we discuss the implications of these regions of concavity along the PMF.
v.2 Mechanical instability
A first consequence of the regions of concavity along the molecular PMF is that they lead to mechanically unstable regions in the - isotherms. To see this consider the configurational partition function of the molecule plus cantilever at extension [Eq. (22)] for large . In this regime, most contributions to the integral will come from the region where . Consequently, can be expanded around this point to give:
Using Eq. (24), the derivative of the force can be obtained from this approximation to the partition function. To zeroth order in it is given by:
That is, an unstable region in the force vs. extension where requires a region of concavity in the PMF.
Note, however, that since the properties that are measured during pulling are those of the molecule plus cantilever, the observed mechanical stability properties will depend on the cantilever stiffness employed. In fact, when employing very soft cantilevers no mechanically unstable region can be observed independent of the molecule that is being pulled, provided that no bond-breaking occurs during the pulling. To see this, consider the soft-spring approximation of the configurational partition function of the system plus cantilever [Eq. (22)]:
where the notation stands for the unbiased (cantilever-free) average of . In writing Eq. (29) we have supposed that in the region of relevant (in which the integrand is non-negligible) and, hence, that the cantilever potential is well approximated by . This approximation is valid provided that that no bond breaking is induced during pulling and permits the introduction of a cumulant expansion Kubo (1962) in the configurational partition function. Specifically, the average can be expressed as:
It then follows that to lowest order in ,
That is, no unstable region in the - curve can arise in the soft-spring limit, irrespective of the specific form of .
Figure 16 illustrates the dependence of the region of mechanical instability on the cantilever stiffness using the rotaxane case discussed in Sec. IV.3.1 as an example. The figure shows the -dependence of the local maximum and minimum in the average force measurements that enclose the region of mechanical instability. The critical forces and show a smooth and strong dependence on the cantilever spring constant. For large the system persistently shows a region of instability in the isotherms. However, as the cantilever spring is made softer, and approach each other and for small the mechanical instability in the - isotherms is no longer present.
Another interesting aspect of single-molecule pulling is that there is a close relationship between the mechanical stability properties and the fluctuations in the force measurements. To see this, note that the average slope of the - curves can be expressed in terms of the fluctuations in the force:
where we have employed Eq. (24). The sign of determines the mechanical stability during the extension and hence Eq. (33) relates the thermal fluctuations in the force measurements with the stability properties of the - curves. Specifically, for the stable regions for which the force fluctuations satisfy
|In turn, in the unstable regions the force fluctuations are larger than in the stable regions and satisfy the inequality|
The fact that the fluctuations are larger around the region of mechanical bistability is noticeable in some of the force-extension curves that we have discussed here (see Figs. 8 and 14). Figure 17 illustrates these general observations in the specific case of the pulling of rotaxane. As shown, for and (where pN/Å) there is a region where the force fluctuations become larger than and, consequently, unstable behavior in - develops. For or less the fluctuations in the force are never large enough to satisfy Eq. (34b) and no critical points in the average - curve develop, as can be confirmed in Fig. 16.
v.3 Dynamical bistability
The observation of force measurements that blink between a high-force and a low force regime for selected extensions requires the composite molecule plus cantilever system to be bistable along the end-to-end distance for some . That is, the effective potential [Eq. (23)] must have a double minimum. Since the molecular PMF is not usually bistable along the bistability must be introduced by the cantilever potential. To see how this bistability arises, consider the effective potential for the rotaxane for selected shown in Figure 18. For in the mechanically stable regions of the force-extension isotherms ( Å and Å) the effective potential exhibits a single minimum along the coordinate. However, when Å a bistability in the potential develops. At this extension the cantilever potential turns the region of concavity in into a barrier between two minimum. The secondary minimum is the cause for the blinking in the force measurements observed at this (cf. Fig. 9).
What are the minimum requirements for the emergence of bistability along ? A necessary condition is that the effective potential is concave for some region along , that is:
For Eq. (35) to be satisfied it is required that both: (i) the PMF of the isolated molecule has a region of concavity where , and (ii) the cantilever employed is sufficiently soft such that
for some . Equation (36) imposes an upper bound on for bistability to be observable. If is very stiff the inequality would be violated for all and bistability would not be manifest. Note, however, that there is no lower bound for that prevents bistability along .
Since the region of concavity in the PMF is a common feature of molecules that have stable folded conformations, this bistability for selected extensions is a common feature of the pulling. We have already encountered this phenomenon during the pulling of the oligorotaxane and the RNA hairpin (Fig. 2). The DNA dimer (Fig. 14) and the polyprotein in Fig. 5 are also expected to show this behavior since their PMF has regions of concavity along . The bistability has also been predicted to be measurable when pulling -stacked molecules Franco et al. (2011a); Kim et al. (2009). For a discussion of the bistability in the isotensional ensemble, see Ref. Kirmizialtin et al. (2005). Figure 19 shows a striking recent experimental demonstration of the bistability obtained during the pulling of a leucine zipper Gebhardt et al. (2010).
The power of single-molecule pulling techniques is that they allow for mechanical control over the molecular conformation while simultaneously making thermodynamic measurements of any mechanically induced unfolding events. The force-extension data often exhibits mechanically unstable regions where the force decreases with increasing extension and bistable regions where, for fixed extension, blinks between a high-force and a low-force regime are observed.
From the force-extension isotherms it is possible to reconstruct the PMF along the extension coordinate using the WHAM. In fact, the pulling process can be seen as an experimental realization of the WHAM methodology. The PMF summarizes the changes in the free energy during folding and determines the elastic properties of the molecule. Quite remarkably, the PMF’s of widely different molecules often shares the same basic structure. It consists of regions of convexity that represent mechanically stable molecular conformations interspersed by regions of concavity that signal molecular unfolding events.
The starting point in the interpretation of single-molecule pulling experiments is Eq. (22) which shows that the pulling process can be viewed as thermal motion along a one-dimensional potential that is determined by the PMF and the cantilever potential . Several important conclusions follow from this simple observation. The first and most obvious one is that the properties that are measured during pulling are those of the molecule plus cantilever and hence that the basic phenomenon observed during pulling depend on the cantilever stiffness used. It also follows that in order for the mechanical instability and dynamical bistability to be observable, the molecular PMF has to have a region of concavity. However, while the dynamical bistability survives for soft cantilever and decays for rigid cantilevers [see Eq. (36)], the opposite is true for the mechanical instability [recall Eq. (32)]. The mechanical stability properties during pulling were also shown to be intimately related to the fluctuations in the force measurements [see Eq. (34)].
One basic challenge in modeling the force spectroscopy using MD is to bridge the several orders of magnitude gap that exists between the pulling speeds used experimentally and those that are computationally feasible. We detailed two possible strategies that can be employed to overcome this difficulty: either strive for pulling speeds that are slow enough such that reversible behavior is recovered or use an indirect approach in which the PMF is first reconstructed and the force-extension behavior is then estimated from the PMF. Other strategies based on nonequilibrium pulling are also possible. From the simulation of the pulling one can gain insights into the conformations encountered during molecular unfolding and into structure-function relations responsible for the elastic behavior of single molecules.
Currently, the accurate reconstruction of the PMF of molecular systems using equilibrium sampling is a serious computational challenge except for relatively modest systems. The recent discovery of nonequilibrium work fluctuation relations has led to a surge of activity in seeking nonequilibrium methods to reconstruct the PMF that can be used in situations where proper equilibrium sampling is unfeasible. The full consequences of this venue of research is still to be determined.
Acknowledgements.This work was supported by the Non-equilibrium Energy Research Center (NERC) which is an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences under Award Number DE-SC0000989. The authors thank Dr. Martin McCullagh for useful discussions.
- To be precise, we actually deal with pseudorotaxanes that lack the steric blocking units at the chain termini that force the cyclophanes to remain attached to the oligomeric chain. For simplicity, we will use the term rotaxane.
- A. Ashkin, Proc. Natl. Acad. Sci. U. S. A. 94, 4853 (1997).
- M. Rief, M. Gautel, F. Oesterhelt, J. M. Fernandez, and H. E. Gaub, Science 276, 1109 (1997).
- S. Weiss, Science 283, 1676 (1999).
- W. E. Moerner and M. Orrit, Science 283, 1670 (1999).
- T. R. Strick, V. Croquette, and D. Bensimon, Nature 404, 901 (2000).
- A. Nitzan and M. A. Ratner, Science 300, 1384 (2003).
- F. Ritort, J. Phys.: Condens. Matter 18, R531 (2006).
- J. P. Camden, J. A. Dieringer, Y. Wang, D. J. Masiello, L. D. Marks, G. C. Schatz, and R. P. Van Duyne, J. Am. Chem. Soc. 130, 12616 (2008).
- J. Liphardt, B. Onoa, S. B. Smith, I. Tinoco, and C. Bustamante, Science 292, 733 (2001).
- A. F. Oberhauser, P. K. Hansma, M. Carrion-Vazquez, and J. M. Fernandez, Proc. Natl. Acad. Sci. U. S. A. 98, 468 (2001).
- E. Evans, Annu. Rev. Biophys. Biomol. Struct. 30, 105 (2001).
- C. Bustamante, J. Liphardt, and F. Ritort, Phys. Today 7, 43 (2005).
- K. C. Neuman and A. Nagy, Nat. Meth. 5, 491 (2008).
- L. B. Freund, Proc. Natl. Acad. Sci. U. S. A. 106, 8818 (2009).
- J. M. Huguet, C. V. Bizarro, N. Forns, S. B. Smith, C. Bustamante, and F. Ritort, Proc. Natl. Acad. Sci. U. S. A. 107, 15431 (2010).
- H. J. Kreuzer and S. H. Payne, Phys. Rev. E 63, 021906 (2001).
- M. Süzen, M. Sega, and C. Holm, Phys. Rev. E 79, 051118 (2009).
- S. Kirmizialtin, L. Huang, and D. E. Makarov, J. Chem. Phys. 122, 234915 (2005).
- N. Gunari, A. C. Balazs, and G. C. Walker, J. Am. Chem. Soc. 129, 10046 (2007).
- M. Rubinstein and R. H. Colby, Polymer Physics (Oxford University Press, New York, 2003).
- A. Halperin and E. B. Zhulina, Europhys. Lett. 15, 417 (1991).
- A. Imparato, F. Sbrana, and M. Vassalli, Europhys. Lett. 82, 58006 (2008).
- C. Jarzynski, Phys. Rev. Lett. 78, 2690 (1997).
- C. Jarzynski, Eur. Phys. J. B 64, 331 (2008).
- C. Jarzynski, Annu. Rev. Condens. Matter Phys. 2, 329 (2011).
- J. Liphardt, S. Dumont, S. B. Smith, I. Tinoco, and C. Bustamante, Science 296, 1832 (2002).
- J. G. Kirkwood, J. Chem. Phys. 3, 300 (1935).
- I. Franco, G. C. Schatz, and M. A. Ratner, J. Chem. Phys. 131, 124902 (2009).
- A. M. Ferrenberg and R. H. Swendsen, Phys. Rev. Lett. 63, 1195 (1989).
- S. Kumar, D. Bouzida, R. H. Swendsen, P. A. Kollman, and J. M. Rosenberg, J. Comp. Chem. 13, 1011 (1992).
- B. Roux, Comput. Phys. Commun. 91, 275 (1995).
- D. Frenkel and B. Smit, Understanding Molecular Simulation (Academic Press, 2002), 2nd ed.
- M. E. J. Newman and G. T. Barkena, Monte Carlo Methods in Statistical Physics (Oxford University Press, New York, 2007).
- J. D. Chodera, W. C. Swope, J. W. Pitera, C. Seok, and K. A. Dill, J. Chem. Theory Comput. 3, 26 (2007).
- J. S. Hub, B. L. de Groot, and D. van der Spoel, J. Chem. Theory Comput. 6, 3713 (2010).
- L. Felberg, I. Franco, M. McCullagh, M. A. Ratner, G. C. Schatz, and M. Carignano, Molpull: a tool for molecular free energy reconstruction along a pulling coordinate, DOI:10254/nanohub-r9583.2 (2010).
- A. Grossfield, The weighted histogram analysis method, http://membrane.urmc.rochester.edu/content/wham.
- D. D. L. Minh, Ferbe free energy reconstruction from biased experiments, https://simtk.org/home/ferbe.
- S. Izrailev, S. Stepaniants, B. Isralewitz, D. Kosztin, H. Lu, F. Molnar, W. Wriggers, and K. Schulten, in Computational Molecular Dynamics: Challenges, Methods, Ideas, edited by P. Deuflhard, J. Hermans, B. Leimkuhler, A. E. Marks, S. Reich, and R. D. Skeel (Springer-Verlag, 1998), vol. 4 of Lecture Notes in Computational Science and Engineering, pp. 39–65.
- B. Isralewitz, M. Gao, and K. Schulten, Curr. Opin. Struct. Biol. 11, 224 (2001).
- B. Heymann and H. Grubmüller, Chem. Phys. Lett. 305, 202 (1999).
- S. Nosé, Mol. Phys. 52, 255 (1984).
- W. G. Hoover, Phys. Rev. A 31, 1695 (1985).
- G. J. Martyna, M. L. Klein, and M. Tuckerman, J. Chem. Phys. 97, 2635 (1992).
- J. M. Thijssen, Computational Physics (Cambridge University Press, New York, 2007), 2nd ed.
- H. J. C. Berendsen, J. P. M. Postma, W. F. Vangunsteren, A. Dinola, and J. R. Haak, J. Chem. Phys. 81, 3684 (1984).
- S. C. Harvey, R. K. Z. Tan, and T. E. Cheatham, J. Comp. Chem. 19, 726 (1998).
- J. C. Phillips, R. Braun, W. Wang, J. Gumbart, E. Tajkhorshid, E. Villa, C. Chipot, R. D. Skeel, L. Kalé, and K. Schulten, J. Comp. Chem. 26, 1781 (2005).
- B. Hess, C. Kutzner, D. van der Spoel, and E. Lindahl, J. Chem. Theory Comput. 4, 435 (2008).
- J. Ponder, Tinker: Software tools for molecular design 4.2, http://dasher.wustl.edu/tinker/ (2004).
- I. Franco, C. B. George, G. C. Solomon, G. C. Schatz, and M. A. Ratner, J. Am. Chem. Soc. 133, 2242 (2011a).
- L. V. Schäfer, E. M. Müller, H. E. Gaub, and H. Grubmüller, Angew. Chem. 119, 2282 (2007).
- E. H. Lee, J. Hsin, E. von Castelmur, O. Mayans, and K. Schulten, Biophys. J. 98, 1085 (2010).
- G. E. Crooks, Phys. Rev. E 60, 2721 (1999).
- G. Hummer and A. Szabo, P. Natl. Acad. Sci. U.S.A. 98, 3658 (2001).
- S. Park and K. Schulten, J. Chem. Phys. 120, 5946 (2004).
- G. Hummer and A. Szabo, Acc. Chem. Res. 38, 504 (2005).
- N. C. Harris, Y. S., and C.-H. Kiang, Phys. Rev. Lett. 99, 068101 (2007).
- D. D. L. Minh and A. B. Adib, Phys. Rev. Lett. 100, 180602 (2008).
- A. Pohorille, C. Jarzynski, and C. Chipot, J. Phys. Chem. B 114, 10235 (2010).
- G. Hummer and A. Szabo, Proc. Natl. Acad. Sci. U. S. A. 107, 21441 (2010).
- C. Jarzynski, Phys. Rev. E 73, 046105 (2006).
- S. Marsili and P. Procacci, J. Phys. Chem. B 114, 2509 (2010).
- I. Franco, M. A. Ratner, and G. C. Schatz, J. Phys. Chem. B 115, 2477 (2011b).
- S. Basu, A. Coskun, D. C. Friedman, M. A. Olson, D. Benítez, E. Tkatchouk, G. Barin, J. Yang, A. C. Fahrenbach, W. A. Goddard, et al., Chem. Eur. J. 17, 2107 (2011).
- N. L. Allinger, Y. H. Yuh, and J. H. Lii, J. Am. Chem. Soc. 111, 8551 (1989).
- M. McCullagh, I. Franco, M. A. Ratner, and G. C. Schatz, J. Am. Chem. Soc. 133, 3452 (2011).
- R. Kubo, J. Phys. Soc. Jpn. 17, 1100 (1962).
- J. S. Kim, Y. J. Jung, J. W. Park, A. D. Shaller, W. Wan, and A. D. Q. Li, Adv. Mater. 21, 786 (2009).
- J. C. M. Gebhardt, T. Bornschlögl, and M. Rief, Proc. Natl. Acad. Sci. U. S. A. 107, 2013 (2010).