Universal Aspects of Barrier Crossing Under Bias

The thermal activation process by which a system passes from one local energy minimum to another is a recurring motif in physics, chemistry, and biology. For instance, biopolymer chains are typically modeled in terms of energy landscapes, with folded and unfolded conformations represented by two distinct wells separated by a barrier. The rate of transfer between wells depends primarily on the height of the barrier, but it also depends on the details of the shape of the landscape along the trajectory. We consider the case of bias due to an external force, analogous to the pulling force applied in optical tweezer experiments on biopolymers. Away from the Arrhenius-law limit and well out of equilibrium, somewhat idiosyncratic behavior might be expected. Instead, we identify universal behavior of the biased activated-barrier-crossing process and demonstrate that data collapse onto a universal curve can be achieved for simulated data over a wide variety of energy landscapes having barriers of different height and shape and for loading rates spanning many orders of magnitude.


I. INTRODUCTION
Thermally activated barrier crossing [1][2][3][4] is a ubiquitous and highly consequential process in physics, chemistry, and biology.An understanding of the factors that influence the rate of barrier crossing [5][6][7][8] is necessary for the interpretation of experiments that attempt to infer barrier height and shape from measurements of the escape rate.An important specialization is the case in which the barrier is diminished by an applied force, with the escape rate enhanced accordingly.
Experimental access to escape rate information in the biochemistry context has been revolutionized by the development of single-molecule force spectroscopy [9][10][11][12][13][14][15][16][17][18], in which a mechanical load is applied across a single molecule using an atomic force microscope or optical tweezers.In the energy landscape picture [19][20][21][22][23][24], molecular motion is viewed as a Brownian diffusion process over a free energy surface [25], parameterized by the conformational degrees of freedom.The landscapes for biologically relevant sequences contain distinct, barrier-separated wells corresponding to various folded and unfolded conformations.The rate of transition [26][27][28] from one well to another depends primarily on the height of the intervening barrier but also depends on its shape.
In pulling experiments, where molecular extension serves as a natural reaction coordinate, the multidimensional energy landscape covering the full comformational space is projected onto an effective one-dimensional energy profile that encodes some features of the full landscape and that reproduces the folding dynamics [29][30][31][32][33]. Numerous studies have been carried out to explore the unfolding process under the application of constant and time-varying pulling forces [34][35][36][37][38][39][40][41][42][43][44].A key experimental goal is to be able to reliably reconstruct the effective one-dimensional free energy profile from measurements of an ensemble of escape events [45][46][47][48][49][50].
There are additional complications that may arise because of the multidimensional nature of the landscape [51][52][53][54][55][56], but we consider the common case in which the landscape can be meaningfully projected onto a one-dimensional effective energy in a well-chosen reaction coordinate [56].
The purpose of this article is to describe universal aspects of the biased activated-barrier-crossing process that we have uncovered in numerical simulations of various one-dimensional potentials.Our work points the way to an alternative data analysis technique that would allow for the determination of otherwise unknown landscape details by overlaying data from multiple experiments and adjusting free parameters until the scattered data align along a common curve.
The concept of universality comes to us from the study of critical phenomena [57].In that context it allows us to understand how phase transitions can be characterized and grouped into families according to common critical exponents, wholly independent of the microscopic details of the underlying models; it also explains the existence of scaling relations that govern how thermodynamic quantities behave in the vicinity of criticality.An important mark of universality is that data from different models or different physical systems can be plotted in reduced variables so that they collapse onto a single universal curve [58][59][60].
Criticality has previously been invoked by Singh, Krishan, and Robinson in the context of the unbiased-activatedbarrier crossing problem [61,62].They considered the non-Markovian crossing of a quadratic barrier, where the frictional term in the Langevin equation includes a memory kernel with a long time scale.The authors proposed a scaling hypothesis, making analogy with the criticality of the Ising model, and were able to derive scaling relations for the reduced rate near a critical value of the memory kernel time scale.
Our approach here is rather different.We focus on the relative change in the escape rate as a function of an applied pulling force-both for uniform pulling ( constant) and steady loading ( =  with  constant).We propose that  exists alongside two other important force scales and that the two independent ratios that can be formed serve as arguments to a scaling function.We have carried out Langevin-type simulations of a particle in a one-dimension energy potential, coupled to a heat bath.Many thousands of instantiations provide us with a large data set that offers good coverage of the model space.What is so striking is the almost unreasonable effectiveness of the scaling ansatz, which appears to be valid over  a huge variety of well shapes and barrier heights and over loading rates spanning many orders of magnitude.

II. SCALING ANSATZ
We argue that the barrier-crossing process is controlled by the relative magnitudes of three intrinsic force scales: the typical thermal force that provides the kick out of the well (  ≈ 1/ ‡ ); the larger applied force required to fully extinguish the barrier (  ≈  ‡  ‡ ); and the pulling force used as an external bias ().In our notation,  is the inverse temperature, and  ‡ ,  ‡ are the barrier distance (see Fig. 1) and effective curvature [63].In particular, we propose that the rescaled, logarithmic relative escape rate (  /  ) ln[ ()/ 0 ], when plotted against the reduced pulling force /  , collapses onto a universal curve.The function  () is the escape rate associated with the potential landscape under bias, and  0 is the corresponding rate in the untilted landscape.
More precisely, the claim is that Up to subleading corrections (characterized by an exponent  > 0), the escape behaviour is controlled by a universal function Y(, ) that satisfies Y(, 0) ∼  + ( 2 ).The implication of Eq. ( 1) is that a plot of (  /  ) versus  = /  should produce data collapse regardless of the microscopic details of the simulation.
In the context of dynamical pulling, there is another useful analysis.A population of particles trapped in the originating well is depleted according to −/ =  ()(), where the right-hand side of the equality is a product of the instantaneous escape rate and the current population.The population's half-life is characterized by ln Here, F =  t is the typical applied force that is in effect during barrier transit, and t is the median elapsed time for escape.It follows from Eq. (1) that F (measured with respect to the thermal force   ) must be a monotonic, universal function of  =  (measured with respect to  0   , a loading rate threshold defined by the thermal processes in the potential well).Hence, there is an additional data collapse analysis that can be used to independently test the validity of the scaling hypothesis.
In our numerical experiments, the external bias is applied in two ways: (i) as a time-invariant pulling force of constant strength and (ii) as a linearly time-varying force with a constant loading rate.In the case of constant pulling, the system is prepared in the equilibrium state of the tilted energy profile [viz.Ũ () =  () − ] and remains in thermal equilibrium throughout the simulation.In the case of steady loading, the system is prepared in the equilibrium state of the unbiased profile ( = 0 for all  ≤ 0), but as time elapses it is driven away from equilibrium ( =  for all  > 0) in proportion to how much  exceeds  0   .
In both cases, the role of  is to gently tilt the landscape (statically or dynamically), as depicted in the upper panel of Fig 1 .There, the purple curve depicts the potential profile in its unbiased state; the green curve shows the profile after application of the external bias.Generically, the force-dependent values of the barrier distance  ‡ () and barrier height Δ ‡ () are monotonic decreasing in , and hence the barrier crossing process becomes energetically less costly (and crossing events more frequent) as the external bias is ramped up.The freeenergy minimum is progressively destabilized and disappears entirely at the threshold for barrier extinction.

A. Locally quadratic approximation
We consider a one-dimensional, double-well energy landscape  () with minima on the left and right, at positions   and   , separated by a barrier at   .A barrier of height Δ ‡ =  (  ) −  (  ) impedes transitions from left to right.By definition  ′ (  ) =  ′ (  ) =  ′ (  ) = 0.In the locally quadratic approximation, we assume where   =  ′′ (  ) and   = − ′′ (  ) are measures of the curvature at the bottom of the well and at the top of the barrier.
With the application of a bias force , the extrema of the tilted landscape Ũ () =  () −  are found as follows: At this level of approximation, the bias-induced shifts in the extrema are linear in .In response to the applied force ( > 0), the well basin moves to the right and the barrier peak moves to the left: The inverse spring constants, 1/  and 1/  , represent the compliance of the reactant and transition state; see Eq. ( 4) of Ref. 55 and the accompanying discussion.The two points eventually coalesce when x = x , i.e., when The particular force value at which Eq. ( 5) holds is the barrier extinction force  ‡  ‡ .We follow the usual practice of decorating with a double-dagger superscript any quantity that is defined with respect to the barrier and the originating well.This includes the barrier distance  ‡ =   −   and the effective curvature In order to find an expression for the barrier height that is consistent with the approximation in Eq. ( 2), we must match the two piecewise quadratic curves.We do so at the point of common slope, where The reference position is a weighted average satisfying   ≤  * ≤   .The height of the barrier in the untilted landscape ( = 0) is estimated to be Formally, the barrier extinction force is given by the derivative of the barrier height with respect to the barrier position.At the level of approximation of Eq. ( 9), we have

B. Higher-order corrections
There are various forms of "extended Bell theory" [64-66] that offer systematic refinements to the transition rate under bias.These are generally structured order by order in the applied force, as per Eq.(1) of Ref. 64.Indeed, we can improve on the derivation in Sec.III A by including further contributions to the energy landscape expansion: As in Eq. ( 2), the upper expression in Eq. ( 5) is for  ≃   ; the lower corresponds to  ≃   .In addition to the two local curvatures,   and   , we have also defined measures of the skew [  = − ′′′ (  ) = − (3) (  ) and   = − (3) (  )] and the kurtosis [  =  (4) (  ) and   = − (4) (  )].
The positions of the shifted extrema are once again determined by 0 =  ′ () − .This demands that the expression vanish for both  = ,  = +1 and  = ,  = −1.Ensuring that it does so leads to and hence to an expression for x − x , the barrier distance in the tilted energy landscape: As a convenience, we have adopted the notation If we then repeat the analysis used previously, we can produce expressions for the barrier height and barrier extinction force that are analogs of Eqs. ( 9) and (10): and It is helpful to distinguish the barrier height expressions in Eqs. ( 9) and ( 16) by the labels Δ ‡ quad and Δ ‡ .Their ratio is simply the shape parameter defined by Dudko, Hummer, and Szabo [67]: That is to say, 1 −  encodes deviations from the behaviour of the purely quadratic model (in which   =   = 0, etc.).Insofar as Eq. ( 18) is a fast-converging power-series in  ‡ , with each subsequent term much smaller than the previous, it makes sense to view the subleading term on the right-hand-side of Eq. ( 18) as a proxy for those deviations: Hence, via Eq.( 17), the extinction force can be approximated by its quadratic-model value [Eq.(10)] up to rescaling by a shape-dependent factor: Typical values for smooth energy profiles (2/3 ≲  ≲ 6/5) suggest 0.5 ≲ (3 − 1)/2 ≲ 1.3, so we expect the true extinction force value to be never more than a factor of two away from  ‡  ‡ .Of course, when  () is known, it is straightforward to compute the exact extinction force numerically.

C. Universality of the biased escape rate
The calculations in this section are meant merely as a motivation for the two-force-scale arguments we make in the paper.We assume Langevin behaviour with moderate to strong friction and ignore the complications of non-Markovian dynamics.Following Kramers, the escape rate from the left well of the untilted energy landscape is The corresponding expression for the tilted case can be produced by substituting  (  ) → Ũ ( x ).If we expand around the  = 0 case and collect terms order by order within the argument of the exponential, we arrive at In this equation, the terms proportional to  come from the exponential in Eq. ( 21), whereas the temperature-independent contributions originate under the radical of the prefactor in Eq. (21).Since Eq. ( 22) can also be expressed as Although the terms replaced here by  (1/, ) may be small-in the limit of very large barrier height (Δ ‡ (0) ≫ 1) or possibly even for particular shapes of the energy landscape-they are generally not negligible.As we will see, a renormalization of the force scales   and   from their "bare" values is necessary to absorb the discrepancy.
An important insight is that the logarithmic relative rate can be written in the form In this expression we have introduced two new dimensionful coefficients (with units of force), defined according to 1 along with a dimensionless constant .Matching the  ( 3 ) terms in Eqs. ( 22) and ( 25) and invoking Eq. ( 19), we identify  = 1 −  + • • • , with the elision hiding additional terms that are shape and temperature dependent but small; specifically, The advantage of the rewriting in Eq. ( 25) is that we have picked out two force scales,   and   , whose magnitude is determined-up to modest renormalization by   and   by  (0)  = 1/ ‡ and  (0)  =  ‡  ‡ .Equation (26) implies To give a more physical picture, we interpret   as the typical thermal force that provides the kick out of the well and   as the applied force required to fully extinguish the barrier: The ratio of the two force scales is A key observation is that if we view the escape rate as a function of a reduced applied force  = /  , measured in units of the barrier extinction force scale, then Eq. ( 25) transforms to Note that the terms at order  and  2 are wholly independent of the details of the system.(/2)  3 is the leading nonuniversal term, but even there [as per Eq. ( 27)] the shape dependence is quite weak and the temperature dependence almost negligible.This means that truly idiosyncratic contributions do not show up until order  4 , and those we expect to be heavily suppressed just by power reduction; in practice,  = /  < 1, since escape almost always precedes complete elimination of the barrier.Moreover, since physical considerations demand that the escape rate increase with , it is legitimate to apply a series acceleration transformation by which Eq. ( 30) is expanded in terms of some function of  that is monotonic increasing but slower growing than the monomial  itself; one might consider  /(1 +  /2) (as in Ref. 63) or ln(1 +  ), say.Then  − (1/2)  2 − (/2)  3 + • • • can be recast as or Since we expect −1/5 ≲  ≲ 1/3 (and often | | ≪ 1), Eqs. ( 32) and ( 31) are close to being universal even up to order three.This leads us to posit that the logarithmic relative escape rate has a form reminiscent of the finite-size scaling ansatz of a critical state: viz., the form given by Eq. ( 1).Then, since   /  ≈ 1/( ‡  ‡ 2 ) ≈ /(2Δ ‡ ) ≪ 1, the quantity (  /  ) () should collapse onto a universal curve when plotted against  = /  : Other combinations of (  /  ) () may bring about an even cleaner coincidence.For example, ln(1 As for the utility, imagine that there is a set of escape rate measurements for which the underlying profile  () is unknown.Then, even without a model or fitting form for  (), we can still engineer graphical collapse of the data onto a common curve by rescaling and careful adjustment of the free parameters   and   .

D. Data collapse of the rupture force
A population () of systems prepared in a well and subject to an escape rate  () is subject to the rate equation  = −.If the pulling force increases linearly in time, with a constant loading rate , then The time t for half the population to escape is given by In the Bell-Evans picture [34,35], which supposes that the biased rate is simply  () =  0 exp  ‡ , Eq. ( 36) becomes Then F, the force at half-life, is On the other hand, if the escape rate is represented using the universal part of Eq. ( 32), via ln[ () The omitted terms [denoted by • • • in the exponent of the last line of Eq. ( 39)] are ones that become negligible at low temperature and large barrier height; in that same limit, we can formally take  ‡  ‡ → ∞, which allows us to recover the Bell-Evans expression,  () →  0 exp  ‡ .
We need not resort to such a limit, however, since the halflife can be solved analytically: This corresponds to an average rupture force A useful resummation is the first term of which is identical to the right-hand-side of Eq. ( 38), up to the renormalization  (0)  →   .In order to put Eq. ( 42) into a scale-invariant form, we define the half-life pulling force with respect to the thermal force scale, f = F/  , and a dimensionless loading rate,   = /(   0 ).This leads to f = ln 1 +   ln 2 In general,   /  ∼ 2Δ ‡ / is not small.But so long as

IV. NUMERICAL RESULTS
We carried out a thorough and comprehensive Langevin simulation study.At the start of each run, the system was prepared in a properly equilibrated state: an initial position and velocity were drawn from the heat bath distribution of the appropriate energy profile, with the constraint that the particle be situated on the originating-well side of the barrier.Forward evolution was carried out with adaptive time steps taken small enough that the discretization error could be shown to be negligible.The simulation made use of a high-quality, long-period pseudorandom number generator that guaranteed the statistical independence of the instantaneous thermal forces (based on Gaussian-distributed noise  () that is unbiased, ⟨ ()⟩ = 0, and uncorrelated except at identical time slices, ⟨ () ( ′ )⟩ = ( −  ′ )).
We considered seven different potentials having shape parameter [67]  = 0.66, 0.75, 0.83, 0.9, 1, 1.1, 1.2; these values step through the full range of possibilities for smooth potentials based on polynomials.This family of energy potentialstranslated and rescaled to coincide at the bottom of the originating well and at the top of the barrier so as to emphasize the shape difference-is displayed in the lower panel of Fig 1 .We also considered eight barrier regimes, with 1/Δ ‡ =   /Δ ‡ taking values 0.25, 0.3, 0.35, 0.40, 0.45, 0.50, 0.55, 0.60, a list that includes temperatures high enough (or, equivalently, barriers low enough) to be outside the range of validity for pure Arrhenius-law behavior.Simulations were carried out in both the constant-force and steady-loading modes, with relative applied forces (/  ) and relative loading rates (/ 0   ) each spanning nearly ten orders of magnitude.For each run, the trajectory leading to barrier traversal was captured and analyzed.
The numerical simulations were carried out using a modified version [68] of the standard Verlet algorithm [69].In each run, a critical time   , taken to be either the first-passage time of the particle over the barrier or the moment at which the barrier vanished, was recorded.For each energy potential profile and simulation mode, 3000 instantiations were generated.
In the constant-force mode, the rate  () was computed from the mean escape time:  () = 1/ avg , where  avg = (1/3000) 3000 =1  ()  .In the steady-loading mode, the linear correspondence   =   gave rise to 3000 critical force values, on the basis of which further analyses were performed.First, the   values were sorted to identify their median value, which corresponds to the half-life force F (the force at which half of a population of independent particles would have escaped the well).Second, the   values were bootstrapped [70] to obtain the cumulative probability distribution (  ) = ∫   0  () and probability density (  ) =  ′ (  ).Finally the value of  () =  (), the instantaneous rate of barrier crossing at a particular bias strength, was obtained using the relation  (  ) =  (  )/(1 − (  )) [67].The next step was to test the universality proposition by graphical means.We found strong evidence in its favor: the data collapse predicted by Eq. ( 1) is revealed in Fig. 2. In order to perform the conversion to reduced variables, each data point was associated with an individualized value of   and   .The former was obtained numerically, simply by solving for the applied force required to extinguish the barrier; the latter was estimated according to 1 Here, we have made use the notation introduced in Sec.III B. Figure 2(a) presents a linear-log plot of (  /  ) () = (  /  ) ln ( ()/ 0 ) versus  = /  .The data points for the steady-loading analysis are colored according to the simulation-specific loading rate, and one can observe the smooth progression of data-point placement, weak loading to strong, tracing out the universal curve from left to right.The constant-force data (black circles) show considerably less  scatter, but the two data sets are remarkably consistent.What makes this result so compelling is that the data collapse holds over a huge diversity of energy profiles and simulation conditions.We also remark that the steady-loading and constantforce approaches require quite different styles of simulation and analysis, but both yield the same underlying curve; Padé approximants fit to one or the other data set produce nearly identical functions.Figure 2(b) shows the same data plotted on a linear scale.This view highlights the behavior at large forces, a regime in which the barrier is already substantially reduced at the time of barrier traversal.Here, the false color emphasizes the diversity in barrier height regimes, and we can see that data collapse holds over a wide range of ratios   /  .
Figure 3 presents a wholly different data collapse scheme, based only on simulations performed in the steady-loading mode.There, the reduced half-life force f = F/  is plotted versus the reduced loading rate   = /( 0   ).It is worth emphasizing again that the complete data set comes from simulations with seven different potential landscapes covering the full range of plausible  values, 1/(Δ ‡ ) ranging from 0.2 to 0.6, and loading rates running from   = 10 −8 to 100.Despite encompassing a large collection of different systems in distinct physical regimes, these data show an astonishing degree of collapse.

V. CONCLUSION
We have argued for universality in the biased activatedbarrier-crossing problem and presented strong numerical evidence in favor of the existence of some underlying scaling function for  () = ln[ ()/ 0 ].Our simulated data show collapse onto a single curve when recast into suitably reduced coordinates.This is true for data generated in simulations operating over a wide range of bath temperatures, applied forces, and loading rates and over a family of potential landscapes with different underlying barrier shapes.
These observations suggest the utility of data collapse as a practical tool for analysis.While the original motivation for this work was the mechanical unfolding of biopolymers, the universality we have identified is widely relevant.It applies to situations across many branches of science, wherever the energy landscape picture is germane and the experimental setup involves barrier traversal assisted by active pulling.Our recommendation is that measurements of well-escape statistics be transformed to identify best values of the intrinsic force scales (from which can be inferred some combination of  ‡ ,  ‡ , Δ ‡ , and ).  and   are to be treated as free parameters and tuned until data collapse is achieved and the universal curve emerges.
To motivate our general approach, and to provide logical and formal scaffolding for the scaling argument, we have relied on analytic expressions for the biased escape rates.In Sec.III, these were computed within the context of Kramer's theory-in particular, the most straightforward version, which presumes moderate, Ohmic friction, does not include explicit finite-barrier corrections [71], and does not attempt a more sophisticated reconsideration of the Kramers-Grote-Hynes transmission factor [72].
Our demonstration of scaling, however, does not depend on Kramer's theory specifically nor on the quality of the analytical approximation.We emphasize that the data collapse in this work is achieved with escape rates determined empirically from Langevin simulations (as described in Sec.IV).That is to say, the scaling viewpoint we advocate is agnostic with respect to the well-escape model.We simply point out that, in principle, it should be possible to plot experimental measurements of the escape rate in reduced coordinates such that the data falls on a single curve; this relies on identifying system-specific values of   and   , the two force scales whose values must be varied to produce the data collapse.Details of the underlying potential can then be inferred from   and   , although that final step does introduce some dependence on how escape from the well is modeled.
We close with one final comment.It has already been established that the barrier vanishes generically according to Δ ‡ ∼ (  − ) 3/2 for bias forces in close vicinity of the barrier extinction force   [73,74].This is an important result.It leads to predictable, standardized behavior in that force regime and thus provides the basis for a kind of force spectroscopy [51].Nonetheless, we view this as asymptotic behavior (as  →   ) and not an example of universality in the traditional sense.The results reported in this paper are applicable to all applied forces less than   and do not dependent on any particular force limit.The data collapse, for instance, appears to be valid for  = /  ranging at least from 10 −7 to 0.4.The universality we advocate for here is an underlying commonality across many orders of magnitude that is revealed by rescaling.

FIG. 1 .
FIG. 1.(a) In the unbiased energy landscape [solid purple line, labelled  ()], a particle escaping (left to right) from the well must traverse a barrier of height of Δ ‡ =  (  ) −  (  ), over a distance  ‡ =   −  , where   is the position of the bottom of the left well and   is the position of the barrier peak.With application of an assistive pulling force ( > 0), the energy landscape tilts (solid green line) to favor the destination well to the right of the barrier.The pulling force causes the extrema to shift; the barrier height Δ ‡ () and barrier distance  ‡ () both decrease.(b) Seven energy profiles are depicted, rescaled so that the bottom of the well and top of the barrier coincide.The color key shows the shape parameter () values for each curve.

FIG. 2 .
FIG.2.The upper two panels show constant-force data only, (a) unscaled and raw and (b) rescaled to produce the predicted data collapse.The lower panels, like panel (b), show plots of (  /  ) () = (  /  ) ln ( ()/ 0 ) versus  = /  .These include all the output generated by our simulations (from both forcing protocols) and offer different views of the same underlying data set.Black circles correspond to simulations executed in the constant-force mode.Colored solid circles denote data from steady-loading runs.The green line is a low-order Padé approximant fit through the data points.(c) The horizontal axis uses a logarithmic scale.Color saturation increases with the relative pulling rate,   = /( 0   ).Numbers on the palette legend refer to the order of magnitude, ln   .(d) The horizontal axis is linear.Colors now represent   /  ≈ 1/ ‡  ‡ 2 ≈ /2Δ ‡ , which characterizes the barrier regime.
FIG.3.These data are from simulations performed with dynamic tilting of the potential at a constant loading rate.Shown here is f  = F/  , the force at half-life measured with respect to the thermal force scale, plotted against the effective pulling rate,   = /( 0   ).Each data point is a solid circle, colored as per the legend according to its   /  value.