How to probe the molecular energy landscape: “force spectroscopy” and application of work-fluctuations-theorems Diego Frezzato, September 2013 PART I Force and torque spectroscopies Technical means to handle (mechanically) single molecules have been developed in the last 20 years in biological contexts (manipulation of myosins, RNA polymerase, kinesins, DNA, etc.) Force spectroscopies: unfolding/refolding along an internal “reaction coordinate” defined by the pair of points (residues) of force application. Experimental devices to apply forces/torques - Laser-trap deflection (tweezers) - Cantilevers - Glass needles - Magnetic traps - Flow fields Example with optical traps bead optical trap (fixed) molecule handle handle pipette (moving) F ∆ Measurable istantaneous displacement ∆ → istantaneous force (known by calibration) bead Laser-trap deflection (laser tweezers) From the middle of the ‘80s: A. Ashkin et al, Opt. Lett. 11, 288 (1986) A. Ashkin, PNAS 94, 4853 (1997) R. Simmons et al, Biophys. J. 70, 1813 (2006) • Use of “gentle” laser radiation in the IR range (to prevent damage of chemicals). • The radiation produces two forces which act on a neutral polarisable particle (typically, a colloidal polystyrene bead): - pressure force along the propagation direction; - gradient force orthogonal to the propagation direction, originated by inhomogeneity of the radiation intensity (it decreases moving away from the beam axis on radial directions…) • A double-beam setup with two beams in opposition minimize the drift due to the pressure force. • Highly-polarizable particles (i.e. with high refractive index) prevail on solvent molecules to occupy the laser spot location. Radial Gaussian profile of the radiation intensity ⇒ radial force of elastic type (Hooke-like): F =κ ∆ The actual displacement of the bead from the spot center is monitored by using a second radiation conveyed to pass through the trap center. The intensity of the emergent radiation is measured as electric current by a photodiode. The current i grows (roughly) linearly with the distance ∆ of the bead from the trap center: ∆ (t ) = βr i (t ) Stiffness constant κ and response factor βr can be determined by calibration; notice that these parameters depend also on the used bead. Calibration is made via spectral (Fourier) analysis of the fluctuations of the current intentisy at thermal equilibrium: νc Figure adapted from http://www.nbi.dk/~tweezer/introduction.htm κ = 2πγ ν c where ν c is the frequency at the turning point of the power spectrum, and γ is the friction coefficient of the bead in the viscous medium (evaluated via Stokes law…) βr = k BT i2 κ Double-beam laser tweezers Figures taken from the talk of Ciro Cecconi at Workshop BIOSTRUCT09, Firenze 2009 www.biostruct09.fi.isc.cnr.it Ultrastable Atomic Force Microscope combined with cantilevers Molecular extension can be hold fixed for hundreds of seconds. G. M. King et al, Nano Lett. 9, 1451 (2009) Also in this case, the applied force has a Hook-like form with stiffness constant to be determined by “thermal tune” (similar to that for optical tweezers); see http://www.veeco.com/pdfs/appnotes/an90_thermaltune_11115_121.pdf Glass needles A. Kishino, T. Yanagida, Nature 334, 74 (1988). Micrometer-sized needles (50-500 µm length, 0.1-1 µm width); higher sensitivity with respect to cantilevers, suitable for biological systems. One end is fixed to a glass slab, the other is tethered to the molecule. Bending of the needle corresponds to an applied force. Not yet a standard tool. Magnetic traps Uses two micromagnets to drive a micrometric magnetic bead (radius < 3 µm) tethered to the molecule via a molecular spacer. Suitable to generate torques by rotating the magnetic field (see in the following…) Flow fields Exploitation of laminar flow to exert an impulse to a bead through controlled transfer of momentum. Table from C. Bustamante et al, Nature Reviews 1, 130 (2000); see references therein. Given the characteristic “stiffness” constant κ : uncertainty on bead localization: δ∆ 2 = k BT κ 2 uncertainty on applied force: δ F = κ k BT AFM Cantilevers are denoted as “rigid devices” working in high-force regimes. Laser tweezers work in the low-force regime (0.1-100 pN); most used devices. Operation protocols with laser tweezers At the zero-time, the molecule is “picked” at thermal equilibrium; then a tranformation protocol is applied. The transformation is repeated several times under same conditions to perform a statistical analysis. × Constant velocity (of the pipette) v = const. x × feedback adjustment (“istantaneous”) of the force to keep x constant ∆ F Extension clamp The bead-bead separation is kept fixed (while the force fluctuates) Y. Suzuki, O. K. Dubko, Phys. Rev. Lett. 110, 158105 (2013). × feedback adjustment (“istantaneous”) to keep the force constant ∆ F × feedback adjustment (“istantaneous”) to fix the force to the target value at time t ∆ F(t) Force clamp The force is kept constant (bead-bead distance fluctuates) K. Visscher et al, Nature 400, 184 (1999) A. Oberhauser et al, PNAS 98, 468 (2001) W. J. Greenleaf et al, Phys. Rev. Lett. 95, 208102 (2005) Force ramp The force is changed at a prescribed rate. M. Rief et al, Science 276, 1109 (1997) E. Evans, K. Ritchie, Biophys. J. 72, 1541 (1997) Theoretical picture of force ramp/clamp experiments (a perspective) Q1 QM molecule handle 1 Q2 handle 2 1 × 2 ∆ x (bead-bead) F feedback adjustment (“istantaneous”) to set the force QM : internal (configurational) coordinates of the molecule Q1 : internal coordinates of the handle 1 Q2 : internal coordinates of the handle 2 Q = (Q M , Q1 , Q 2 ) collective coordinates x = x(Q) : bead-bead distance as function of the global internal coordinates (fluctuations of bead 2 out of the pipette-trap axis are ignored …) Energetics (mean-field configurational potential, possibly time-dependent): V (Q, t ) = V0 (Q) − F (t ) x(Q) Vext (Q, t ) intrinsic internal potential external potential corresponding to the applied force: − ∂Vext (Q, t ) ∂x ≡ F (t ) x = x (Q ) Stocastic dynamics of Q can be described by means of - Generalized Langevin Equation: single-trajectory view to generate Q(t) - Non-stationary Fokker-Planck equation: probabilistic view to generate p (Q, t ) Generalized Langevin Equation ɺ = ∂ ⋅ D(Q) − D(Q) ∂V (Q, t ) + [ 2 D(Q) ]1/ 2 η(t ) Q ∂Q ∂Q D(Q) : configuration-dependent diffusion matrix η(t ) : array (column-vector) of independent sorces of white noise ηi (t ) = 0 ηi (t1 )η j (t2 ) = δ i, j δ (t1 − t2 ) force-dependent potential. Force-ramp experiments brings the non-stationary character Fokker-Planck ∂ ∂ ∂ + β V (Q , t ) / k B T p (Q, t ) = ⋅ D(Q)e−V (Q,t ) / k BT e p (Q, t ) ∂t ∂Q ∂Q Initial states: the system is “picked” from the pool of configurations at thermal equilibrium distribution in the asbsence of force (mainly, the native “folded” state for biopolymers…). p(Q, 0) = peq (Q, F = 0) = e−V0 (Q) / k BT ∫ dQe−V0 (Q) / kBT Modelling of the system (internal friction/energetics) D(Q) MATCHING ? adjustment of the system parameters up to get a good matching simulation of an ensenmble of Ntr trajectories Q(t) V0 (Q) distribution function p (Q, t ) Calculation of any average “temporal trace” yielded by the experiment. For example, a time-dependent mean of a function f (Q) : N 1 tr f t = dQ f (Q) p(Q, t ) ≃ f (Q(t )i ) Ntr i =1 Application of specific forces (in force-clamp) or the whole “force history” (in force-ramp) enhance the amount of information… ∫ ∑ Determination of temporal traces and their average Experiment performed with the chosen protocol A crude simplification: reduce the real multidimensional dynamics of Q to effective one-dimensional fluctuations of the end-to-end distance x (which can be directly measured during a single run). Effective one-dimensional potential: Veff ( x, t ) = Veff ,0 ( x) − F (t ) x Paradigmatic case: bistable profile adopted for biopolymers with Native (folded) / Unfolded “basins” of configurations. Two wells separated by an energy barrier. Veff ,0 ( x) U N − F (t ) x Veff ( x, t ) The external bias F > 0 heavily affects the profile (location and height of the energy barrier). - N→U transition is favoured; - N←U transition is penalized x Simple key-equation can be obtained for force-clamp and force-ramp protocols Force-clamp protocol (F = const) Phenomenological Bell’s model for kinetic constants N ↔ U [G. I. Bell, Science 200, 618 (1978)] *   Fx ku ( F ) = cM ku0 exp  u   k BT  *   − Fx  f  0 k f ( F ) = cM k f exp   k T  B  unfolding k 0f cM : dimensionless factor dependent on the overall machinery setup (kind of beads, molecular handles, etc.) ku0 , k 0f : intrinsic kinetic constants for the unfolding/folding processes in the absence of external forces folding k 0 u N U This is nothing but the application of Kramers theory under the assumptions that: - locations of wells’ minima and barrier are not sensibly changed by the bias; - curvatures at minima/maxima are also unchanged * u x x * f x In force-clamp, the molecule is initially at equilibrium and the given force is suddenly applied (force-jump in unfolding or force-decrease in refolding). In any run, the dwell-time (time passed until the first transition N↔U is detected) is recorded. From the distributions of the dwell-times, constructed by hystograms from counts, the force-dependent kinetic constants ku(F) and kf(F) are extracted. Monoexponential decays are normally observed. [see J. Liphardt et al, Science 292, 733 (2001)] Bell’s model ln k k f (F ) U→N ln(cM k 0f ) dwell-times repetition for several force values x * f N→U xu* ku ( F ) ln(cM ku0 ) dwell-times F “Distances” of the “transition state” from N and U wells’ minima are obtained! Very rough approach, but it gives reasonable outcomes! Some examples… Study of the four-α α-helix acyl-CoA binding protein [P. O. Heidarsson et al, J. Chem. Am. Soc. 134, 17068 (2012)] Extension along 1-86 xu* = (5.5 ± 0.6) nm x*f = (6.5 ± 0.3) nm (figure adapted from the original one) Study of the Prion Protein PrP – elongation along C-terminus, N-terminus [H. Yu er al, PNAS 109, 14452 (2012)] xu* = (10.4 ± 0.6) nm (figures adapted from the original ones) Force-ramp protocol (F(t) according to a prescribed protocol) A linear increse/decrease of the force with rate r = dF / dt = const is normally applied. -In the unfolding runs, the molecule is picked at equilibrium and elongated under force F (t ) = r t up to a certain Fc - In the refolding runs, the molecule is picked at equilibrium (with force Fc) and the force is reduced as F (t ) = Fc − r t Typical Force-elongation curves unfolding/refolding Bacteriophage T4 lysozyme Figures taken from the talk of C. Cecconi at Workshop BIOSTRUCT09, Firenze 2009 Quick changes of extension (“rips”) in the stretching runs are related to transitions N → U. Notice that the values of the forces at the unfolding event are distributed! Key-quantities which can be achieved from raw force-extension curves 1) Force distribution functions pu ( F | r ) , p f ( F | r ) pu ( F | r )δ F : probability that the unfolding transition occurs an applied force between F and F+δF , under force-rate r p f ( F | r )δ F : probability that the re-folding transition occurs at an applied force between F and F+δF , under force-rate r (starting from equilibrium at a given initial force F0) Normalizations: Fc ∫ p ( F | r )dF = 1 u 0 Fc ∫p 0 f ( F | r )dF = 1 2) Populations of Native and Unfolded states ϕu ( F | r ) , ϕ f ( F | r ) ϕu ( F | r ) : fraction of molecules still present in the folded state, during the unfolding, at applied F under force-rate r ϕ f (F | r) : fraction of molecules still present in the unfolded state, during the re-folding, at applied F under force-rate r (starting from equilibrium at a given initial force F0) These populations follow from integration of the associated force-distributions: Fc ϕu ( F | r ) = 1 − ∫ pu ( F ' | r )dF ' 0 Fc ϕ f ( F | r ) = 1 + ∫ p f ( F ' | r )dF ' F Analysis in terms of populations is simple; linearized relations are found: Fxu* ln [ − r ln ϕu (r , F ) ] = ln ( cM k k BT / x ) + k BT * Fx f ln  − r ln ϕ f (r , F )  = ln ( cM k 0f k BT / x*f ) − k BT 0 u * u (*) [compare with Eqs. 3 and 4 in JACS 134, 17068 (2012)] Populations ϕu (r , F ) and ϕ f (r , F ) can be obtained by counting the number of times that a molecule is still folded/undolded at force F in the unfolding/refolding direction (or they can be obtained by integrating the force-distribution-functions, as shown…) From data-fit with eqs (*) one can extract, as in force-clamp experiments, the location of the transition state for the N↔U conversion. [figures adapted from P. O. Heidarsson et al, J. Chem. Am. Soc. 134, 17068 (2012)] r = 15 pN/s r = 45 pN/s r = 8.5 pN/s r = 15 pN/s xu* = (5.2 ± 0.2) nm x*f = (6.7 ± 0.1) nm After all, force-clamp (jump) and force-ramp yield comparable results. Mechanical Tanford β-value (here for the unfolding direction): xu* mβTu := * xu + x*f Usually it ranges between 0.05 and 0.1; for “pliant” globular proteins stabilized by short-range interactions it may reach higher values, ~0.3 A more subtle modelling (beyond Bell’s model) yields improved expressions to interpret force-clamp and force-ramp data. H. Yu et al, PNAS 109, 14452 (2012) O. K. Dudku et al, Phys. Rev. Lett. 96, 108101 (2006) 1 1 −  0  * ε ku ( F ) xu F    ku ku ( F )  pu (r , F ) ∝ exp  * − * 1 − ε  * r xu r  ∆V    xu r   ε: shape parameter of the (unbiased) potential at the energy barrier ∆V * : energy barrier with respect to the well minimum of the native state N ∆V U * 1 −1  ∆V * ε   xF 0 ku ( F ) = cM ku 1 − ε  exp  * ∆ V    k BT * u For ε = 1 one recovers the previous expressions. 1/ ε *      xu F 1 − 1 − ε   * ∆ V      Force-ramp on a different kind of molecule: single RNA fragment from the P5abc domain of Tethahymena thermophila – group 1 intron [a seminal paper from J. Liphardt et al, Science 296, 1832 (2002)]. Unfolding/Refolding at various force-rates r (from slow “almost reversible” to fast protocols …) 2-5 pN/s 52 pN/s Mention to the extension-clamp protocol (still a novelty…) [Theoretical pespectives in Y. Suzuki, O. K. Dubko, Phys. Rev. Lett. 110, 158105 (2013)] xM Fu ( x) Ff ( x) Force is expected to fluctuate between two average values… Here, the transition rates depend parametrically on the fixed extension x : ku ( x ) k f ( x) These are the experimentally achievable data (from dwell-times distributions) ku ( x ) k f ( x) modelling - essential traits of the effetive potential V ( xM | x) ( ∆V * and shape parameters…) - xM* ,u and xM* , f - average forces Fu ( x) , Ff ( x) - ratios ∆V * / ∆V *  F = Fu ( x)  , ∆V * / ∆V *  F = Ff ( x)  experimentally achievable also force-clamp measures are required to get these values! [see Eq. (3) of Phys. Rev. Lett. 110, 158105 (2013) for an expression of the dwelltimes] From data fit one expects to get, at last, the same parameters which are available from force-clamp and force-ramp experiments (mainly, the locations of the transition state from the folded and unfolded conformations). But here the theory and analysis seems to be more hard to manage! Still a novely, no applications yet. Example of “torque spectroscopy” with magnetic tweezers: Rotor Bead Tracking (RBT) technique Inspection of supercoiling in short DNA sequences-of-interest (SOI). F. C, Obestrass et al, Phys. Rev. Lett. 110, 178103 (2013). [See also F. C. Oberstrass et al, PNAS 109, 6106 (2012) for technical details] “Static RBT” -from the applied rotation θ of the magnetic bead, and measuring the rotation ψ of the rotor bead, the difference α is obtained; - the torque τ on the SOI is obtained as τ = KTDα where the proportionality factor KTD is determined by calibration with a reference segment. [figure adapted from the original one, PRL 110, 178103 (2013)] DNA strand of 50 repeated base-pairs G-C (Z50). Two linear regimes of torque-twist 1) low-twist: right-handed B-DNA persists 2) high-twist: left-handed Z-DNA is more stable transition B→Z [figure adapted from the original one, PRL 110, 178103 (2013)] Bulk vs. single-molecule spectroscopies. Let us discuss… Conventional “bulk” spectroscopies Force/torque spectroscopies - Bulk response (sample average) - Single-molecule response - Excellent statistics - Need of many repetitions - Small to large molecules can be inspected - Only “large enough” molecules can be handled - The external mean does not affect the energetics (say, only transitions amongst unperturbed quantum states are induced) - No molecular engineering is required (sometimes, optical or magnetic probes need to be attached) - The external mean heavily perturbs system’s energetics (stiffness constant κ, molecular handles that become part of the global system, etc…) - Molecular engineering is required to link molecule and device - Standard experimental setup/methods - Experimental setup/methods require expertises (not yet standard tools) - All stochastic motions affect temporal response/spectral features (lineshapes, relaxation times…). It may be hard to disentangle internal dynamics from overall motions. - Complex molecules (many degrees of freedom) are untractable in detailed way. -Internal dynamics are probed. Change of points of force/torque application allows one to select specific directions. However, the information really achievable is little. - Complex molecules (many degrees of freedom) are untractable in detailed way. Let’s make the point… 1) Mechanical “spectroscopies” offer the chance to probe internal energetics of a flexible molecule along selected directions, once conformational fluctuations at fixed parameters (fixed load, applied rate of load change, fixed extension, etc.) are properly accounted for. 2) Drastic 1D simplification need to be adopted to let the modelling be tractable on theoretical/computational grounds; fluctuations reduce to transitions Native↔Unfolded under the external bias. 3) Only geometric localization of the transition states N↔U can be likely assessed. At best, estimates of energy barriers can be attemped; effects of the external apparatus and molecular handles (recall the ominous “machine factors” cM) bring uncertainty. 4) The simplest force-clamp measures seem to be enough to yield in transparent way such a geometrical information. The achieved target is too limited, despite of the powerful of the mechanical manipulations tools… Question: Is there a way to exploit mechanical manipulations of single molecules to go beyond such a too much limited level and construct, directly from raw data (i.e., without need of modelling a priori), the energy profile along a controlled internal coordinate? The answer is yes, thanks to the formulation (second half of ’90s) of so-called “Work Fluctuations Theorems” (WFTs) like Jarzynski’s Equality and Crooks Fluctuation Theorem. Still the end-end elongation of the prion protein PrP. [H. Yu er al, PNAS 109, 14452 (2012)] Reconstruction of internal free-energy profile from force-ramp experiments with the use of Jarzynski’s equality as elaborated by Hummer & Szabo [G. Hummer, A. Szabo, PNAS 98, 3658 (2001)] See next chapter… PART II Work Fluctuations Theorems (Thermodynamics of fluctuating systems) Abstract formulation: an externally (partially) contrained system subject to thermal fluctuations fixed parameter/set of parameters (imposed constraints) sys Λ T (termal bath) tem x fluctuating configutational variables of the system Example: - fixed parameter: end-to-end distance of a protein controlled in a position-clamp experiment with optical tweezers - “relevant” fluctuating variables x: set of dihedral angles VΛ (x) : mean-field configurational potential at imposed constraints A( Λ) : Helmholtz free-energy at imposed constraints The key-relation [from the theory of the statistical “canonical” ensemble: matching with macroscopic Thermodynamics requires this link between Helmholtz freeenergy and ensemble average of configuration-dependent Boltzmann factors]: e 1 β= k BT − β A( Λ ) from now on. = ∫ dx e − βVΛ ( x ) “configurational” partition function Remark: energies are implicitly expressed apart of an additive constant…. Is there a way to achieve the profile of A( Λ) versus Λ directly from experiments? The answer is yes, thanks to Work-Fluctuations-Theorems. We shall deal here with main ones: the Jarzynski Equality (JE) and the Crooks Fluctuation Theorem (CFT) [1] C. Jarzynski, Phys. Rev. E 56, 5018 (1997) [2] G. E. Crooks, Phys. Rev. E 60, 2721 (1999); ibid, Phys. Rev. E 61, 2361 (2000) [3] Excellent reviews on Work-Fluctuations-Theorems: C. Jarzynski, Annu, Rev. Condens. Matter Phys. 2, 329 (2011); C. Jarzynski, Eur. Phys. J. B 64, 331 (2008) a configuration “picked” Steered the pool at equilibrium, transformations from i.e. drawn from Boltzmann distribution peq ,Λ0 (x) Λ0 x(t) x T Λ1 Λ(t) Λ0 x(ts) finite-time non-equilibrium transformation system at equilibrium with the thermal bath; the controllable parameters are fixed. time 0 ts duration of the “steered” transformation • Freely chosen deterministic protocol for the fully controlled parameters: Λ (0) = Λ 0 → Λ (t ) → Λ (ts ) = Λ1 • Fluctuations of the uncontrolled variables: x(t ) stochastic trajectory “Work” is the amount of energy exchanged between external “device” and system through controlled intervent on the variables Λ [Such an identification of “work” follows from macroscopic Thermodynamics; see discussion in G. E. Crooks, J. Chem. Phys. 130, 107101 (2009)]. Work done along the i-th (stochastic) trajectory: ts wi = ∫ dt 0 ∂VΛ (t ) (x) ∂t ∂VΛ (x) ɺ = ∫ dt Λ (t ) ⋅ x = x ( t )i ∂ Λ Λ=Λ (t ) 0 ts x = x ( t )i ɺ = dΛ ( t ) / dt Λ contribution δwi in the time-interval dt: work performed to change the internal energy from VΛ (x) to VΛ+δ Λ (x) if the uncontrolled variables are kept “frozen”. For completeness, “heat” is the amount of energy exchanged with the thermal reservoir through uncontrolled change (fluctuations) of variables x. Its infinitesimal amount reads: ∂VΛ ( t ) (x) δ qi = δ txɺ × ∂x x = x ( t )i xɺ = xɺ ( t )i The uncontrolled nature of trajectory x(t) makes that wi is a stochastic variable! ts wi = ∫ dt 0 ∂VΛ (t ) (x) ∂t x = x ( t )i work spread of work values time work ts slower protocol time ts If we imagine to repeate the same transformation an infinite number of times (under same conditions), we end up with a distribution of work values: pprot. ( w) normalized as ∫ dwp prot. ( w) = 1 at given initial/final state-parameters Λ0 , Λ1 , it depends on - kind of dynamics x(t) - kind of protocol Λ(t) This is a quite different scenario compared to macroscopic Thermodynamics! In the macroscopic world, if a chose a “path” of transformation (for example, a schedule to compress a gas by regulating the external pressure…) the amount of work is always the same! Here, system’s fluctuations during the transformation produce a distribution of performed works! The messagge: Thermodynamics at the molecular scale must consider fluctuations. Expected transition from micro- to macro- systems Molecular-scale systems Macroscopic systems Finite-”size” system (few degrees of freedom x): fluctuations a relevant! Very many degrees of frredom (fluctuations of negligible amplitude) pprot. (w) Dirac’s Delta-function pprot. ( w) w wmacro w δ ( w − wmacro ) By increasing the system’s “size”, the work-distribution-function shrinks: at given schedule of transformation, always the same work is performed! Is there some case where such a shrink happens also for microscopic systems? Yes, it happens if the change Λ0 → Λ1 is infinitely slow: system and external device are always (almost) at equilibrium. Molecular-scale systems finite-time protocol (irreversible transf.) infinitely slow protocol (reversible transf.) pprot. (w) w pprot. ( w) ∆A w = A( Λ1 ) − A( Λ 0 ) δ ( w − ∆A) Intuitively: the trasformation is so slow that the system explores (by means of fluctuations) all possible internal conformations x at any actual state of the controlled parameters! By repeating the transformation, the performed work is always the same (equal to ∆A), and it does not even depend on the particular path to go from Λ0 to Λ1 . Theoreticians’ target: is there a way to exctract ∆A from a work-distributionfunction pprot. ( w) also if I adopt a finite-time protocol? This is practical need: we always operate in a finite time! Can we get out ∆A also in this case? If so, an equilibrium information (difference of free-energy) would be obtained from non-equilibrium transformations! This turns out to be possible, if we make the following assumptions about the uncontrolled dynamics on x: 1) Dynamics are of Markov type (memory-less process); Brownian motion belongs to this category. 2) If the protocol is stopped at a general Λ* reached at time t*, these dynamics make that the non-equilibrium distribution p (x, t*) relaxes to the “underlying” equilibrium distribution: lim p ( x, t ) = peq ,Λ* ( x) t *< t →∞ “uderlying” Boltzmann distribution On these basis, it can be demonstrated that the following identity holds exactly: e − β [ A ( Λ1 )− A ( Λ 0 )] JE = ∫ dw p prot. ( w) e −βw Jarzynski Equality [C. Jarzynski, Phys. Rev. E 56, 5018 (1997)] For a simple proof, see sect. 6.2 of C. Jarzynski Eur. Phys. J. B 64, 331 (2008). A( Λ1 ) = A( Λ 0 ) − β −1 ln ∫ dw pprot. ( w) e − β w Notice that the end-point Λ1 can be freely changed: a whole free-energy profile can be constructed! Crucial remark: Recovering of the Second Principle of Thermodynamics at molecolar level Exploit the Jensen inequality: for any convex function f ( x) (i.e. d 2 f ( x) / dx 2 ≥ 0) : Apply it to the convex function e e − β ∆A JE = −β w : f ( x) ≥ f ( x ) (averages on any distribution on variable x) −βw −β w −β w dw p ( w ) e ≡ e ≥ e ∫ prot. w ≥ ∆A Second Principle (under isothermal conditions) if the work is intended as the average work. By letting growing the system’s size, w → wmacro … By taking into account system’s fluctuations one gets an equality (the JE) rather than an inequality (Second Principle)! The inequality comes out if we miss the details of the spread of work values, and focus only on the average work. Back to applicative level. Other important remarks: 1) Strength of the JE: the result does not depend on the kind of protocol! (In fact, the free-energy is a state function: the free-energy difference is unaffected by the path/protocol to go from Λ0 to Λ1). If we change the protocol (path and local speed), the profile of pprot. ( w) changes but A( Λ1 ) − A( Λ 0 ) is always the same. 2) Something more subtle… Also the specific kind of (Markov) stochastic dynamics x(t) does not affect the outcome: dynamical aspects are disentangled by thermodynamical properties [you would realize this by looking at the demonstration of the JE…]. This fact is exploited in simulations when experiments are replaced by “steered molecular dynamics” or “steered Langevin/Monte Carlo dynamics”: fluctuations are artificially accelerated to get a more efficient sampling of the configurational space (see below…) How can I get the pprot. ( w) to insert into the JE? A single run is not sufficient! The distribution must be constructed (by hystograms) from an ensemble of repeated runs generated at same conditions (system initially at equilibrium, same temperature, same protocol). We need an ensemble of replicas of non-equilibrium transformations! A closer look at the JE: an intrinsic pitfall pprot. ( w) e − β [ A ( Λ1 )− A ( Λ 0 )] JE = −βw dw p ( w ) e ∫ prot. w low-work region Low-work trajectories dominate in the integral (they contribute with high e but these trajectories are rarely sampled (tail of the distribution!). −β w ), The idea is to sharpen the work-distribution-function (around ∆A) by slowering the protocol (i.e., by letting the transformation “less irreversible”, better: the trajectories “less dissipative”…) In the practice, the construction of pprot. ( w) is skipped: transformation is repeated a number Ntr of times (under same conditions/same protocol) and one gets directly:  1 −1 A( Λ1 ) ≃ A( Λ 0 ) − β ln   N tr N tr ∑e i =1 − β wi    The accuracy of the outcome depends on the kind of protocol and on the number of repetitions: as Ntr → ∞ the result is exact. For a finite number Ntr , “rare” low-work trajectories which contribute largely, may be not encountered! This gives a systematic bias to the outcome: ∆A is overestimated! To reduce such a bias, the protocol should be slowered… Think about this… For errors estimate in JE applications see: - C. Jarzynski, Phys. Rev. E 56, 5018 (1997) - D. M. Zuckerman,T. B. Woolf, Phys. Rev. Lett. 89, 180602 (2002). Crooks Fluctuation Theorem (CFT) A “bidirectional” transformation is considered: forward (F) and reverse (R) processes one with a “conjugate” (time-reversed) protocol of the other. Forward: equilibrium at Λ 0 → Λ F (t ) → Λ1 Reverse: equilibrium at Λ1 → Λ R (t ) = Λ F (ts − t ) → Λ 0 Under the same conditions for the JE validity, Crooks obtained pF ( w) = pR (− w) e − β  A( Λ1 )− A( Λ 0 )− w G. E. Crooks, Phys. Rev. E 60, 2721 (1999) Notice that the JE follows as a corollary of the CFT: pF ( w) = pR (− w) e − β  A( Λ1 )− A( Λ0 )− w ⇒ pF ( w)e − β w = pR (− w) e − β  A( Λ1 )− A( Λ 0 ) By integrating both members on w: −β w dw p ( w ) e =e ∫ F − β  A( Λ1 )− A( Λ 0 )  ∫ dw pR (−w) = e =1 − β  A( Λ1 )− A( Λ 0 )  This is the JE! The crossing point of the two distributions falls at wcross ≡ A ( Λ1 ) − A ( Λ 0 ) . It gives directly the free-energy difference! Computer simulation of unfolding/refolding of deca-alanine by means of “steered” Molecular Dynamics [Figure adapted from P. Procacci et al, J. Chem. Phys. 124, 164101 (2006)] Use of the bidirectional CFT gives better results than the unidirection JE method: more accurate results. free energy By changing the end-state (Λ1 here is simply the end-to-end sistance), the full profile of free-energy is obtained: [Figure adapted from P. Procacci et al, J. Chem. Phys. 124, 164101 (2006)] In the essence: JE and CFT yield the free energy profiles versus Λ (with respect to some reference state Λ0) as a direct experimental outcome from repeated transformations on the molecular system. No modeling is needed! What we need is a way to quantify the amount of work per trajectory… Devices for molecular manipulations allow to get it. For example, the work performed can be obtained by integrating the force vs. extension in force-ramp runs… A problem: in this way we should end up with the free-energy profile of the global system made of molecule + external apparatus! A “reweighting procedure” is needed to extract only the system’s free-energy profile. The external apparatus (eg., a cantilever, a tweezer) perturbs the system by inserting a “bias” term in the potential. The effective potential is V eff , Λ ( x ) = V ( x ) + V bias , Λ ( x ) unconstrained system We can only control mechanical parts of the apparatus, not of the molecule! Hence, the Λ are parameters of the apparatus. For example, a single λ may be the distance pipette-trap center, elongated at constant velocity v: V eff ,λ ( x ) = V ( x ) + 1 2 κ (d (x) − λ ) , λ ( t ) = λ 0 + vt 2 controlled distance real end-to-end distance of a biopolymer The work can be obtained from the real end-to-end trajectory (but also from t force-elongation curves): s wi = κ v ∫ dt [ d (x(t )i ) − vt ] 0 In this way we will end up with the global free-energy profile Aeff (λ ) , not with the intrinsic A( d ) that we want! A “reweighting procedure” is needed to get it… G. Hummer, A. Szabo, PNAS 98, 3658 (2001). A seminal work A computational study (Steered Molecular Dynamics and application of the JE) Energetics of glycerol through aquaglyceroporin (Glpf) M. Jensen et al, PNAS 99, 6731 (2002) Back to the end-end elongation of the prion protein PrP. [H. Yu er al, PNAS 109, 14452 (2012)] A recent work Integration of force-elongation curves works use of JE with reweighting procedure (Hummer-Szabo) With such uni-dimensional free-energy profile, one can construct a model for stochastic dynamics N↔ ↔U (this is half of the story: remember that also information about friction is needed in Langevin or Fokker-Planck equations!)