Optimal discrete control: minimizing dissipation in discretely driven nonequilibrium systems
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.
December 29, 2018
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 . For instance, the rotary ATP-synthase motor produces the cellular energy currency ATP , 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 ATPase motor, for instance, can operate at speeds up to 350 rotations per second , 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  and the reversal of magnetization in an Ising magnet , 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 . While this applies well to the single-molecule experimental paradigm, the biomolecular machines that we hope to understand 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 discrete perturbations to a thermodynamic system, as opposed to a continuous driving process.
In this article, we develop a theoretical framework for nonequilibrium control using discrete control parameter changes. Our central result (10mopq) quantifies the nonequilibrium energetic costs associated with discretely driving a microscopic system. 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 5. This work complements and generalizes previous results on the entropy production associated with discrete processes . In the continuum limit, our formalism reduces to previously known results, namely the thermodynamic-length formalism introduced by Crooks  and the entropy-differential metric of Burbea and Rao , and is related to the generalized friction coefficient of Sivak and Crooks . Finally, we demonstrate in two model systems the quantitative implications of this theory for discretely driven systems.
We consider a system in contact with a heat bath, subject to a set of experimentally controlled parameters , such that at equilibrium the distribution over microstates is
where is the inverse temperature, is the system Hamiltonian, and is the equilibrium free energy at control parameter vector .
A control protocol is a particular time-dependent perturbation applied to the control parameter vector to transform it between an initial and final in a prescribed time . For a given control protocol, the system responds stochastically. Across the entire control protocol , the average amount of excess work (supplied by an external source)—or work required above and beyond the equilibrium free energy difference between the initial and final control parameter values—is , where indicates an average of system responses to the control protocol .
Here, we consider discrete control protocols which consist of a series of instantaneous perturbations for consecutive control parameter values and . The system spends a prescribed time at each control parameter . Thus, each protocol is defined by a set of control parameter values and the associated times spent at them:
Previous work has considered energetic flows in discrete-stepping processes [22, 25, 26, 27], 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 ) associated with a particular discrete perturbation that transforms the control parameter from to is
where is the (generally nonequilibrium) distribution over system microstates at the time that the control parameter changes from to .
We consider protocols which start in equilibrium at initial control parameter value , equivalent to taking the time spent at to infinity, . A particular control protocol begins with the first control parameter change , and finishes when the control parameter arrives at its terminal value . The protocol duration 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 control parameter values , the protocol duration is
and the total control parameter displacement for fixed control parameter endpoints is
For given control parameter endpoints and (hence given protocol displacement ) and duration , minimizing work involves choosing intermediate control parameter values and associated dwell times .
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 (1). The average work required to discretely change the control parameter vector from to is (3), with the equilibrium initial distribution (1):
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, . The average work from (6) can then be written solely in terms of equilibrium distributions,
where is the difference between the equilibrium free energy at control parameter and . The integral in (7) is the relative entropy (or Kullback-Leibler divergence)  between the equilibrium distributions before () and after () the control parameter change,
A protocol consists of such control parameter steps , for . If at each control parameter value the system fully equilibrates, then the average work associated with any step 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,
where is the equilibrium free energy change between the initial and final control parameter values. Thus the average excess work is
For sufficiently small control parameter steps , the relative entropy in (10b) can be Taylor expanded about its current value to yield
where is the equilibrium covariance of conjugate forces at control parameter (see A for details) . Throughout, we employ the Einstein summation notation, where repeated indices are implicitly summed over.
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 . For a system at equilibrium for control parameter , the energy at the next control parameter value in the protocol can be linearly approximated as
where is the vector of conjugate forces with elements . When the control parameter instantaneously changes from to , the time-dependent relaxation of towards its equilibrium value at is, under the linear-response approximation,
where indicates an average over the instantaneous nonequilibrium system distribution after relaxing (under the Hamiltonian ) for a time starting from the equilibrium distribution at . Both the force autocovariance and the average force on the RHS are taken over the equilibrium ensemble at fixed control parameter value . B provides a detailed derivation of (10mn).
If the system relaxes for a time at control parameter before the next control parameter change , this step requires average work (3)
for infinite-time work from §3 and linear-response correction due to incomplete system relaxation.
(10moc) is only strictly valid if the system was at equilibrium for before the step to , so more generally the work required for the control parameter change 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. Within this limit, the total average excess work required to perform a discrete control protocol is
where (10mopc) uses the small-step limit of the infinite-time excess work (10l). The assumption that the system begins in equilibrium at ensures that the term does not contribute to the total nonequilibrium excess work. By setting , both sums in (10mop) can be taken over the same index range, leading to a more compact form for the excess work:
This captures the combined effects of the control parameter step sizes and time allocations on the excess work during a discrete control protocol . 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 step-size factor).
Equation (10mopq) generalizes the near-equilibrium expression for the dissipation of a discrete control protocol in , 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 as well as the allocation of times . 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, the linear-response correction to the excess work recovers the generalized friction formalism from .
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 . 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
Interpretation is most immediate in the continuous-protocol limit, where each is the distance along an infinitesimal segment of the control protocol , measured with respect to the metric ; therefore, the sum over all steps gives the thermodynamic length between the initial and final equilibrium system macrostates [23, 22].
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:
The lower bound is saturated if and only if the are identical,
Along an optimal protocol (indicated by the superscript ), the condition (10moprt) implies that , and thus equal excess work is done during each step of the protocol.
For a single control parameter with fixed endpoints and a given set of time allocations , the condition (10moprt) implies the optimal placement of control parameter values through the proportionality . Thus optimal control parameter placement tends to avoid regions with large force variance and slowly decaying force autocovariance, subject to the quadratic cost 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 , our predictions reduce to previous calculations by Nulton et al.  of the optimal placement of discrete steps. In particular, for a single control parameter, the condition (10moprt) implies that optimal protocols have the proportionality . Furthermore, in the continuous-protocol limit, the infinite-time thermodynamic length between the initial and final control parameters (measured with respect to ) converges to that of Crooks .
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:
Using Lagrange multipliers, the optimal allocation of time among a fixed set of control parameter values, subject to the protocol duration constraint (4), is
where (see E for a detailed derivation). This result was also found by Nulton et al. in . In the long-duration limit, where , the second RHS term in (10moprv) is negligible, and the optimal allocation of time takes on the simple form
Intuitively, this implies that along minimum-work protocols, more time is allocated to regions where the integral relaxation time  is larger.
For more general scenarios, even with one control parameter, analytic optimization methods become cumbersome, and numerical methods are most practical (E provides more details).
6 Harmonic trap
We now focus on a system diffusing in a one-dimensional harmonic trap defined by the potential
Here is the trap strength and the control parameter is the time-dependent trap minimum . The work required to perform an -step discrete control protocol , taking the control parameter from its initial value to , can be calculated exactly.
6.1 Infinite-time limit
Thus the infinite-time work for a discrete protocol of steps is
Based on the convexity of this expression, equal step sizes minimize the infinite-time work,
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 , the average excess work (in this case equal to the total work since ) for the th step is
Summing (10moprab) over the entire protocol gives
For this simple system, the normalized force autocovariance (the force autocorrelation) is , so the approximate excess work within the linear-response regime (10mopq) is
which is satisfied in the limit of long times spent at each control parameter.
For a protocol of protocol duration and consisting of control parameter steps, each of uniform size , with uniform time allocations , the exact excess work is
the infinite-time work is (10moprz), and the linear-response work is
In each case, for a fixed duration allocated to each control parameter value, the work scales as . Figure 1a and b show the average excess work for and , and the difference between the average work and the infinite-time limit as a function of the step duration. For sufficiently large protocol step duration, the exact result (10mopracag) converges to the linear-response prediction (10mopracah) and exponentially approaches the infinite-time limit.
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 term reduces the exact protocol work (10mopracad) to the infinite-time limit (10moprz). Figure 1c shows, for the particular step-duration , the scaling (for large ) of the average protocol work.
7 Periodic potential
Now we consider a single diffusing particle in a one-dimensional energy landscape consisting of two components: a control-parameter-independent periodic potential, and a control-parameter-dependent harmonic trapping potential:
where is the harmonic trap strength, is the energetic barrier height between adjacent minima on the periodic potential, and 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 .
Figure 2b shows numerical estimation of the autocovariance from equilibrium simulations at several fixed control parameters evenly spaced over a single period of the underlying potential (10mopracai). Using the force autocovariance as input to the linear-response approximation (10mopq), we minimize the average excess work during a discrete control protocol with a fixed number 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 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 (4), protocol endpoints , and number of steps . 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 and control parameter step-sizes , relative to their naive values and , as a function of the control parameter value over a single period, for a protocol with steps per potential period.
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 (10moprw), 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 for the three distinct protocol classes: time-optimized, space-optimized, and fully optimized. In the limit of , time optimization yields no gain over naive protocols, while spatial optimization and full optimization are indistinguishable. For intermediate durations (), time optimization has maximum effect, and full optimization significantly improves upon spatial optimization. In the long-duration limit, 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 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 .
However, there is a significant difference between the discrete and continuous control protocols at low step numbers [21, 15]. In particular, relative to an optimal continuous protocol, fully optimized discrete protocols allocate a smaller fraction of their duration at , near the energy barrier.
Furthermore, the continuous protocols allocate time symmetrically about the energy barrier because the generalized friction maintains the same symmetries as the underlying energetic landscape (10mopracai) . As a result, a continuous optimal protocol traverses the same path in both the forward and reverse directions. Discrete protocols break this symmetry because of the infinite-time contribution (§3) as, in general, . In the small-step limit this asymmetry persists; in (10mopq) the excess work during the control parameter step is a function of the force variance at the current control parameter (and independent of the force variance at the destination control-parameter value ). This produces a directional asymmetry as the excess work for the control parameter step is generally different than the excess work for step .
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, reminiscent 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). In the small-step limit, where each step is sufficiently small and hence each perturbation is sufficiently weak, our derivation reduces to previously known results [24, 23]. 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  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 (10mopracai). 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 for small step numbers and short protocol durations. Significant reduction persists even for intermediate durations () 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 , and future investigation on the synthesis of these ideas 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.
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).
-  Howard J 2001 Mechanics of Motor Proteins and the Cytoskeleton (Sinauer Associates)
-  Yoshida M, Muneyuki E and Hisabori T 2001 Nat. Rev. Mol. Cell Biol. 2 669–677
-  Hirokawa N and Noda Y 2009 Nat. Rev. Mol. Cell Bio. 10 682–696
-  Sellers J R 2000 BBA-Mol. Cell Res. 1496(1) 3–22
-  Ueno H, Suzuki T, Kinosata K and Yoshida M 2004 Proc. Natl. Acad. Sci. 102 1333–1338
-  Liu S, Chistol G, Hetherington C L, Tafoya S, Aathavan K, Schnitzbauer J, Grimes S, Jardine P J and Bustamante C 2014 Cell 157 702–713
-  Yasuda R, Noji H, Kinosata K J and Yoshida M 2001 Cell 93 1117–1124
-  Berndsen Z T, Keller N, Grimes S, Jardine P J and Smith D E 2014 Proc. Natl. Acad. Sci. U.S.A. 111 8345–8350
-  Bialek W 2012 Biophysics: Searching for Principles (Princeton University Press)
-  Niven J E and Laughlin S B 2008 J. Exp. Biol. 211 1792–1804
-  Chen Y J, Groves B, Muscat R A and Seelig G 2015 Nat. Nanotechnol. 10 748–760
-  Peng H, Li X F, Zhang H and Le C 2017 Nat. Comm. 8 14378
-  Zulkowski P R and DeWeese M R 2014 Phys. Rev. E. 89
-  Rotskoff G M and Crooks G E 2015 Phys. Rev. E 92 060102(R)
-  Schmiedl T and Seifert U 2007 Phys. Rev. Lett. 98
-  Gomez-Marin A, Schmiedl T and Seifert U 2008 J. Chem. Phys. 129 024114
-  Esposito M, Kawai R, Lindenberg K and Van den Broeck C 2010 Phys. Rev. Lett. 105 150603
-  Zulkowski P R, Sivak D A, Crooks G E and DeWeese M R 2012 Phys. Rev. E. 86
-  Zulkowski P R and DeWeese M R 2015 Phys. Rev. E. 92
-  Sivak D A and Crooks G E 2016 Phys. Rev. E. 94
-  Sivak D A and Crooks G E 2012 Phys. Rev. Lett. 108 190602
-  Nulton J, Salamon P, Andresen B and Amin Q 1985 J. Chem. Phys. 83 334–338
-  Crooks G E 2007 Phys. Rev. Lett. 99 100602
-  Burbea J and Rao C R 1982 Journal of Multivariate Analysis 12 575–596
-  Andreson B, Salamon P and Berry R S 1984 Phys. Today 37 62–70
-  Salamon P and Berry R S 1983 Phys. Rev. Lett. 51 1127–1130
-  Salamon P, Nulton J and Ihring E 1984 J. Chem. Phys. 80 436–437
-  Cover T M and Thomas J A 2006 Elements of Information Theory 2nd ed (Wiley)
-  Chandler D 1987 Introduction to Modern Statistical Mechanics (Oxford Univeristy Press)
-  Garanin D A 1996 Phys. Rev. E 54 3250–3256
-  Reimann P 2002 Phys, Rep. 361 57–265
-  Large S J, Chetrite R and Sivak D A 2018 Europhys. Lett. 124 20001
-  Gardiner C 2009 Stochastic Methods: A Handbook for the Natural and Social Sciences 4th ed (Springer-Verlag)
-  Winston W L 1978 Operations Research: Applications and Algorithms (Duxbury Press)
Appendix A Expansion of the relative entropy
The relative entropy (Kullback-Leibler divergence) between two continuous probability distributions and is defined as 
In the context of the present work, the equilibrium distribution is parameterized by the control variable . The integrand of the relative entropy for two consecutive equilibrium distributions at and is
For small changes in the control parameter, we Taylor expand Eq. (10mopracak) about ,
where is the partial derivative of with respect to the th component of the control parameter , indicates that the argument is evaluated at , and we have made use of the Einstein summation notation, where repeated indices are summed over.
The first term in (A) is
The derivative on the RHS of (A), evaluated at is
and the second-derivative term in (A) is
Equation (10mopracan) can be simplified by noting that the equilibrium probability distribution is normalized, , and partial differentiation commutes with integration, so substituting (10mopracan) into the relative entropy expression (10mopracaj), gives
so this term does not contribute to the overall relative entropy. This results from the relative entropy being a convex function with a minimum at . In analogy with (10mopracap), the second term on the RHS of (10mopracao) also vanishes,
where the integral is the Fisher information matrix at control parameter . In the small-step limit, the term is negligible, so for a discrete control protocol