Optimal discrete control: minimizing dissipation in discretely driven nonequilibrium systems

TL;DR

This paper develops a theoretical framework to quantify and minimize energy dissipation in systems driven by discrete, instantaneous perturbations, with applications to molecular machines and nonequilibrium systems.

Contribution

It generalizes existing results to include finite-time, discrete control protocols and provides a method to design minimum-dissipation control strategies.

Findings

Derived the energetic cost of discrete, instantaneous driving out of equilibrium.

Compared theoretical formalism with an exactly solvable model system.

Demonstrated dissipation reduction in a multistable molecular machine model.

Abstract

Microscopic machines utilize free energy to create and maintain out-of-equilibrium organization in virtually all living things. Often this takes the form of converting the free energy stored in nonequilibrium chemical potential differences into useful work, via a series of reactions involving the binding, chemical catalysis, and unbinding of small molecules. Such chemical reactions occur on timescales much faster than the protein conformational rearrangements they induce. Here, we derive the energetic cost for driving a system out of equilibrium via a series of such effectively instantaneous (and hence discrete) perturbations. This analysis significantly generalizes previously established results, and provides insight into qualitative, as well as quantitative, aspects of finite-time, minimum-dissipation discrete control protocols. We compare our theoretical formalism to an exactly solvable model system and also demonstrate the dissipation reduction achievable in a simple multistable model for a discretely driven molecular machine.

Equations

$$\pi (x|\mathbf{\lambda })=e^{-\beta E(x,\mathbf{\lambda })+\beta F(\mathbf{\lambda })},$$

$$\Lambda \equiv \{{\lambda }_i,\Delta t_i\}.$$

$$⟨\beta W{⟩}_{{\mathbf{\lambda }}_i\to {\mathbf{\lambda }}_{i+1}}=\beta \int \left[E(x,{\mathbf{\lambda }}_{i+1})-E(x,{\mathbf{\lambda }}_i)\right]p_{\Lambda }(x,t_{i,i+1})\mathrm{d}x,$$

$$\tau \equiv \sum _{i=1}^{N-1}\Delta t_i,$$

$$\Delta {\mathbf{\lambda }}_{\mathrm{tot}}\equiv \sum _{i=0}^{N-1}\Delta {\mathbf{\lambda }}_{i,i+1}={\mathbf{\lambda }}_N-{\mathbf{\lambda }}_0.$$

$$⟨\beta W{⟩}_{{\mathbf{\lambda }}_0\to {\mathbf{\lambda }}_1}=\beta \int \left[E(x,{\mathbf{\lambda }}_1)-E(x,{\mathbf{\lambda }}_0)\right]\pi (x|{\mathbf{\lambda }}_0)\mathrm{d}x.$$

$$⟨\beta W{⟩}_{{\mathbf{\lambda }}_0\to {\mathbf{\lambda }}_1}=\int \ln \left[\frac{\pi (x|{\mathbf{\lambda }}_0)}{\pi (x|{\mathbf{\lambda }}_1)}\right]\pi (x|{\mathbf{\lambda }}_0)\mathrm{d}x+\beta \Delta F_{0,1},$$

$$⟨\beta W{⟩}_{{\mathbf{\lambda }}_0\to {\mathbf{\lambda }}_1}=D[\pi (x|{\mathbf{\lambda }}_0)\Vert \pi (x|{\mathbf{\lambda }}_1)]+\beta \Delta F_{0,1}.$$

$$⟨\beta W{⟩}_{\Lambda }=\sum _{i=0}^{N-1}D\left[\pi (x||{\mathbf{\lambda }}_i)\Vert \pi (x|{\mathbf{\lambda }}_{i+1})\right]+\beta \Delta F_{\mathrm{tot}},$$

$$⟨\beta W_{\mathrm{ex}}{⟩}_{\Lambda }$$

$$=\sum _{i=0}^{N-1}D[\pi (x|{\mathbf{\lambda }}_i)\Vert \pi (x|{\mathbf{\lambda }}_{i+1})].$$

$$D\left[\pi (x|{\mathbf{\lambda }}_i)\Vert \pi (x|{\mathbf{\lambda }}_{i+1})\right]\approx \frac{1}{2}{\beta }^2⟨\delta f_j\delta f_k{⟩}_{{\mathbf{\lambda }}_i}\Delta {\lambda }_{i,i+1}^j\Delta {\lambda }_{i,i+1}^k,$$

$$⟨\beta W_{\mathrm{ex}}{⟩}_{\Lambda }\approx \frac{1}{2}{\beta }^2\sum _{i=0}^{N-1}⟨\delta f_j\delta f_k{⟩}_{{\mathbf{\lambda }}_i}\Delta {\lambda }_{i,i+1}^j\Delta {\lambda }_{i,i+1}^k,$$

$$E(x,{\mathbf{\lambda }}_i)$$

$$=E(x,{\mathbf{\lambda }}_{i-1})-\mathbf{f}{|}_{{\mathbf{\lambda }}_{i-1}}\cdot \Delta {\mathbf{\lambda }}_{i-1,i}$$

$$⟨f_j(\Delta t){⟩}_{{\mathbf{\lambda }}_{i-1},{\mathbf{\lambda }}_i}=⟨f_j{⟩}_{{\mathbf{\lambda }}_i}+\beta ⟨\delta f_j(0)\delta f_k(\Delta t){⟩}_{{\mathbf{\lambda }}_i}\Delta {\lambda }_{i-1,i}^k,$$

$$⟨\beta W{⟩}_{{\mathbf{\lambda }}_i\to {\mathbf{\lambda }}_{i+1}}$$

$$\approx \beta \Delta {\lambda }_{i,i+1}^j⟨f_j{⟩}_{{\mathbf{\lambda }}_i}+{\beta }^2\Delta {\lambda }_{i,i+1}^j⟨\delta f_j(0)\delta f_k(\Delta t_i){⟩}_{{\mathbf{\lambda }}_i}\Delta {\lambda }_{i-1,i}^k$$

$$=⟨\beta W^{\infty }{⟩}_{{\mathbf{\lambda }}_i\to {\mathbf{\lambda }}_{i+1}}+⟨\beta W_{\mathrm{ex}}^{\mathrm{neq}}{⟩}_{{\mathbf{\lambda }}_i\to {\mathbf{\lambda }}_{i+1}}$$

$$⟨\beta$$

$$\approx ⟨\beta W_{\mathrm{ex}}^{\infty }{⟩}_{\Lambda }+{\beta }^2\sum _{i=1}^{N-1}\Delta {\lambda }_{i,i+1}^j⟨\delta f_j(0)\delta f_k(\Delta t_i){⟩}_{{\mathbf{\lambda }}_i}\Delta {\lambda }_{i-1,i}^k$$

$$\approx {\beta }^2\sum _{i=0}^{N-1}\frac{1}{2}⟨\delta f_j\delta f_k{⟩}_{{\mathbf{\lambda }}_i}\Delta {\lambda }_{i,i+1}^j\Delta {\lambda }_{i,i+1}^k+{\beta }^2\sum _{i=1}^{N-1}\Delta {\lambda }_{i,i+1}^j⟨\delta f_j(0)\delta f_k(\Delta t_i){⟩}_{{\mathbf{\lambda }}_i}\Delta {\lambda }_{i-1,i}^k,$$

$$⟨\beta W_{\mathrm{ex}}{⟩}_{\Lambda }={\beta }^2\sum _{i=0}^{N-1}⟨\delta f_j\delta f_k{⟩}_{{\mathbf{\lambda }}_i}\Delta {\lambda }_{i,i+1}^j\Delta {\lambda }_{i,i+1}^k\left[\frac{1}{2}+\frac{⟨\delta f_j(0)\delta f_k(\Delta t_i){⟩}_{{\mathbf{\lambda }}_i}}{⟨\delta f_j\delta f_k{⟩}_{{\mathbf{\lambda }}_i}}\frac{\Delta {\lambda }_{i-1,i}^k}{\Delta {\lambda }_{i,i+1}^k}\right].$$

$$⟨\beta W_{\mathrm{ex}}{⟩}_{\Lambda }$$

$$\equiv \sum _{i=0}^{N-1}{\Gamma }_{jk}({\mathbf{\lambda }}_{i-1},{\mathbf{\lambda }}_i,{\mathbf{\lambda }}_{i+1},\Delta t_i)\Delta {\lambda }_{i,i+1}^j\Delta {\lambda }_{i,i+1}^k$$

$$\equiv \sum _{i=0}^{N-1}{\mathcal{D}}_i^2,$$

$${\Gamma }_{jk}({\mathbf{\lambda }}_{i-1},{\mathbf{\lambda }}_i,{\mathbf{\lambda }}_{i+1},\Delta t_i)\equiv {\beta }^2⟨\delta f_j\delta f_k{⟩}_{{\mathbf{\lambda }}_i}\left[\frac{1}{2}+\frac{⟨\delta f_j(0)\delta f_k(\Delta t_i){⟩}_{{\mathbf{\lambda }}_i}}{⟨\delta f_j\delta f_k{⟩}_{{\mathbf{\lambda }}_i}}\frac{\Delta {\lambda }_{i-1,i}^k}{\Delta {\lambda }_{i,i+1}^k}\right],$$

$$⟨\beta W_{\mathrm{ex}}{⟩}_{\Lambda }\ge \frac{1}{N}{\left(\sum _{i=0}^N{\mathcal{D}}_i\right)}^2.$$

$${\mathcal{D}}_i=\mathcal{D}.$$

$$\Delta {\lambda }_{i,i+1}^{\ast }\propto \frac{1}{\sqrt{\Gamma ({\lambda }_{i-1}^{\ast },{\lambda }_i^{\ast },{\lambda }_{i+1}^{\ast },\Delta t_i)}},$$

$$\frac{⟨\delta f(0)\delta f(\Delta t_i){⟩}_{{\lambda }_i}}{⟨\delta f^2{⟩}_{{\lambda }_i}}=e^{-\Delta t_i/{\tau }_{\mathrm{R}}({\lambda }_i)}.$$

$${\Delta t_i}^{\ast }={\tau }_{\mathrm{R}}({\lambda }_i)\left[\frac{\tau }{{\sum }_{n=1}^{N-1}{\tau }_{\mathrm{R}}({\lambda }_n)}-\frac{{\sum }_{n=1}^{N-1}{\tau }_{\mathrm{R}}({\lambda }_n)\ln \left({\mathcal{P}}_n/{\mathcal{P}}_i\right)}{{\sum }_{n=1}^{N-1}{\tau }_{\mathrm{R}}({\lambda }_n)}\right],$$

$${\Delta t_i}^{\ast }\propto {\tau }_{\mathrm{R}}({\lambda }_i).$$

$$E_{\mathrm{t}}(x,{\lambda }_i)=\frac{1}{2}{k}_{\mathrm{t}}(x-{\lambda }_i{)}^2.$$

Full text

,

1 Introduction

At all scales, biological systems exhibit a striking degree of organization and coordination. The continuous flow of information, energy, and material within and between biological cells preserves the ordered structure necessary for their proper functioning. At microscopic scales, a variety of molecular machines are largely responsible for maintaining this order, performing a broad range of intracellular tasks [1]. For instance, the rotary ${\rm F}_{\rm o}{\rm F}_{1}$ ATP-synthase motor produces the cellular energy currency ATP [2], whereas transport motors such as kinesin, myosin, and dynein are responsible for the directed trafficking of material within the cell [3, 4].

Operationally, molecular machines function by coupling the free energy stored in nonequilibrium environmental conditions (often imbalances of chemical potential) to mechanical motion. The ${\rm F}_{1}$ ATPase motor, for instance, can operate at speeds up to $\sim$350 rotations per second [5], presumably far from thermodynamic equilibrium, by coupling rotational torque of a central crankshaft to the hydrolysis or synthesis of ATP.

A variety of single-molecule experimental techniques have made possible the detailed study of these molecular machines [6, 2, 5, 7, 8], but their out-of-equilibrium operation complicates a theoretical understanding of energetic flows into, within, and out of these systems. Under the hypothesis that selective pressures favor efficient cellular machinery [9, 10], a theory which elucidates how energetic flows depend on operational parameters promises to deepen our understanding of the fundamental operational constraints facing evolved molecular machines. In practice, this would aid in the creation of de novo molecular machines, perhaps accelerating their implementation for next-generation nanomedicine [11, 12]

Much previous work on the properties of minimum-dissipation control (applied to a variety of model systems including the erasure of a classical bit [13] and the reversal of magnetization in an Ising magnet [14], among others [15, 16, 17, 18, 19, 20]) has assumed that the system of interest is subjected to a controlling apparatus that can be manipulated in a continuous manner [21]. While this applies well to the single-molecule experimental paradigm, many interesting microscopic systems, such as biomolecular machines, often drive their mechanical motion via a sequence of chemical reactions. The time scales of chemical reaction and mechanical response can differ by several orders of magnitude; as such, the driving process is well approximated by a series of stochastically timed discrete perturbations to a thermodynamic system, as opposed to a continuous driving process.

Linear-response theory has been a vital tool in the investigation of nonequilibrium relaxation, from relaxation kinetics [22] to chemical dynamics [23, 24, 25] to nonequilibrium dissipation [26, 27, 28]. As a first-order theory, we sacrifice applicability far from the considered limits in order to gain tractability and generality. Moreover, linear-response theory has proven surprisingly effective in describing dissipation, even in contexts that clearly extend beyond its regime of rigorous validity [29].

In this article, as a first step toward adapting the control-protocol framework to chemically driven systems, and to isolate the effect of discrete control parameter changes, we develop a theoretical framework for nonequilibrium control using deterministic discrete control parameter changes. Our central result (10mopq) quantifies the nonequilibrium energetic costs associated with discretely driving a microscopic system. In particular, we assume the system is within the linear-response regime and is subject to sufficiently weak perturbations that a low-order approximation of the energetic costs of discrete steps is warranted. Within these limitations, this framework allows for straightforward optimization of discrete driving protocols that accounts for both the effects of the size of discrete steps as well as the local relaxation times, which leads to novel characteristics of discrete protocols not observed in continuously driven systems (Fig. 5). This work complements and generalizes previous results on the entropy production associated with discrete processes [30]. In the continuous-driving limit, our formalism reduces to previously known results, namely the thermodynamic-length formalism introduced by Crooks [31] and the entropy-differential metric of Burbea and Rao [32], and is related to the generalized friction coefficient of Sivak and Crooks [21].

The paper proceeds as follows: §2 introduces the relevant theoretical background information; §3 derives an exact expression for the energetic cost of a series of discrete steps applied to an equilibrium system, and then investigates a small-perturbation approximation of this expression to compare to established results; §4 extends this analysis to systems which are out of equilibrium, quantifying within the linear-response regime the excess energetic cost due to incomplete relaxation; §5 discusses general implications of optimal, minimum-work protocols in discretely driven systems, under the assumptions laid out in the previous sections; finally, §6 and §7 explore the quantitative implications of this theoretical framework in two model systems: a translating harmonic trap and a periodic potential.

2 Background

We consider a system in contact with a heat bath, subject to a set of experimentally controlled parameters $\bm{\lambda}$, such that at equilibrium the distribution over microstates is

[TABLE]

where $\beta\equiv(k_{\rm B}T)^{-1}$ is the inverse temperature, $E(x,\bm{\lambda})$ is the system Hamiltonian, and $F(\bm{\lambda})$ is the equilibrium free energy at control parameter vector $\bm{\lambda}$.

A control protocol $\Lambda:\bm{\lambda}_{0}\to\bm{\lambda}_{N}$ is a particular time-dependent perturbation applied to the control parameter vector $\bm{\lambda}$ to transform it between an initial $\bm{\lambda}_{0}$ and final $\bm{\lambda}_{N}$ in a prescribed time $\tau$. For a given control protocol, the system responds stochastically. Across the entire control protocol $\Lambda$, the average amount of excess work (supplied by an external source)—or work required above and beyond the equilibrium free energy difference $\beta\Delta F_{\rm tot}$ between the initial and final control parameter values—is $\langle\beta W_{\rm ex}\rangle_{\Lambda}\equiv\langle\beta W\rangle_{\Lambda}-\beta\Delta F_{\rm tot}$, where $\langle\cdots\rangle_{\Lambda}$ indicates an average of system responses to the control protocol $\Lambda$.

Here, we consider discrete control protocols which consist of a series of instantaneous perturbations $\Delta\bm{\lambda}_{i,i+1}\equiv\bm{\lambda}_{i+1}-\bm{\lambda}_{i}$ for consecutive control parameter values $\bm{\lambda}_{i}$ and $\bm{\lambda}_{i+1}$. The system spends a prescribed time ${\Delta t_{i}}$ at each control parameter $\bm{\lambda}_{i}$. Thus, each protocol is defined by a set of control parameter values and the associated times spent at them:

[TABLE]

Previous work has considered energetic flows in discrete-stepping processes [30, 33, 34, 35], but these efforts typically focused on the continuum limit. More recent investigations of driven nonequilibrium systems [21, 15, 18, 13, 19, 20] have focused on control protocols which are continuous functions of time. In contrast to these previous works, and motivated by the chemically driven paradigm characteristic of microscopic machines, we consider a control protocol as a series of discrete steps of substantial size.

The average work (divided by $k_{\rm B}T$) associated with a particular discrete perturbation that transforms the control parameter from $\bm{\lambda}_{i}$ to $\bm{\lambda}_{i+1}$ during the protocol $\Lambda$ is

[TABLE]

where $p_{\Lambda}(x,t_{i,i+1})$ is the (generally nonequilibrium) distribution over system microstates $x$ at the time $t_{i,i+1}$ that the control parameter changes from $\bm{\lambda}_{i}$ to $\bm{\lambda}_{i+1}$ during the protocol $\Lambda$, and angle brackets $\langle\cdots\rangle_{\bm{\lambda}_{i}\to\bm{\lambda}_{i+1}}$ indicate a nonequilibrium average as the control parameter changes from $\bm{\lambda}_{i}$ to $\bm{\lambda}_{i+1}$ during protocol $\Lambda$.

We consider protocols which start in equilibrium at initial control parameter value $\bm{\lambda}_{0}$, equivalent to taking the time spent at $\bm{\lambda}_{0}$ to infinity, ${\Delta t_{0}}\to\infty$. A particular control protocol begins with the first control parameter change $\bm{\lambda}_{0}\to\bm{\lambda}_{1}$, and finishes when the control parameter $\bm{\lambda}$ arrives at its terminal value $\bm{\lambda}_{N}$. The protocol duration $\tau$ is the time taken to complete a given protocol, not counting the time taken to equilibrate at the initial control parameter value. Thus, for a protocol with $N+1$ control parameter values $\bm{\lambda}_{0},\bm{\lambda}_{1},\ldots,\bm{\lambda}_{N}$, the protocol duration is

[TABLE]

and the total control parameter displacement for fixed control parameter endpoints $\bm{\lambda}_{0},\bm{\lambda}_{N}$ is

[TABLE]

For given control parameter endpoints $\bm{\lambda}_{0}$ and $\bm{\lambda}_{N}$ (hence given protocol displacement $\Delta\bm{\lambda}_{\rm tot}$) and duration $\tau$, minimizing work involves choosing intermediate control parameter values $\bm{\lambda}_{i}$ and associated dwell times ${\Delta t_{i}}$.

3 Infinite-time work

We first consider the work associated with making a single discrete change to the control parameter of a system which is initially at equilibrium with the control parameter $\bm{\lambda}_{0}$ (1). The average work required to discretely change the control parameter vector from $\bm{\lambda}_{0}$ to $\bm{\lambda}_{1}$ is (3), with the equilibrium initial distribution $p_{\Lambda}(x,t_{0,1})=\pi(x|\bm{\lambda}_{0})$ (1):

[TABLE]

From the definition of the equilibrium ensemble at a fixed control parameter (1), the energy can be expressed in terms of the equilibrium distribution and free energy, $\beta E(x,\bm{\lambda}_{i})=-\ln\pi(x|\bm{\lambda}_{i})+\beta F(\bm{\lambda}_{i})$. The average work from (6) can then be written solely in terms of equilibrium distributions,

[TABLE]

where $\beta\Delta F_{0,1}=\beta F(\bm{\lambda}_{1})-\beta F(\bm{\lambda}_{0})$ is the difference between the equilibrium free energy at control parameter $\bm{\lambda}_{0}$ and $\bm{\lambda}_{1}$. The integral in (7) is the relative entropy (or Kullback-Leibler divergence) $D[p(x)||q(x)]\equiv\int\ln[p(x)/q(x)]\,p(x){\rm d}x$ [36] between the equilibrium distributions before ($p(x)=\pi(x|\bm{\lambda}_{0})$) and after ($q(x)=\pi(x|\bm{\lambda}_{1})$) the control parameter change,

[TABLE]

A protocol consists of $N$ such control parameter steps $\Delta\bm{\lambda}_{i,i+1}$, for $i=0,\ldots,N-1$. If at each control parameter value $\bm{\lambda}_{i}$ the system fully equilibrates, then the average work associated with any step $\Delta\bm{\lambda}_{i,i+1}$ is of the same form (8), and the work to complete the entire protocol is the sum of the work associated with each individual step,

[TABLE]

where $\beta\Delta F_{\rm tot}\equiv\beta F(\bm{\lambda}_{N})-\beta F(\bm{\lambda}_{0})=\sum_{i=0}^{N-1}\beta\Delta F_{i,i+1}$ is the equilibrium free energy change between the initial and final control parameter values. Thus the average excess work is

[TABLE]

For sufficiently small control parameter steps $\Delta\bm{\lambda}_{i,i+1}$, the relative entropy in (10b) can be Taylor expanded about its current value $\bm{\lambda}_{i}$ to yield

[TABLE]

where $\langle\delta f_{j}\delta f_{k}\rangle_{\bm{\lambda}_{i}}$ is the equilibrium covariance of conjugate forces $f_{j}\equiv-\partial_{\lambda^{j}}E$ at control parameter $\bm{\lambda}_{i}$ (see A for details) [21]. Throughout, we employ the Einstein summation notation, where repeated indices are implicitly summed over.

Substituting (10k) in (10b) gives the average excess work required to perform the infinite-time discrete protocol $\Lambda$,

[TABLE]

which in the continuous-protocol limit is equivalent to the Burbea-Rao entropy differential metric [32] and the thermodynamic metric derived by Crooks [31].

4 Nonequilibrium excess work

To consider the more general situation of finite-time protocols, where at each control parameter value the system does not fully equilibrate, we appeal to static linear-response theory [37]. For a system at equilibrium for control parameter $\bm{\lambda}_{i-1}$, the energy at the next control parameter value $\bm{\lambda}_{i}$ in the protocol $\Lambda$ can be linearly approximated as

[TABLE]

where $\bm{f}|_{\bm{\lambda}_{i-1}}$ is the vector of conjugate forces with elements $f_{j}$, evaluated at $\bm{\lambda}_{i-1}$. When the control parameter instantaneously changes from $\bm{\lambda}_{i-1}$ to $\bm{\lambda}_{i}$, the time-dependent relaxation of $f_{j}$ towards its equilibrium value at $\bm{\lambda}_{i}$ is, under the linear-response approximation,

[TABLE]

where $\langle f_{j}({\Delta t})\rangle_{\bm{\lambda}_{i-1},\bm{\lambda}_{i}}$ indicates an average over the instantaneous nonequilibrium system distribution after relaxing (under the Hamiltonian $E(x,\bm{\lambda}_{i})$) for a time ${\Delta t}$ starting from the equilibrium distribution at $\bm{\lambda}_{i-1}$. Both the force autocovariance and the average force on the RHS are taken over the equilibrium ensemble at fixed control parameter value $\bm{\lambda}_{i}$. B provides a detailed derivation of (10mn). The second RHS term in (10mn) is a linear-response relaxation function, which have a long history of usage to understand the dynamic properties of nonequilibrium systems, for instance in chemical relaxation kinetics [22] and chemical dynamics [25, 24, 23].

If the system relaxes for a time ${\Delta t_{i}}$ at control parameter $\bm{\lambda}_{i}$ before the next control parameter change $\bm{\lambda}_{i}\to\bm{\lambda}_{i+1}$, this step requires average work (3)

[TABLE]

for infinite-time work $\langle\beta W^{\infty}\rangle=\beta\Delta\lambda_{i,i+1}^{j}\langle f_{j}\rangle_{\bm{\lambda}_{i}}$ from §3 and linear-response correction $\langle\beta W_{\rm ex}^{\rm neq}\rangle$ due to incomplete system relaxation.

(10moc) is only strictly valid if the system was at equilibrium for $\bm{\lambda}_{i-1}$ before the step to $\bm{\lambda}_{i}$, so more generally the work required for the control parameter change $\bm{\lambda}_{i}\to\bm{\lambda}_{i+1}$ includes contributions from all previous steps. However, (10moc) approximates the work when the force autocovariance associated with the most recent step is the largest time-dependent contribution to the excess work. This limit is reached when the time spent at each control parameter value is long compared to the relaxation time of the conjugate forces. As discussed in detail in C, this approximation is fundamentally distinct from approximations made in deriving the continuous-protocol formalism [21].

Within this limit, the total average excess work required to perform a discrete control protocol $\Lambda$ is

[TABLE]

where (10mopc) uses the infinite-time excess work for small steps (10l). The assumption that the system begins in equilibrium at $\bm{\lambda}_{0}$ ensures that the $i=0$ term does not contribute to the total nonequilibrium excess work. By setting $\bm{\lambda}_{-1}=\bm{\lambda}_{0}$, both sums in (10mop) can be taken over the same index range, leading to a more compact form for the excess work:

[TABLE]

This captures the combined effects of the control parameter step sizes $\Delta\bm{\lambda}_{i,i+1}$ and time allocations ${\Delta t_{i}}$ on the excess work during a discrete control protocol $\Lambda$. The time-independent, infinite-time contribution penalizes large control parameter steps departing from regions with large force covariance. The time-dependent linear-response correction penalizes steps that are particularly quick (reflected by the force autocorrelation factor) and/or large (reflected by the $\Delta\lambda_{i,i+1}^{j}\Delta\lambda_{i-1,i}^{k}$ step-size factor).

Equation (10mopq) generalizes the near-equilibrium expression for the dissipation of a discrete control protocol in [30], because the time dependence captured by the conjugate-force autocovariance in our approach allows for non-exponential relaxation kinetics, and the explicit form permits simultaneous optimization of both the placement of control parameter values $\bm{\lambda}_{i}$ as well as the allocation of times ${\Delta t_{i}}$. Moreover, C provides an alternative derivation of the nonequilibrium excess work contribution (10moc) using dynamic linear response theory to show that in the continuous-protocol limit, a linear-response correction to the excess work can recover the generalized friction formalism from [21], but only if approximations are made which are incompatible with those used to derive (10mopq).

5 Minimum-work protocols

The nonequilibrium excess work (10mopq) provides a relatively simple expression, within the linear-response approximation, for the energetic cost required to perform discrete control protocol $\Lambda$. Although the specific form of a minimum-work protocol depends on the particular system, there are two special cases admitting analytic solutions: the case of the infinite-time limit (§3) where the time-dependent term in (10mopq) is negligible, and the case where there is a single dominant exponential relaxation mode.

In general, the excess work can be approximated as

[TABLE]

and

[TABLE]

where, in this case, the $k$ index within the brackets is not summed over. Interpretation is most immediate in the continuous-protocol limit, where each $\mathcal{D}_{i}$ is the distance along an infinitesimal segment ${\rm d}\bm{\lambda}$ of the control protocol $\Lambda$, measured with respect to the metric $\Gamma_{jk}(\bm{\lambda}_{i-1},\bm{\lambda}_{i},\bm{\lambda}_{i+1},{\Delta t_{i}})$; therefore, the sum $\sum_{i=0}^{N}\mathcal{D}_{i}$ over all steps gives the thermodynamic length between the initial and final equilibrium system macrostates [31, 30].

For a positive semidefinite force-autocovariance matrix, the total excess work of a particular control protocol can be lower-bounded via the Cauchy-Schwarz inequality:

[TABLE]

The lower bound is saturated if and only if the $\mathcal{D}_{i}$ are identical,

[TABLE]

Along an optimal protocol (indicated by the superscript $*$), the condition (10moprt) implies that $\langle\beta W_{\rm ex}^{*}\rangle_{\bm{\lambda}_{i}\rightarrow\bm{\lambda}_{i+1}}=\mathcal{D}^{2}$, and thus equal excess work is done during each step of the protocol.

For a single control parameter with fixed endpoints $\lambda_{0},\lambda_{N}$ and a given set of time allocations ${\Delta t_{i}}$, the condition (10moprt) implies the optimal placement of control parameter values through the proportionality

[TABLE]

but the implicit dependence of $\Gamma$ on the step size through $\lambda_{i-1}^{*},\lambda_{i}^{*}$, and $\lambda_{i+1}^{*}$ complicates the practical use of this bound for deriving optimal protocols. However, the proportionality (10mopru) can give useful qualitative guidance into the general properties of protocols which saturate the lower bound (10moprs). In particular, optimal control parameter placement tends to avoid regions with large force variance and slowly decaying force autocovariance, subject to the quadratic cost $\Delta\lambda_{i,i+1}^{j}\Delta\lambda_{i,i+1}^{k}$ on step sizes. For more than one control parameter, the qualitative insights gained from the lower bound (10moprs) and the equality (10moprt) can provide a way to derive the optimal time-schedule along a particular path in control parameter space, but unfortunately they do not generally provide a constructive means to identify a path that saturates the bound.

In the infinite-time limit, where $\Gamma_{jk}(\bm{\lambda}_{i-1},\bm{\lambda}_{i},\bm{\lambda}_{i+1},{\Delta t_{i}})\to\Gamma_{jk}^{(\infty)}(\bm{\lambda}_{i})=\frac{1}{2}\beta^{2}\langle\delta f_{j}\delta f_{k}\rangle_{\bm{\lambda}_{i}}$, our predictions reduce to previous calculations by Nulton et al. [30] of the optimal placement of discrete steps. In particular, for a single control parameter, the condition (10moprt) implies that optimal protocols have the proportionality $\Delta\lambda_{i,i+1}^{*}\propto 1/\sqrt{\langle\delta f^{2}\rangle_{\lambda_{i}^{*}}}$. Furthermore, in the continuous-protocol limit, the infinite-time thermodynamic length between the initial and final control parameters (measured with respect to $\Gamma_{jk}^{(\infty)}(\bm{\lambda}_{i})$) converges to that of Crooks [31].

The optimal time allocation is analytically solvable for a single control parameter when the time dependence in (10mopq) takes a simple, control-parameter-dependent, exponential form:

[TABLE]

Using Lagrange multipliers, the optimal allocation of time among a fixed set of control parameter values, subject to the protocol duration constraint (4), is

[TABLE]

where $\mathcal{P}_{n}\equiv\beta^{2}\Delta\lambda_{n,n+1}\Delta\lambda_{n-1,n}\langle\delta f^{2}\rangle_{\lambda_{n}}/\tau_{\rm R}(\lambda_{n})$ (see E for a detailed derivation).

In the long-duration limit, where $\tau\gg\sum_{s=1}^{N-1}\tau_{\rm R}(\lambda_{s})$, the second RHS term in (10moprw) is negligible, and the optimal allocation of time takes on the simple form

[TABLE]

Intuitively, this implies that along minimum-work protocols, more time is allocated to regions where the integral relaxation time [38] is larger.

The special case of a single control parameter and exponential relaxation kinetics produces (10moprw), which recovers the result of Nulton et al. in [30]. However, our more general framework (10mopq) can be applied to a broader class of problems (as we detail in §6 and §7), in particular to cases with multiple control parameters, non-exponential relaxation kinetics, and optimization of control parameter placements. For more general scenarios, even with one control parameter, analytic optimization methods become cumbersome, and no simple analogs of (10moprw) can be found. Nevertheless, (10mopq) provides a relatively simple expression that can be minimized using numerical methods (E provides more details).

6 Harmonic trap

We now focus on a system diffusing in a one-dimensional harmonic trap defined by the potential

[TABLE]

Here $k_{\rm t}$ is the trap strength and the control parameter is the time-dependent trap minimum $\lambda_{i}$. The work required to perform an $N$-step discrete control protocol $\Lambda$, taking the control parameter from its initial value $\lambda_{0}$ to $\lambda_{N}$, can be calculated exactly.

6.1 Infinite-time limit

In the infinite-time limit (§3), the excess work for a single step $\lambda_{i}\to\lambda_{i+1}$ is the relative entropy (8) between the equilibrium distributions (1) at $\lambda_{i}$ and $\lambda_{i+1}$:

[TABLE]

Thus the infinite-time work for a discrete protocol of $N$ steps is

[TABLE]

Based on the convexity of this expression, equal step sizes $\Delta\lambda_{i,i+1}=\Delta\lambda_{\rm tot}/N$ minimize the infinite-time work,

[TABLE]

which scales with the number of steps as $1/N$ [32]. For this simple system (10mopry), the small-step approximation of the relative entropy (10l) is exact, for arbitrary step sizes.

6.2 General solution: finite-time work

Finite-duration control protocols feature both the infinite-time excess work and the time-dependent contribution (§4). For a system initially in equilibrium at the initial control parameter $\lambda_{0}$, the average excess work (in this case equal to the total work since $\Delta F=0$) for the $(i+1)$th step is

[TABLE]

where

[TABLE]

($\xi_{-1}=0$ as there are no terms in that summation.) D provides a detailed derivation of (10moprac,10moprada).

Summing (10moprac) over the entire protocol $\Lambda$ gives

[TABLE]

For this simple system, the normalized force autocovariance (the force autocorrelation) is $\langle\delta f(0)\delta f({\Delta t_{i}})\rangle_{\lambda_{i}}/\langle\delta f^{2}\rangle_{\lambda_{i}}=\exp(-\beta Dk_{\rm t}{\Delta t_{i}})$, so the approximate excess work within the linear-response regime (10mopq) is

[TABLE]

which is equivalent to truncating $\xi_{i-1}$ from (10mopradb) after the first term: $\xi_{i-1}\approx\Delta\lambda_{i-1,i}$. From (10moprada), it follows that the linear-response approximation (10mopq) holds when

[TABLE]

which is satisfied in the limit of long times ${\Delta t_{i}}\gg 1/(\beta Dk_{\rm t})$ spent at each control parameter.

For a protocol $\Lambda$ of protocol duration $\tau$ and consisting of $N$ control parameter steps, each of uniform size $\Delta\lambda_{i,i+1}=\Delta\lambda_{\rm tot}/N$, with uniform time allocations ${\Delta t_{i}}=\tau/(N-1)\equiv{\Delta t}_{\rm step}$, the exact excess work is

[TABLE]

the infinite-time work is (10mopraa), and for large $N$ the linear-response work is

[TABLE]

In each case, for a fixed duration ${\Delta t}_{\rm step}$ allocated to each control parameter value, the work scales asymptotically ($N\to\infty$) as $1/N$. Figure 1a and b show the average excess work for $N=10$ and $\Delta\lambda_{\rm tot}/N=1$, and the difference between the average work and the infinite-time limit as a function of the step duration. For sufficiently large step duration, the exact result (10mopradah) converges to the linear-response prediction (10mopradai) and exponentially approaches the infinite-time limit. The three curves have a fixed ordering: the exact solution (10mopradah) has a series of positive terms added beyond the linear-response expression (10mopradai), which in turn has an extra positive term added beyond the infinite-time limit (10moprab).

For protocol durations sufficiently long that the time spent at each control parameter significantly exceeds the relaxation time, the linear-response approximation and the exact result converge. Furthermore, neglecting the $\xi_{i-1}$ term reduces the exact protocol work (10mopradae) to the infinite-time limit (10mopraa). Figure 1c shows, for the particular step-duration ${\Delta t}_{\rm step}=(\beta Dk_{\rm t})^{-1}$, the $1/N$ scaling (for large $N$) of the average protocol work.

7 Periodic potential

Now we consider a single diffusing particle in a one-dimensional energy landscape $E(x,\lambda_{i})$ consisting of two components: a control-parameter-independent periodic potential, and a control-parameter-dependent harmonic trapping potential:

[TABLE]

where $k_{\rm t}$ is the harmonic trap strength, $E^{\ddagger}$ is the energetic barrier height between adjacent minima on the periodic potential, and $\ell$ is the period (Fig. 2). This potential represents a system with a sequence of metastable states, such as those often found in models of molecular machines [39].

Figure 2b shows numerical estimation of the autocovariance $\langle\delta f(0)\delta f({\Delta t})\rangle_{\lambda}$ from equilibrium simulations at several fixed control parameters evenly spaced over a single period of the underlying potential (10mopradaj). Using the force autocovariance as input to the linear-response approximation (10mopq), we minimize the average excess work during a discrete control protocol $\Lambda$ with a fixed number $N$ of steps. F gives details on the equilibrium simulations and numerical optimization of the excess work.

We consider three different protocol optimization schemes in order to isolate the effects of the optimal allocation of times ${\Delta t_{i}}$ to a fixed ‘naive’ sequence of control parameter values (a ‘time-optimized’ protocol), the optimal placement of control parameters for a fixed ‘naive’ set of time allocations (‘space-optimized’), and the simultaneous optimization of time allocations and control parameter placements (‘fully optimized’). In all cases, protocols are constrained by having fixed protocol duration $\tau$ (4), protocol endpoints $\lambda_{0},\lambda_{N}$, and number of steps $N$. In order to minimize the effect of the boundary conditions, we consider control protocols which traverse several periodic repetitions of the underlying potential.

For such a discrete protocol, Fig. 3 shows the time allocations ${\Delta t_{i}}$ and control parameter step-sizes $\Delta\lambda_{i,i+1}$, relative to their naive values ${\Delta t}_{\rm naive}\equiv\tau/(N-1)$ and $\Delta\lambda_{\rm naive}\equiv\Delta\lambda_{\rm tot}/N$, as a function of the control parameter value $\lambda$ over a single period, for a protocol with $N_{\rm p}=6$ steps per potential period and protocol duration $\tau=4(N-1)\overline{\tau_{\rm relax}}$, where $\overline{\tau_{\rm relax}}$ is the mean integral relaxation time over a single period of the underlying potential.

In each case, the behavior predicted by our theoretical analysis of simplified systems in §5 is borne out. In particular, time-optimized protocols allocate a larger fraction of the protocol duration to regions where the force is slowly relaxing (10moprx), while space-optimized protocols take step sizes which are largest in regions where the force variance is small and rapidly relaxing. In the fully optimized protocols, both behaviors are present.

Figure 4 shows the theoretically expected excess work for these minimum-work protocols, specifically the predicted excess work ratio $\langle W_{\rm ex}^{*}\rangle_{\Lambda}/\langle W_{\rm ex}^{\rm naive}\rangle_{\Lambda}$ for the three distinct protocol classes: time-optimized, space-optimized, and fully optimized. For short durations $(\widetilde{{\Delta t}}\equiv\tau/(N\overline{\tau_{\rm relax}})<1)$, time optimization yields no gain over naive protocols, while spatial optimization and full optimization are indistinguishable. For intermediate durations ($\widetilde{{\Delta t}}\approx 2$), time optimization has maximum effect, and full optimization significantly improves upon spatial optimization. For longer durations ($\widetilde{{\Delta t}}\gtrapprox 8$), time optimization again gives no benefit over the naive protocol, as the time-dependent term in (10mopq) becomes negligible.

Figure 5 shows that as the number $N_{\rm p}$ of steps per periodic image increases, the time allocation for fully optimized discrete control protocols converges to that of the optimal continuous protocol derived from the generalized friction coefficient [21].

However, there is a significant difference between the discrete at low step numbers $N_{\rm p}$ and continuous control protocols [21, 15]. In particular, relative to an optimal continuous protocol, fully optimized discrete protocols allocate a smaller fraction of their duration at $\lambda=\ell/2$, near the energy barrier (and in fact at the lowest step numbers completely avoid this region). This ability of low step-number protocols to entirely avoid regions of control parameter space with high force variance and slow relaxation (generally speaking, near energy barriers) represents a qualitatively distinct optimization strategy that is simply not available to continuous protocols.

Furthermore, the continuous protocols allocate time symmetrically about the energy barrier because the generalized friction maintains the same symmetries as the underlying energetic landscape (10mopradaj) [21]. As a result, a continuous optimal protocol traverses the same path in both the forward $\lambda_{0}\to\lambda_{\rm N}$ and reverse $\lambda_{\rm N}\to\lambda_{0}$ directions. Discrete protocols break this symmetry because of the infinite-time contribution (§3) as, in general, $D[p(x)||q(x)]\neq D[q(x)||p(x)]$. For small steps (10l), this asymmetry persists; in (10mopq) the excess work during the control parameter step $\Delta\lambda_{i,i+1}$ is a function of the force variance at the current control parameter $\langle\delta f^{2}\rangle_{\lambda_{i}}$ (and independent of the force variance $\langle\delta f^{2}\rangle_{\lambda_{i+1}}$ at the destination control-parameter value $\lambda_{i+1}$). This produces a directional asymmetry as the excess work for the control parameter step $\Delta\lambda_{i+1,i}$ is generally different than the excess work for step $\Delta\lambda_{i,i+1}$. However, as the number of control parameter steps increases and the distance between those steps becomes sufficiently small, the difference between the force variance at consecutive control parameter values becomes negligibly small, and the asymmetry between forward and reverse protocols vanishes (Fig. 5).

8 Discussion

In this article, we derived the work required to drive a microscopic system out of equilibrium via a discrete control protocol. Such a control protocol transforms the energy landscape through a series of discrete intermediate states, capturing the discrete nature of the chemical reaction sequences that drive many biological molecular motors. The central result is the linear-response expression for excess work (10mopq), which quantifies the near-equilibrium work of a particular control protocol, as a function solely of the equilibrium system properties.

We deduced a general expression for the work required to make a discrete change in the control parameter vector of a system in equilibrium (10b) and used this to exactly quantify the work required to perform a discrete protocol in the infinite-time limit (§3). When each step is sufficiently small and hence each perturbation is sufficiently weak, our derivation reduces to previously known results [32, 31]. Our primary contribution is to generalize these analyses outside of the infinite-time limit, where we use a linear-response approximation to derive the leading-order time-dependent contribution to the excess work (§4). Theoretically, this work goes significantly beyond previous efforts [30] to quantify energy flows in discretely driven nonequilibrium systems, by incorporating the effects of relaxation kinetics that are non-exponential and that vary across control parameter space, and by simultaneously optimizing both the placement of control parameter values as well as the allocation of times.

We investigated the correspondence between our linear-response approximation and an exact solution for a harmonically trapped Brownian particle driven by a series of discrete steps of equal size. We also studied the optimal allocation of time and placement of control parameter values that minimize the work for protocols traversing many repetitions of a periodic energy landscape (10mopradaj). We find that fully optimized discrete control protocols have qualitatively distinct features when compared to their continuous-protocol analogs. In particular, discrete protocols do not obey the same directional symmetry that continuous protocols do, and in the context of the periodic potential (§7), discrete protocols allocate a smaller fraction of their total duration near the energy barrier. More generally, minimum-work protocols allocate more time to regions where the force has a smaller variance and is more slowly decaying. Finally, we quantified the reduction in excess work, relative to a naive protocol, achieved by these minimum-work discrete control protocols. In particular, the theoretical excess work reduction (relative to a naive protocol) of a fully optimized protocol exceeds $50\%$ for small step numbers and short protocol durations. Significant reduction persists even for intermediate durations ($\tilde{\tau}\approx 2$) when fully optimized, highlighting the benefits of both optimized placement of control parameter values and the allocation of time among them.

The paradigm of a discretely driven nonequilibrium system is motivated, in part, due to its resemblance to the chemical driving in many biomolecular machines. Stochastic protocol ensembles have been recently considered in related work [40], and future investigation on how the stochastic properties of chemical driving affect the dissipation in discretely driven systems promises a more robust framework within which to compare the theoretical predictions for optimal operation to experimental results on the natural operation of molecular machines.

9 Acknowledgements

The authors thank John Bechhoefer and Emma Lathouwers (SFU Physics) and Miranda Louwerse (SFU Chemistry) for insightful comments on the article. This work is supported by Natural Sciences and Engineering Research Council of Canada (NSERC) CGS Masters and Doctoral Fellowships (S.J.L.), an NSERC Discovery Grant (D.A.S.), a Tier-II Canada Research Chair (D.A.S.), and WestGrid (www.westgrid.ca) and Compute Canada Calcul Canada (www.computecanada.ca).