Efficient two-dimensional control of barrier crossing

Modern advances in single-molecule biophysics make possible the precise spatial and temporal control of biological systems. Optical tweezers can be used to probe the conformational and energetic properties of biopolymers (DNA and RNA molecules) [1–6] and molecular machines (ATP synthase [7–9], kinesin [10–13], and myosin [14–17]). Additionally, computer simulations such as steered molecular dynamics have the freedom to fully control the molecular and trapping potentials. Despite the relative freedom of control, experiments and simulations rarely exploit the possibility of optimized control protocols, and the few that do are generally limited to optimization of a single control parameter [18–20]. In this letter, we design minimum-dissipation protocols for harmonic trapping potentials under two-dimensional control of trap center and stiffness allowing for the specification of the time-dependent mean and variance and resulting in significantly reduced dissipation compared to control over the trap center only. We find generic features that can be readily applied to biophysical experiments and simulations.

We are interested in describing micro- and nano-scale thermodynamic systems. Due to the small scale, thermal fluctuations play a significant role in the dynamics, and the systems are best described by stochastic thermodynamics. Stochastic thermodynamics typically describes the nonequilibrium transformation of heat, work, and entropy of small-scale fluctuating systems [21,22]. This field has seen significant growth over the last couple of decades, stemming from important results such as the Jarzynski equality [23] and Crooks fluctuation theorem [24].

One important implication of the Jarzynski equality is that it allows for the determination of equilibrium free-energy differences from nonequilibrium work measurements. However, the accuracy of the free-energy estimate decreases with increasing dissipation [25,26]. Designing protocols that dissipate less energy can therefore improve the accuracy of free-energy estimates [27]. Additionally, in order to maintain complex nonequilibrium order, molecular machines must operate rapidly, potentially incurring large energetic costs [28]. The dissipation incurred from rapid driving can be mitigated by designing less dissipative protocols [29–31], resulting in improved performance.

Linear-response theory has been used to derive a thermodynamic-geometry framework to guide design of single and multidimensional minimum-dissipation protocols in general stochastic-thermodynamic systems [29]. This framework has been applied to DNA-hairpin pulling experiments, demonstrating that one-dimensional designed pulling protocols significantly reduce dissipation [32]. The DNA-pulling process can be reasonably well described as a driven barrier crossing [33], where the barrier to be overcome is the transition state between folded and unfolded conformational ensembles. One-dimensional minimum-dissipation control of driven barrier crossing slows down as the trap crosses the energy barrier [34].

We explore two-dimensional control (of both trap center and stiffness) of driven barrier crossing. This greater control allows specification of both the time-dependent mean and variance of the position distribution, and results in a qualitatively distinct designed protocol (fig. 1). Such a designed protocol has jumps at the start and end that decrease in size as the duration increases, and slows down and tightens as it crosses the barrier, approximately linearly driving the mean and maintaining roughly constant variance throughout the protocol. For any duration, the designed protocols significantly improve performance in terms of both dissipation and flux compared to naive and one-dimensional control (fig. 2).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Consider a system obeying Fokker-Planck dynamics,

$$\begin{equation} \beta\gamma\frac{\partial p_{\bm{\Lambda}}(x,t)}{\partial t} = -\frac{\partial }{\partial x}\left\{\beta f_{x}[x,\bm{\lambda}(t)]-\frac{\partial }{\partial x}\right\}p_{\bm{\Lambda}}(x,t), \end{equation} \tag{ 1 }$$

describing the overdamped motion of a continuous degree of freedom x with damping coefficient γ driven by a time-dependent conservative force $f_{x}[x,\bm{\lambda}(t)] = -\partial V_\mathrm{tot}[x,\bm{\lambda}(t)]/\partial x$, for internal energy $V_\mathrm{tot}(x,\bm{\lambda})$. $p_{\bm{\Lambda}}(x,t)$ is the probability distribution over microstates x at time t given the control protocol $\bm{\Lambda}$. The system is in contact with a heat bath at temperature T such that the equilibrium probability distribution over x at fixed control parameters $\bm{\lambda}$ is $\pi(x|\bm{\lambda}) \equiv \exp\{\beta[F(\bm{\lambda})-V_\text{tot}(x,\bm{\lambda})]\}$, for free energy $F(\bm{\lambda}) \equiv -k_\mathrm{B}T\, \ln\sum_{x}\exp[-\beta V_\mathrm{tot}(x,\bm{\lambda})]$, where $\beta \equiv (k_\mathrm{B}T)^{-1}$ for Boltzmann's constant $k_\mathrm{B}$. The average excess work $W_\text{ex} \equiv W-\Delta F$ by an external agent changing control parameters $\bm{\lambda}$ according to protocol Λ is

$$\begin{equation} \langle W_\mathrm{ex} \rangle_{\bm{\Lambda}} = -\int_{0}^{\Delta t }\mathrm{d}t \,\frac{\mathrm{d}\bm{\lambda}^\mathrm{T}}{\mathrm{d}t} \langle \delta \textbf{f}_{\bm{\lambda}}(t) \rangle_{\bm{\Lambda}}, \end{equation} \tag{ 2 }$$

where superscript $\mathrm{T}$ denotes transpose. $\textbf{f}_{\bm{\lambda}} \equiv -\partial V_\mathrm{tot}/\partial\bm{\lambda}$ are the forces conjugate to the control parameters, and $\delta \textbf{f}_{\bm{\lambda}} \equiv \textbf{f}_{\bm{\lambda}} - \langle \textbf{f}_{\bm{\lambda}} \rangle_{\bm{\lambda}}$ the deviations from the equilibrium averages. Angle brackets $\langle \cdots\rangle_{\bm{\Lambda}}$ denote a nonequilibrium ensemble average given the control protocol $\bm{\Lambda}$ and $\langle \cdots\rangle_{\bm{\lambda}}$ an equilibrium average given control-parameter values $\bm{\lambda}$.

In the quasistatic (infinitely slow) limit, the probability distribution remains at equilibrium throughout the protocol, and the excess work approaches zero. For long-but-finite protocol duration, linear-response (LR) theory yields the leading-order contribution to the excess work [29],

$$\begin{equation} \langle W_\mathrm{ex}\rangle_{\bm{\Lambda}} \approx \int_{0}^{\Delta t}\mathrm{d}t\frac{\mathrm{d}{\bm{\lambda}}^\mathrm{T}}{\mathrm{d}t}\zeta[\bm{\lambda}(t)]\frac{\mathrm{d}{\bm{\lambda}}}{\mathrm{d}t}, \end{equation} \tag{ 3 }$$

in terms of the generalized friction tensor with elements

$$\begin{equation} \zeta_{j \ell}(\lambda) \equiv \beta \int_0^{\infty} \mathrm{d}t \, \langle \delta f_{\lambda_{j}}(t) \delta f_{\lambda_{\ell}}(0)\rangle_{\bm{\lambda}}. \end{equation} \tag{ 4 }$$

$\zeta_{j \ell}$ is the Hadamard product $\beta \langle\delta f_{\lambda_{j}} \delta f_{\lambda_{\ell}}\rangle_{\bm{\lambda}} \circ \tau_{j \ell}$ of the conjugate-force covariance (the force fluctuations) and the integral relaxation time

$$\begin{equation} \tau_{j \ell} \equiv \int_0^{\infty} \mathrm{d}t \, \frac{\langle \delta f_{\lambda_{j}}(t) \delta f_{\lambda_{\ell}}(0)\rangle_{\bm{\lambda}}}{\langle \delta f_{\lambda_{j}} \delta f_{\lambda_{\ell}}\rangle_{\bm{\lambda}}}, \end{equation} \tag{ 5 }$$

the characteristic time for these fluctuations to die out.

For overdamped dynamics, the friction can be calculated directly from the total energy as [35]

$$\begin{equation} \zeta_{j\ell}(\bm{\lambda}) = \int_{-\infty}^{\infty} \mathrm{d}x \, \frac{\partial_{\lambda_{j}}\Pi_\mathrm{eq}(x,\bm{\lambda})\partial_{\lambda_{\ell}}\Pi_\mathrm{eq}(x,\bm{\lambda})}{\pi_\mathrm{eq}(x,\bm{\lambda})} \, \end{equation} \tag{ 6 }$$

where $\Pi_\mathrm{eq}(x,\bm{\lambda}) \equiv \int_{-\infty}^{x}\mathrm{d}x'\pi_\mathrm{eq}(x',\bm{\lambda})$ is the equilibrium cumulative distribution function and $\partial_{\lambda_{j}}$ is the partial derivative with respect to $\lambda_{j}$.

Within the linear-response approximation, the excess work is minimized by a protocol with constant excess power [29]. For a single control parameter, this amounts to proceeding with velocity $\mathrm{d}\lambda^\mathrm{LR}/\mathrm{d}t \propto \zeta(\lambda)^{-1/2}$, which when normalized to complete the protocol in a fixed allotted time $\Delta t$, gives

$$\begin{equation} \frac{\mathrm{d}\lambda^\mathrm{LR} }{\mathrm{d}t}= \frac{\Delta \lambda}{\Delta t}\frac{\overline{{\zeta}^{1/2}}}{\sqrt{\zeta(\lambda)}}, \end{equation} \tag{ 7 }$$

where the overline denotes the spatial average over the naive (linear) path between the control-parameter endpoints.

For multidimensional control, the minimum-dissipation protocol solves the Euler-Lagrange equation

$$\begin{equation} \zeta_{j \ell}\frac{\mathrm{d}^2\lambda_{\ell}}{\mathrm{d}t^2}+ \frac{\partial\zeta_{j \ell}}{\partial \lambda_{m}} \frac{\mathrm{d}\lambda_{\ell}}{\mathrm{d}t}\frac{\mathrm{d}\lambda_{m}}{\mathrm{d}t} = \frac{1}{2}\frac{\partial\zeta_{\ell m}}{\partial \lambda_{j}} \frac{\mathrm{d}\lambda_{\ell}}{\mathrm{d}t}\frac{\mathrm{d}\lambda_{m}}{\mathrm{d}t}, \end{equation} \tag{ 8 }$$

where we have adopted the Einstein convention of implied summation over all repeated indices. We directly calculate the friction matrix from (6) and find geodesics by numerically solving (8) with specified initial and final control parameters, as described in refs. [36,37].

In the fast limit, the excess work approaches that of an instantaneous protocol, which spends no time relaxing towards equilibrium and requires excess work proportional (up to a factor of $k_\mathrm{B}T$) to the relative entropy $D(\pi_\mathrm{i}||\pi_\mathrm{f})$ between the initial and final equilibrium distributions [38]. Spending a short duration $\Delta t$ relaxing towards equilibrium throughout the protocol results in saved work $W_\text {save} \equiv k_\mathrm{B}T D(\pi_\mathrm{i}||\pi_\mathrm{f}) -W_\mathrm{ex}$ compared to an instantaneous protocol, which can be approximated as [38]

$$\begin{equation} \langle W_\mathrm{save} \rangle_{\bm{\Lambda}} \approx \int_{0}^{\Delta t}\mathrm{d}t \, \textbf{R}_{\bm{\lambda}_\mathrm{i}}^\mathrm{T}[\bm{\lambda}(t)] \, [\bm{\lambda}_{\mathrm{f}} - \bm{\lambda}(t)], \end{equation} \tag{ 9 }$$

in terms of the initial force-relaxation rate (IFRR)

$$\begin{equation} \textbf{R}_{\bm{\lambda}_\mathrm{i}}[\bm{\lambda}(t)] \equiv \frac{\mathrm{d}\langle \textbf{f}_{\bm{\lambda}} \rangle_{\bm{\lambda}_\mathrm{i}}}{\mathrm{d}t}\bigg|_{\bm{\lambda}(t)}, \end{equation} \tag{ 10 }$$

the rate of change of the initial mean conjugate forces at the current control-parameter values.

The saved work is maximized by the short-time efficient protocol (STEP) which spends the entire duration at the control-parameter value that maximizes the short-time power savings

$$\begin{equation} P_\text{save}^\mathrm{st}(\bm{\lambda}) \equiv \textbf{R}_{\bm{\lambda}_\mathrm{i}}^\mathrm{T}(\bm{\lambda})(\bm{\lambda}_{\mathrm{f}} - \bm{\lambda}), \end{equation} \tag{ 11 }$$

satisfying

$$\begin{equation} \frac{\partial P_\text{save}^\mathrm{st}(\bm{\lambda})}{\partial \bm{\lambda}}\bigg|_{\bm{\lambda}^\text{STEP}} = 0. \end{equation} \tag{ 12 }$$

The STEP achieves this by two instantaneous control-parameter jumps: one at the start from the initial value to the optimal value $\bm{\lambda}^{\mathrm{STEP}}$, and one at the end from $\bm{\lambda}^{\mathrm{STEP}}$ to the final value.

For overdamped dynamics, the short-time power saving from the STEP are conveniently expressed as (see the Supplementary Material Supplementarymaterial.pdf (SM), first section),

$$\begin{align} &P_\mathrm{save}^\mathrm{st}(\bm{\lambda}) = \nonumber\\ &\frac{1}{\gamma}\left\langle [f_{x}(x,\bm{\lambda}_\mathrm{f})-f_{x}(x,\bm{\lambda})][f_{x}(x,\bm{\lambda})-f_{x}(x,\bm{\lambda}_\mathrm{i})]\right\rangle_{\bm{\lambda}_\mathrm{i}}, \end{align} \tag{ 13 }$$

which is maximized if

$$\begin{equation} \left\langle \frac{\partial f_{x}(x,\bm{\lambda})}{\partial\bm{\lambda}}\left[f_{x}(x,\bm{\lambda})-\frac{f_{x}(x,\bm{\lambda}_\mathrm{i})+f_{x}(x,\bm{\lambda}_\mathrm{f})}{2}\right]\right\rangle_{\bm{\lambda}_\mathrm{i}} = 0. \end{equation} \tag{ 14 }$$

This is achieved by control-parameter values which for all x satisfy $\partial f_{x}(x,\bm{\lambda})/\partial\bm{\lambda} = \bm{0}$ or $f_{x}(x,\bm{\lambda}) = [f_{x}(x,\bm{\lambda}_\mathrm{i})+f_{x}(x,\bm{\lambda}_\mathrm{f})]/2$. In what follows we will enforce the second condition, although it is more stringent than (14) which only constrains a single average over the entire system distribution.

Here we consider a model system relevant to DNA-hairpin experiments: a Brownian bead driven by a time-dependent quadratic trapping potential with center and stiffness modulated by the focus and intensity of the laser. This model is also typical of steered molecular-dynamics simulations, which use a time-dependent quadratic potential to drive reactions [18]. The total potential $V_{\text{tot}} = V_{\mathrm{hp}} + V_{\text{trap}}$ is the sum of the static hairpin potential and time-dependent trap potential (shown schematically in fig. 1). The hairpin potential is modeled as a static double well (symmetric for simplicity) with the two minima at $x = 0$ and $x = \Delta x_\mathrm{m}$ representing the folded and unfolded states [5,6,33,34],

$$\begin{equation} V_\mathrm{hp}(x) = E_\mathrm{B} \left[\left(\frac{2x-\Delta x_\mathrm{m}}{\Delta x_\mathrm{m}}\right)^2-1\right]^{2}, \end{equation} \tag{ 15 }$$

for barrier height $E_\mathrm{B}$, distance $x_\mathrm{m}$ from the minimum to barrier, and distance $\Delta x_\mathrm{m}= 2x_\mathrm{m}$ between the minima. The system is driven by a quadratic trap

$$\begin{equation} V_\mathrm{trap}[x,x^\mathrm{c}(t),k(t)]= \frac{k(t)}{2}\left[x^\mathrm{c}(t)-x\right]^2, \end{equation} \tag{ 16 }$$

with time-dependent stiffness k(t) and center $x^\mathrm{c}(t)$.

The total work can be separated into two components, $W = W_{\rm c} + W_{\rm s}$, one for each control parameter. The trap-center component

$$\begin{equation} \langle W_{\rm c} \rangle_{\bm{\Lambda}} = \int_{0}^{\Delta t} \mathrm{d}t\, k(t) \left[ x^{\rm c}(t)-\langle x\rangle_{\bm{\Lambda}}\right]\frac{\mathrm{d}x^{\rm c}(t)}{\mathrm{d}t} \end{equation} \tag{ 17 }$$

is analogous to "force-distance" work. The stiffness component

$$\begin{equation} \langle W_{\rm s} \rangle_{\bm{\Lambda}} = \frac{1}{2}\int_{0}^{\Delta t} \mathrm{d}t \left\langle \left[x^{\rm c}(t)-x\right]^2\right\rangle_{\bm{\Lambda}} \frac{\mathrm{d}k(t)}{\mathrm{d}t} \end{equation} \tag{ 18 }$$

resembles "pressure-volume" work, i.e., the stiffness controls the effective volume available to the system, and the variance contributes to an effective pressure resisting changes in trap stiffness.

We consider protocols that drive the system between the two minima $x^{\rm c}_{\rm i} = 0$ and $x^{\rm c}_{\rm f} = \Delta x_{\rm m}$ with equal initial and final stiffness $(k_{\rm i} = k_{\rm f})$. This allows us to directly compare control over only trap center to control over both trap center and stiffness. We speculate that for unequal initial and final stiffness the qualitative features of the designed protocols would remain similar.

In the slow limit, the two-dimensional linear-response (2D LR) protocol minimizes dissipation by tightening the trap and slowing down as it traverses the barrier (fig. 2). Tightening the trap when approaching the barrier (fig. 2(b)) helps the system maintain roughly constant variance throughout the protocol and approximately linearly changes the quantiles of the position distribution (fig. 1), which is a generic property for minimum-dissipation protocols in optimal transport under full control [39–41]. When the trap does not tighten (e.g., naive and 1D LR protocols in fig. 2(b)), the variance increases as the system crosses the barrier and the quantiles do not change linearly. Slowing down while crossing the barrier (previously observed for one-dimensional (constant-stiffness) barrier crossing [34]) allows time for thermal fluctuations to kick the system over the barrier (fig. 2(a)).

In the SM, secound section we show that the amount that the trap tightens and slows down depends on the initial stiffness. We call the initial stiffness large (small) when the initial energy of the trap at the barrier $x_{\rm m}$ is significantly larger (smaller) than $E_{\rm b}$; i.e., $k_{\rm i} \gg E_{\rm b}/x_{\rm m}^2\, (k_{\rm i} \ll E_{\rm b}/x_{\rm m}^2)$. Physically, large (small) initial stiffness ensures that the initial equilibrium distribution is unimodal (bimodal). Throughout, we compare the stiffness to the scaled barrier height $E_{\rm B}/x_{\rm m}^2$, essentially comparing the initial energy of the trap and hairpin potential at the barrier. The 2D LR protocol changes stiffness most when the initial stiffness is comparable to the scaled barrier height $(k_{\rm i} \sim E_{\rm B}/x_{\rm m}^2)$, and leaves stiffness virtually unchanged when the initial stiffness is either large $(k_{\rm i} \gg E_{\rm B}/x_{\rm m}^2)$ or small $(k_{\rm i} \ll E_{\rm B}/x_{\rm m}^2)$.

The minimum-dissipation protocol in the fast limit (the STEP) maximizes the short-time power saving (13) by jumping from and to the control-parameter endpoints to spend the entire duration at control-parameter values $x^{\text{STEP}} = (x^{\rm c}_{\rm i}+x^{\rm c}_{\rm f})/2$, and $k^{\text{STEP}} = k_{\rm i}$. This result is independent of the hairpin potential since $f_{x}(x,x^\text{STEP},k^\text{STEP}) = [f_{x}(x,x^{\rm c}_{\rm i},k_{\rm i})+f_{x}(x,x^{\rm c}_{\rm f},k_{\rm i})]/2$ maximizes the short-time power saving (14) for all x independent of V(x) (note $k_{\rm f} = k_{\rm i}$).

Given theory describing minimum-dissipation control in both the slow and fast limits, we develop a simple interpolation scheme to design protocols that reduce dissipation at all driving speeds. Similar to the one-dimensional case [38], we choose the interpolated protocol to have an initial jump $(\lambda^\text{STEP}-\lambda_{\rm i})/(1+\Delta t/\tau)$ and a final jump $(\lambda_{\rm f}-\lambda^\text{STEP})/(1+\Delta t/\tau)$, and follow the original linear-response path between them,

$$\begin{equation} \lambda^{\rm interp}(t) = \frac{1}{1+\frac{\Delta t}{\tau}}\lambda^{\rm STEP} + \left(1-\frac{1}{1+\frac{\Delta t}{\tau}}\right)\lambda^{\rm LR}(t), \end{equation} \tag{ 19 }$$

with τ the crossover duration. This guarantees that the protocol approaches the minimum-dissipation protocol in both the fast and slow limits. For system timescale, we choose (primarily for its simplicity) the diffusion time $\tau_{\rm D} \equiv \Delta x_{\rm m}^2/(2D)$ between wells, for $D \equiv (\beta\gamma)^{-1}$. More sophisticated measures of relaxation time [37,39] could yield improved performance of the interpolated protocol.

For comparison to an ideal process, we evaluate the performance of optimal-transport (OT) theory under full control [40–42]. We use optimal transport to calculate the minimum work required —assuming complete control over the potential— to move probability from an initial to final distribution within a fixed duration as

$$\begin{equation} W_{\rm OT} = F_{\rm f} - F_{\rm i} + \frac{\beta\gamma}{\Delta t}\int_{0}^{1}\mathrm{d}y \left[Q_{\rm f}(y)-Q_{\rm i}(y)\right]^2, \end{equation} \tag{ 20 }$$

where $Q_{\rm f}$ and $Q_{\rm i}$ are the final and initial quantile functions (inverse cumulative distribution functions). We then perform a second optimization over the final probability distribution, subject to constrained initial and final control-parameter endpoints. This yields the minimum work to drive the trap between the two endpoints if we had full control over the potential (rather than just one- or two-dimensional parametric control).

Figure 2 compares naive, one-dimensional linear-response (1D LR), two-dimensional linear-response (2D LR), and interpolated protocols. We measure performance by the average work (the direct target of the protocol design) and the probability $p_{\bm{\Lambda}}(x<x_{\rm m})$ that the system remains in its initial well (related to the average flux) which was not directly considered in the design. The barrier height $E_{\rm B} = 4k_{\rm B} T$ is intermediate in the context of DNA hairpins [5,6], and the initial and final stiffness are comparable to the scaled barrier height, $k_{\rm i} = k_{\rm f} = E_{\rm B}/x_{\rm m}^2$. For a 1 μm bead in water at standard temperature and pressure, with inter-well distance $\Delta x_{\rm m} \approx 20\,\text{nm}$, the initial and final stiffness correspond to $k_{\rm i} = k_{\rm f} \approx 0.16\,\text{pN/nm}$, the 2D LR protocol reaches a maximum stiffness of $k_{\rm max} \approx 1.12\,\text{pN/nm}$, and the diffusion time between wells is $\tau_{\rm D} \approx 0.4\,\text{s}$.

For long duration $(\Delta t \gg \tau_{\rm D})$, the 1D LR protocol requires $(\sim1.6\times)$ less work than the naive; however, the 2D LR and interpolated protocols most significantly reduce work (fig. 2(c)); ${\sim}5.6\times$ less than naive, $\sim3.5\times$ less than 1D LR, and within 1% of full control). Intermediate-duration designed protocols give the largest-magnitude work reduction $\langle W \rangle_{\rm des} - \langle W \rangle_{\rm naive}$: 2D LR and interpolated protocols save $\sim2.7\,k_{\rm B} T$, whereas 1D LR only saves $\sim0.4\,k_{\rm B} T$ (fig. 2(e)).

1D LR protocols often reduce dissipation but as a side effect also decrease flux, as seen in fig. 2(d). 2D LR and interpolated protocols have the opposite effect, decreasing dissipation while increasing flux. For intermediate duration, the 2D LR protocol drives up to 78% and the interpolated up to 17% more probability to the destination well, compared to naive; the 1D LR drives 19% less.

For long duration $(\Delta t \gg \tau_{\rm D})$, two-dimensional control provides significant advantages over one-dimensional control for both average work and flux; however, for short duration the 2D LR protocol can perform worse than 1D LR and naive (figs. 2(c) and (e); similar behavior has been observed for multidimensional control of the Ising model [37]). For short duration, the system cannot keep up with the rapid changes in the trap potential, and the linear-response approximation breaks down. Although the increased stiffness of the 2D LR protocol results in the strongest driving and hence the greatest flux of the protocols considered here (fig. 2), it does so at the cost of increased dissipation for short duration. Indeed, the minimum-dissipation protocol for short duration ($\Delta t \ll \tau_{\rm D}$, the STEP) is monotonic and discrete. Our interpolated protocol asymptotes to the STEP in the short-duration limit, resulting in reduced dissipation and increased flux at any duration. In terms of dissipation, the interpolated protocol achieves within 1% of the minimum work under full control for short and long duration and remains within 30% of full control at intermediate duration (fig. 2(c)).

We have shown that multidimensional control protocols can significantly outperform their one-dimensional counterparts, improving both work and flux. For a system undergoing driven barrier crossing, one-dimensional control of only the trap center limits the control over the position distribution, with a large increase in variance as the protocol crosses the barrier (fig. 1). Control over both the trap center and stiffness makes approximately linear driving of the position mean and variance between specified endpoints possible, consistent with optimal-transport protocols that minimize work under full control [40–42]. This significantly reduces the work required to drive the system between the two wells and increases the flux compared to naive and one-dimensional control protocols (fig. 2). The main shortcoming of the multidimensional linear-response protocols is that they can perform worse than naive for short duration; however, we remedy this issue by combining linear-response and STEP frameworks to give interpolated protocols that reduce dissipation at any duration. For the model system and parameters we explored, the largest reduction in dissipation occurs from one- to two-dimensional control, and the dissipation in the two-dimensional interpolated protocol is within 30% of full control for intermediate duration and within 1% for short and long duration.

The model system closely resembles DNA-hairpin experiments, and we explore experimentally relevant parameters [5,6]. Our results reveal general design principles for driven barrier crossing that can be readily implemented experimentally: the designed protocols 1) slow down and tighten the trap as it crosses the energy barrier, thereby driving the mean position between the two wells at constant rate while maintaining constant variance; and 2) jump at the beginning and end of the protocol, with larger jumps for faster protocols. Recent experimental protocols implement one-dimensional control and demonstrate significant work reductions from designed protocols [32]. We show that adding an additional control parameter (trap stiffness) can dramatically improve the performance over the one-dimensional counterpart (up to $3.5{}\times$ less work and 80% increased probability of reaching the target well). Although multidimensional control is more difficult to implement, the performance gains can be significant.

This work is supported by an SFU Graduate Deans Entrance Scholarship (SB), an NSERC Discovery Grant and Discovery Accelerator Supplement (DAS), and a Tier-II Canada Research Chair (DAS), and was enabled in part by support provided by WestGrid (www.westgrid.ca) and Compute Canada Calcul Canada (www.computecanada.ca). The authors thank Miranda Louwerse (SFU Chemistry) for sharing code used to determine geodesics of the friction metric and enlightening feedback on the manuscript.

Data availability statement: The data that support the findings of this study are available upon reasonable request from the authors.