Shortcuts in stochastic systems and control of biophysical processes

The biochemical reaction networks that regulate living systems are all stochastic to varying degrees. The resulting randomness affects biological outcomes at multiple scales, from the functional states of single proteins in a cell to the evolutionary trajectory of whole populations. Controlling how the distribution of these outcomes changes over time—via external interventions like time-varying concentrations of chemical species—is a complex challenge. In this work, we show how counterdiabatic (CD) driving, first developed to control quantum systems, provides a versatile tool for steering biological processes. We develop a practical graph-theoretic framework for CD driving in discrete-state continuous-time Markov networks. Though CD driving is limited to target trajectories that are instantaneous stationary states, we show how to generalize the approach to allow for non-stationary targets and local control—where only a subset of system states are targeted. The latter is particularly useful for biological implementations where there may be only a small number of available external control knobs, insufficient for global control. We derive simple graphical criteria for when local versus global control is possible. Finally, we illustrate the formalism with global control of a genetic regulatory switch and local control in chaperone-assisted protein folding. The derived control protocols in the chaperone system closely resemble natural control strategies seen in experimental measurements of heat shock response in yeast and E. coli.

A fundamental dichotomy for biological processes is that they are both intrinsically stochastic and tightly controlled.The stochasticity arises from the random nature of the underlying biochemical reactions, and has significant consequences in a variety of contexts: gene expression [1], motor proteins [2], protein folding [3], all the way up to the ecological interactions and evolution of entire populations of organisms [4,5].Theories for such systems often employ discrete state Markov models (or continuum analogues like Fokker-Planck equations) which describe how the probability distribution of system states evolves over time.On the other hand, biology utilizes a wide array of control knobs to regulate such distributions, most often through time-dependent changes in the concentration of chemical species that influence state transition rates.In many cases these changes occur due to environmental cues-either threatening or beneficial-and the system response must be sufficiently fast to avoid danger or gain advantage.
The interplay of randomness and regulation naturally leads us to ask about the limits of control: to what extent can a biological system be driven through a prescribed trajectory of probability distributions over a finite time interval?Beyond curiosity over whether nature actu-ally tests these limits in vivo, this question also arises in other contexts.In synthetic biology [6] one may want to precisely specify the probabilistic behavior of genetic switches or other regulatory circuit components in response to a stimulus.Control of a system is generally easiest to describe and quantify if the perturbation is applied slowly (adiabatically).The advantage of this assumption is that, at each moment of the control protocol, the approximate form of the state probability distribution is known from equilibrium thermodynamics.However in natural settings, responses to rapid environmental changes may entail sharp changes in the concentrations of biochemical components.For instance, an ambient temperature increase of even a few degrees can significantly increase the probability that proteins misfold and aggregate.In response to such "heat shock", cells quickly upregulate the number of chaperones-specialized proteins that facilitate unfolding or disaggregating misfolded proteins [7][8][9][10][11][12].
There is no guarantee that the quasi-equilibrium assumption holds throughout such a process, and thus the standard tools of equilibrium or near-equilibrium thermodynamics (i.e.linear response theory) do not necessarily apply.
If the system is driven over a finitetime interval, subject to fluctuations that take it far from equilibrium, can there still be a degree of control?We can pose this question more concretely, as illustrated in Fig. 1A.A biological system is typically part of a larger FIG. 1.
Schematic of the biological control problem: A) A system of interest (here a membrane receptor protein) is described via a network of transition rates between discrete states.Certain rates may be influenced by factors external to the system, which we denote as control parameters.For biochemical systems these are often concentrations of chemical species (ligands, ATP, etc.) or environmental factors like temperature.B) We consider two types of control: global control, where we require the probability of every state to follow a target trajectory over a finite time interval; and local control, where we impose this requirement on only a subset of states.C) In either case, the goal is to find whether control is possible for a given target, and if so calculate the control parameter protocol that forces the system to follow the target trajectory.
complex of interacting components.If we focus on a system of interest, and describe it via a discrete state Markov model, the transition rates between states may depend on factors external to the system, like concentrations of ligands that bind to the system, or energy molecules like ATP that are required to fuel certain reactions.In certain experimental or synthetic biology contexts these factors may be under direct human control, but in natural contexts they are often the product of autonomous processes outside the system of interest, like the temperature fluctuations that lead to heat shock.In either case we will denote these external factors for simplicity as control parameters, and investigate their influence on the system dynamics.We will consider a specific question of controllability: can one find a control protocol (a time-dependent function of the parameters) such that the probability distribution of system states follows a certain sequence of target distributions over a finite time interval?We can define two types of control (Fig. 1B): global control, where we demand the probability of every state in the system follow a chosen trajectory; and local control, where we only care about a subset of states following the target, and allow the remainder to have arbitrary dynamics.Ideally, we would like criteria for what kinds of target trajectories are achievable in a given system, and a procedure to calculate the protocol if the target is possible (Fig. 1C).Answering these questions would not only give us new tools to precisely manipulate biological systems in experiments, but shed light on dynamical constraints in vivo.For example, by exploring how controllability depends on the duration of the target trajectory, one can investigate limits on how quickly a system can alter its state distribution in response to an external environmental change.
Interestingly, for one particular class of target trajectories-forcing the system to mimic quasiequilibrium behavior-the situation strongly resembles questions from quantum control and quantum thermodynamics [13], where a new line of research has been dubbed "shortcuts to adiabaticity".In recent years a great deal of theoretical and experimental work has been dedicated to mathematical tools and practical schemes to suppress nonequilibrium excitations in finite-time, nonequilibrium processes.To this end, a variety of techniques have been developed: the use of dynamical invariants [14], the inversion of scaling laws [15], the fastforward technique [16][17][18][19][20][21][22][23], optimal protocols from optimal control theory [24][25][26][27], optimal driving from properties of quantum work statistics [28], "environment" assisted methods [29,30], using the properties of Lie algebras [31], and approximate methods such as linear response theory [32][33][34][35], fast quasistatic dynamics [36], or time-rescaling [37,38], to name just a few.See Refs.[39,40] and references therein for comprehensive reviews of these techniques.
Among this plethora of different approaches, counterdiabatic (CD) or transitionless quantum driving stands out, since it is the only method that suppresses excitations away from the adiabatic manifold at all instants.In this paradigm [41][42][43][44] one considers a timedependent Hamiltonian H 0 (t) with instantaneous eigenvalues { n (t)} and eigenstates {|n(t) }.In the adiabatic limit no transitions between eigenstates occur [45], and each eigenstate acquires only a time-dependent phase that can be separated into a dynamical and a geometric contribution [46].In other words, if we start in a particular eigenstate |n(0) at t = 0, we remain in the corresponding instantaneous eigenstate |n(t) at all later times, up to a phase.The goal of CD driving is to make the system follow the same target trajectory of eigenstates as in the adiabatic case, but over a finite time.
To accomplish this, a CD Hamiltonian H(t) can be constructed, such that the adiabatic approximation associated with H 0 (t) is an exact solution of the dynamics generated by H(t) under the time-dependent Schrödinger equation.It is reasonably easy to derive that time-evolution under [41][42][43], (1) maintains the system on the adiabatic manifold.Note that it is the auxiliary Hamiltonian H 1 (t) that enforces evolution along the adiabatic manifold of H 0 (t): if a system is prepared in an eigenstate |n(0) of H 0 (0) and subsequently evolves under H(t), then the term H 1 (t) effectively suppresses the non-adiabatic transitions out of |n(t) that would arise in the absence of this term.
To date, a few dozen experiments have implemented and utilized such shortcuts to adiabaticity to, for instance, transport ions or load BECs into an optical trap without creating parasitic excitations [40].However, due to the mathematical complexity of the auxiliary Hamiltonian (1), counterdiabatic driving has been restricted to "simple" quantum systems.Note that in order to compute H 1 (t) one requires the instantaneous eigenstates of the unperturbed Hamiltonian, which is practically, conceptually, and numerically a rather involved task.
On the other hand, the scope of CD driving is not limited to the quantum realm.Because of the close mathematical analogies between classical stochastic systems and quantum mechanics, it was recently recognized that the CD paradigm can also be formalized for classical scenarios [23,44,[47][48][49][50][51][52][53].The classical analogue of driving a system along a target trajectory of eigenstates is a trajectory of instantaneous stationary distributions.Last year, our group and collaborators developed the first biological application of CD driving: controlling the distribution of genotypes in an evolving cellular population via external drug protocols [54].This type of "evolutionary steering" has various potential applications, most notably in designing strategies to combat drug resistance in bacterial diseases and tumors.The CD formalism in this case was built around a multi-dimensional Fokker-Planck model, generalizing the one-dimensional Fokker-Planck approach of Ref. [50].
Our current work generalizes these initial results in two significant ways: i) We provide a universal framework capable of handling the wide diversity of stochastic models used in biology, taking advantage of graph theory to construct general algorithms that can be applied to discrete state systems of arbitrary complexity.The discrete state formalism presented here includes the continuum Fokker-Planck theory as a special case.ii) Our earlier results were limited to target trajectories that were instantaneous equilibrium distributions (CD driving) defined for all states (global control).Here we relax both those assumptions: we allow arbitrary target distributions, and the possibility for defining targets on only a subset of states (local control).Thus our new formalism includes for example the case of fast-forward driving [53], where the target trajectory begins and ends in equilibrium, but can be arbitrary in between.The usefulness of our method is of course not confined to biology, but is relevant to other classical systems described by Markovian transitions between states.However biology provides a singularly fascinating context in which to explore driving, both because it sheds light on the possibility of control in complex stochastic systems with many interacting components, and provides an accessible platform for future experimental tests of these ideas.
Outline: In Sec.I we start with the most basic version of the theory, formulating CD driving for any discrete state Markov model.By looking at the properties of the probability current graph associated with the master equation of the model, we can express CD solutions in terms of spanning trees and fundamental cycles of the graph.Beyond its practical utility, the graphical approach highlights the degeneracy of CD driving: the potential existence of many distinct, physically realizable CD protocols that drive a system through the same target trajectory of probability distributions.The graphical approach is schematically summarized in Fig. 2, highlighting the components in the most general form for CD solutions, Eq. (31).
In Sec.II we show how the formalism can be generalized to arbitrary (non-CD) target trajectories and local control.This discussion allows us to derive simple graphical criteria for when global versus local control is possible.The criteria can help us determine what types of target trajectories are achievable in individual biological systems, based solely on the structure of the corresponding Markovian networks.
In Sec.III we apply our formalism to two biological examples, illustrating global and local control respectively.The first is a repressor-corepressor genetic regulatory switch, and the second a chaperone protein that catalyzes the unfolding of a misfolded protein in response to a heat shock.The switch provides perhaps the simplest example where the parameters have been experimentally characterized and various driving solutions can be directly tested in vitro.For the chaperone system, we highlight the qualitative similarities between local control protocols for rapidly suppressing misfolded proteins and experimental measurements of heat shock response in yeast and E. coli.
Sec. IV concludes with connections to other areas of nonequilibrium thermodynamics and questions for future work.

I. GENERAL THEORY OF COUNTERDIABATIC DRIVING IN DISCRETE STATE MARKOV MODELS
A. Setting up the counterdiabatic driving problem

Master equation and the CD transition matrix
Consider an N -state Markov system described by a vector p(t) whose component p i (t), i = 1, . . ., N , is the probability of being in state i at time t.The distribution p(t) evolves under the master equation [55,56], ( The off-diagonal element Ω ij (λ t ), i = j, of the N × N matrix Ω(λ t ) represents the conditional probability per unit time to transition to state i, given that the system is currently in state j.The diagonal elements Ω ii (λ t ) = − j =i Ω ji (λ t ) ensure each column of the matrix sums to zero [55].The transition rates Ω ij (λ t ) depend on a control protocol: a set of time-varying external parameters, denoted collectively by λ(t) ≡ λ t .Ω(t) plays the role of the Hamiltonian H 0 (t) in the classical analogy.
The instantaneous stationary probability ρ(λ t ) associated with Ω(λ t ) is the right eigenvector with eigenvalue zero, Ω(λ t )ρ(λ t ) = 0. ( When λ t has a non-constant time dependence, ρ(λ t ) in general is not a solution to Eq. ( 2), except in the adiabatic limit when the control parameters are varied infinitesimally slowly, ∂ t λ t → 0. The sequence of distributions ρ(λ t ) as a function of λ t defines a target trajectory for the system, analogous to the eigenstate trajectory |n(t) in the quantum version of CD.Given an instantaneous probability trajectory ρ(λ t ) defined by Eq. ( 3), we would like to find a counterdiabatic (CD) transition matrix Ω(λ t , λt ) such that the new master equation, evolves in time with state probabilities described by ρ(λ t ).Here λt ≡ ∂ t λ t .We are thus forcing the system to mimic adiabatic time evolution, even when λt is nonzero.As we will see below, Ω(λ t , λt ) will in general depend both on the instantaneous values of the control parameters λ t and their rate of change λt .In the limit of adiabatic driving we should recover the original transition matrix, Ω(λ t , λt → 0) = Ω(λ t ).Solving for the CD protocol corresponds to determining the elements of the Ω(λ t , λt ) matrix in Eq. ( 4) given a certain ρ(λ t ).This corresponds to finding the CD Hamiltonian H(t) of Eq. (1) in the quantum case.
We can look at the counterdiabatic problem as a special case of a more general question: given a certain time-dependent probability distribution that is our target, what is the transition matrix of the master equation for which this distribution is a solution?In effect, this is the inverse of the typical approach for the master equation, where we know the transition matrix and solve for the distribution.

Representing the system via an oriented current graph
To facilitate finding CD solutions, we start by expressing the original master equation of Eq. ( 2) equivalently in terms of probability currents between states, where the current from state j to i is given by: We can interpret any pair of states (i, j) where either Ω ij (λ t ) = 0 or Ω ji (λ t ) = 0 at some point during the protocol as being connected via an edge on a graph whose vertices are the states i = 1, . . ., N .
Let E be the number of edges in the resulting graph.Define a numbering α = 1, . . ., E and an arbitrary orientation for the edges such that each α corresponds to a specific edge and choice of current direction.For example, if edge α was between states (i, j), and the choice of direction was from j to i, then we can define current J α (t) ≡ J ij (t) for that edge.Alternatively if the choice of direction was from i to j, then J α (t) ≡ J ji (t) = −J ij (t).We denote rates Ω ij (λ t ) oriented parallel to the edge direction as forward rates, and those oriented opposite as backward rates.In this way we associate the master equation with a directed graph, a simple example of which is illustrated in Fig. 3. Eq. ( 5) can be rewritten in terms of the oriented currents J α (t) as where J (t) is an E-dimensional vector with components J α (t), and ∇ is an N × E dimensional matrix known as the incidence matrix of the directed graph [57] (closely related to the stoichiometric matrix of Ref. [58]).∇ is given by where the components of the two matrices ∇ ± are defined as: The αth column of ∇ contains a single 1 and a single −1, since each edge must have an origin and a destination state.With these definitions, Eq. ( 6) can be recast as a relation between the vectors J (t) and p(t), where the E × N dimensional matrix G(λ t ) is given by Here diag(k + (λ t )) is an E × E dimensional diagonal matrix where the diagonal is k + (λ t ), the vector of forward rates associated with each edge.For example k + α (λ t ) = Ω ij (λ t ) if the αth edge is oriented from j to i. Similarly k − (t) is the vector of backward rates.Comparing Eq. (2) to Eqs. ( 7)- (11), we see that the matrix Ω(λ t ) = ∇G(λ t ).In the special case where the matrix G(λ t ) has a right singular vector with singular value zero, we say that the rates in the system satisfy instantaneous detailed balance.We will refer to k ± (λ t ) as the "original" rate protocol for the system, since they determine the original transition matrix Ω(λ t ) and hence the target ρ(λ t ) via Eq.(3).Throughout the text we will use "original" to consistently describe quantities associated with Ω(λ t ).On the other hand quantities associated with the CD matrix Ω(λ t , λt ), like the CD rates k ± (λ t ) described below, will always have tildes to distinguish them from the original case.
Conservation of probability is enforced by summing over rows in Eq. ( 7), since i ∇ iα = 0, and so N i=1 ∂ t p i (t) = 0. Since any given row of Eq. ( 7) is thus linearly dependent on the other rows, it is convenient to work in the reduced representation of the equation where we leave out the row corresponding to a certain reference state (taken to be state N ), Here p(t) = (p 1 (t), . . ., p N −1 (t)) and the (N − 1) × E dimensional reduced incidence matrix ∇ is equal to ∇ with the N th row removed.Our focus will be on systems where there is a unique instantaneous stationary probability vector ρ(λ t ) at every t.In this case the master equation necessarily corresponds to a connected graph in the oriented current picture [55].By a well known result in graph theory, both the full and reduced incident matrices ∇ and ∇ of a connected, directed graph with N vertices have rank N − 1 [57].This means that all N − 1 rows of ∇ are linearly independent for the systems we consider.
Having described the original master equation of Eq. (2) in terms of oriented currents, we can do the same for Eqs.(3) and (4).Let us define the oriented stationary current J α (t) for the distribution ρ(λ t ) as follows: if the αth edge is oriented from j to i then In vector form, analogous to Eq. ( 10), the current is given by The reduced representation of Eq. (3) corresponds to If the system rates satisfy instantaneous detailed balance, J (t) = 0, since ρ(t) is the right singular vector of G(λ t ) with singular value zero.However our CD approach works for the more general case where J (t) can be nonzero but Eq. ( 15) is satisfied.In fact, as we will see in Sec.II.A, we can also generalize our theory to completely arbitrary (non-CD) target trajectories ρ(t) where Eq. ( 15) no longer holds.For the CD master equation, Eq. ( 4), we define the oriented current The time dependence of J α is explicitly through λ t and λt , but we write it in more compact form as J α (t) to avoid cumbersome notation.The analogue of Eq. ( 10) is where G(t) has the same structure as Eq. ( 11) but with forward/backward rate vectors k ± (t) corresponding to the CD rates Ω ij (λ t , λt ).Finally, Eq. ( 4) can be expressed as 3. Counterdiabatic current equation Subtracting Eq. ( 15) from Eq. ( 18) we find where δJ (t) ≡ J (t) − J (λ t ) is the difference between the CD and stationary current vectors.For the CD problem, we are given the original matrix elements Ω ij (λ t ) and thus also have the corresponding stationary distribution values ρ i (λ t ) and stationary currents J α (λ t ).What we need to determine, via Eq.( 19), are the CD currents J (t).The following Secs.I.B through I.D detail the procedure for finding these currents.Once we know J (t), Sec.I E shows how to use Eq.(17) to solve for k ± (t), or equivalently the CD matrix transition rates Ω ij (λ t , λt ).By construction, these satisfy Eq. ( 4), and hence define a CD protocol for the system.As a first step, let us consider the invertibility of Eq. (19) to solve for δJ (t).The (N − 1) × E dimensional matrix ∇ is generally non-square: N (N − 1)/2 ≥ E ≥ N − 1 for a connected graph.Only in the special case of tree-like graphs (no loops) do we have E = N − 1 and a square (N − 1) × (N − 1) matrix ∇.Since the rank of ∇ is N − 1, as mentioned above, for tree-like graphs ∇ is invertible and Eq. ( 19) can be solved without any additional complications: As described in the next section, the elements of ∇ −1 for a tree-like graph can be obtained directly through a graphical procedure, without the need to do any explicit matrix inversion.
In the case where E > N − 1, the solution procedure is more involved, but the end result has a relatively straightforward form: the most general solution δJ (t) can always be expressed as a finite linear combination of a basis of CD solutions.How to obtain this basis, and its close relationship to the spanning trees and fundamental cycles of the graph, is the topic we turn to next.2. Overview of the graphical approach for deriving CD solutions.We start with a Markov model defined by a transition matrix Ω(λt) dependent on the control protocol λt.Associated with this is a graph with N states, E edges, and a target trajectory ρ(λt) consisting of instantaneous stationary states of Ω(λt).The eventual goal is to find the CD transition matrix Ω(λt, λt) where ρ(λt) is the solution to the associated master equation, Eq. (4).To facilitate this, we must first find the CD currents J (t), the main goal of the graphical approach.The most general form of the solution for J (t) is given by Eq. ( 31), and consists of two components: (i) a spanning tree CD solution δJ (1) (t), given by Eq. ( 23) and derived via the procedure outlined in Sec.I B; (ii) a linear combination of the fundamental basis cycle vectors c (γ) , γ = 1, . . ., ∆, where ∆ = E − N + 1, as described in Sec.I D. The coefficient functions Φγ(t) are arbitrary.Once the CD currents J (t) are known, we can use Eq.(34) to solve for the CD transition rates k ± (t) that determine Ω(λt, λt).

B. General graphical solution for the counterdiabatic protocol
The graphical procedure described in this and the following two sections, culminating in the general solution of Eq. (31), is summarized in Fig. 2. To illustrate the procedure concretely, we will use the two-loop system shown in Fig. 3A as an example, where N = 4, E = 5.The solution for this case is relevant to the biophysical model for chaperone-assisted protein folding discussed i n Sec.III.B. Fig. 3A shows the rates k ± α (λ t ) that determine the transition matrix Ω(λ t ), and Fig. 3B labels the oriented stationary currents J α (t), α = 1, . . ., E.
Every connected graph has a set of spanning trees: subgraphs formed by removing ∆ ≡ E − N + 1 edges such that the remaining N − 1 edges form a tree linking together all the N vertices.The number T of such spanning trees is related to the reduced incidence matrix through Kirchhoff's matrix tree theorem [57], T = det ∇ ∇ T .For the current graph of Fig. 3B, this matrix is and the number of trees is thus T = 8.Let us select one spanning tree to label as the reference tree.The choice is arbitrary, since any spanning tree can be a valid starting point for constructing the basis.The left side of Fig. 3C shows one such tree chosen for the two-loop example.Here ∆ = 2, so we have removed two edges: J 3 and J 5 .From this reference tree we can derive ∆ other distinct spanning trees using the following method: 1) Take one of the ∆ edges that were removed to get the reference tree, and add it back to the graph.2) This creates a loop in the graph, known as a fundamental cycle (highlighted in green in Fig. 3C) [57].3) Remove one of the other edges in that loop (not the one just added), such that the graph returns to being a spanning tree.This new tree is distinct from the reference because it contains one of the ∆ edges not present in the reference tree.For example, in the top right of Fig. 3C, we added back edge 3, forming the fundamental cycle on the left states and E = 5 edges.A) The black arrows correspond to entries in the transition matrix Ω(λt): transition rates k ± i (λt) that depend on an external protocol λt.B) The red arrows labeled α correspond to the oriented stationary currents Jα(λt), defined in Eq. ( 13).C) On the left, one of the spanning trees of the graph, chosen to be a reference for constructing the tree basis.Edges deleted to form the tree are shown in faint red.On the right, two trees in this set derived from the reference one.Each such derived tree has a one-to-one correspondence with a fundamental cycle of the graph (highlighted in green).
loop.We then delete edge 1 from this loop, creating spanning tree 2. A similar procedure is used to construct tree 3.
We denote the ∆ + 1 trees (one reference + ∆ derived trees) constructed in this manner as the tree basis.We will label the trees in the basis set with γ = 1, . . ., ∆ + 1, where γ = 1 corresponds to the reference.In general, this basis is a subset of all possible trees, since T ≥ ∆+1.To every tree in the basis, we will associate a CD solution as follows.Let δJ (γ) (t) be a current difference vector that satisfies Eq. ( 19), but with the constraint that at every edge α that is not present in the γth tree, we have δJ (γ) α (t) = 0. We call this a fixed current constraint, since it corresponds to not being able to perturb the current associated with that edge via external control parameters.For example, imposing the restriction Ω ij = Ω ij and Ω ji = Ω ji for the pair (i, j) associated with edge α would make make δJ (γ) α (t) = 0. To find δJ (γ) (t), consider the (N −1)×E-dimensional reduced incidence matrix ∇ of the original graph; for example, Eq. ( 21) in the case of the two-loop graph of Fig. 3B.For a given spanning tree γ, we can construct an (N −1)×(N −1) submatrix ∇ (γ) from ∇ by choosing the N − 1 columns in ∇ that correspond to edges present in γ.This submatrix ∇ (γ) is equal to the reduced incidence matrix of the spanning tree γ.Hence we know that it has rank N − 1 and there exists an inverse [ ∇ (γ) ] −1 .Let us now construct a "stretched inverse": an S where the rows are populated by the following rule.If the row corresponds to one of the ∆ edges that was removed from the original graph to get the tree γ, it is filled with zeros; otherwise, it is filled with the corresponding row of [ ∇ (γ) ] −1 .For the three trees in Fig. 3C, labeled γ = 1, 2, 3 clockwise from left, the matrices [ ∇ (γ) ] −1 S have the following form: Moreover, it turns out one does not have to explicitly write down or invert ∇ (γ) in order to find the elements of [ ∇ (γ) ] −1 S .We can take advantage of a known graphical procedure for constructing inverse reduced incidence matrices of connected tree-like graphs [59,60].To determine the ith column of the matrix [ ∇ (γ) ] −1 S , start at the reference state (the state removed when constructing the reduced incidence matrix ∇, which in our case is always state N ).Among the edges of the spanning tree γ, there is a unique path that connects state N to state i.Following that path, if you encounter the current arrow J α oriented in the direction of the path, put a +1 in the row of [ ∇ (γ) ] −1 S corresponding to J α .Similarly if the current arrow is oriented opposite to the path, put a −1.All other entries in the ith column (current arrows not on the path, or not in the spanning tree) are set to zero.For example, consider the second column of [ ∇ (1) ] −1 S in Eq. ( 22).This corresponds to the path from state N = 4 to state 2 in the tree on the left of Fig. 3C.This includes edges 4 and 2, with the arrows along those edges all oriented opposite to the path.Hence the column has a −1 at the 4th and 2nd rows, and all other entries are set to zero.
By construction, each matrix [ We can now write down a solution for δJ (γ) (t), If we act from the left on both sides by ∇, we see that this form satisfies Eq. ( 19).The αth row of of [ ∇ (γ) ] −1 S is zero if edge α corresponds to a fixed current constraint (edge not present in the tree γ).Thus δJ Not only do the vectors δJ (t) associated with the tree basis constitute ∆ + 1 solutions to Eq. ( 19), they are also linearly independent from one another.To see this, note that because of the procedure to construct derived trees (adding back a distinct edge that was removed in the reference tree), a tree with γ ≥ 2 will have non-zero entry in δJ (γ) (t) in a position where every other tree (reference or derived) has a zero because of constraints.Hence the δJ (γ) (t) vector for each derived tree is linearly independent from all the other vectors in the basis.
We also know that any linear combination of solutions to Eq. ( 19) can be scaled by an overall normalization factor (to make the coefficients sum to one) so that it is also a solution to Eq. ( 19).Hence the following linear combination of basis vectors is a valid solution: Here w γ (t) are any real-valued functions where ∆+1 γ=1 w γ (t) = 1 at each λ t and λt .As we argue in the next section, the tree basis is complete: any CD solution δJ (t) can be expressed in the form of Eq. ( 24).Note that Eq. ( 20) is a special case of Eq. ( 24).When the original graph is tree-like, ∆ = 0 and there is only one spanning tree (γ = 1), equivalent to the original graph.In this case [ ∇ (1) ] −1 S = ∇ −1 and the sole coefficient function w 1 (t) = 1 by normalization.

C. Completeness of the tree basis
To prove that any CD solution can be expressed as a linear combination of tree basis solutions δJ (γ) (t), let us first introduce ∆ vectors of the following form: for γ = 2, . . ., ∆ + 1.Since both basis vectors on the right-hand side of Eq. ( 25) satisfy Eq. ( 19), we know that Hence V (γ) (t) is a vector in the null space of ∇.Moreover since the basis vectors δJ (γ) (t) are linearly independent, the set V (γ) (t) constitutes ∆ linearly independent null vectors of ∇.We can find the dimension of the null space, nullity( ∇), using the rank-nullity theorem: rank( ∇) + nullity( ∇) = E, where E is the number of columns in ∇.Since rank( ∇) = N − 1 for a connected graph, as described earlier, we see that nullity( ∇) = E − (N − 1) = ∆.Thus the ∆ linearly independent vectors V (γ) (t) span the whole null space.If there existed a vector δJ (t) that satisfied Eq. ( 19) but could not be expressed as a linear combination of basis vectors, then the corresponding vector V(t) = δJ (t) − δJ (1) (t) would be a null vector that is linearly independent of all the V (γ) (t).But since the latter span the whole null space, this is impossible.Hence every CD solution δJ (γ) (t) satisfying Eq. ( 19) must be expandable in the form of Eq. ( 24).

D. General solution in the cycle basis
The discussion in the previous section also allows us to rewrite the expansion in Eq. ( 24) in an alternative form that is convenient in practical applications.Using the fact that ∆+1 γ=1 w γ (t) = 1, Eq. ( 24) can be equivalently expressed as Since the vectors V (γ) (t) form a basis for the null space of ∇, the second term in the last line of Eq. ( 27), with its arbitrary coefficient functions w γ (t), is general enough to describe any vector function in the null space.With no loss of generality, we can rewrite this second term in another basis for the null space instead.A convenient choice is the fundamental cycle basis corresponding to some reference spanning tree (we need not choose the same reference as used to find δJ (1) (t)).The ∆ fundamental cycles were identified in the procedure to construct derived trees.If we assign an arbitrary orientation to the cycles (clockwise or counterclockwise), then the E-dimensional cycle vector c (γ) , associated with the derived tree γ + 1, is defined as follows: a ±1 at every row whose corresponding edge in the original graph belongs to the fundamental cycle, with a +1 if the edge direction is parallel to the cycle orientation, −1 if anti-parallel.All edges not belonging to the fundamental cycle are zero.For the reference tree in Fig. 3C the fundamental cycles are highlighted in green on the right of the panel.Here the two cycle vectors are: In general, the ∆ fundamental cycle vectors form a basis for the null space of ∇ [57].
In terms of the cycle vectors, Eq. ( 27) can be written as where v γ (t) for γ = 1, . . ., ∆ are another set of arbitrary coefficient functions.The convenience of Eq. ( 29) over Eq. ( 24) is that we only need to find one spanning tree solution δJ (1) (t).Both have the same number of degrees of freedom: in the first case ∆ coefficient functions w γ (t) for γ = 2, . . ., ∆ + 1 (since w 1 (t) depends on the rest through the normalization constraint); in the second case ∆ coefficient functions v γ (t).
Finally we note that because of Eq. ( 15), the oriented stationary current vector J (t) corresponding to the original protocol is in the null space of ∇.Hence it can also be expanded in terms of the cycle vectors as with some coefficient functions u γ (t).Since the CD currents J (t) = J (t) + δJ (t), we can combine Eqs. ( 29) and ( 30) to get the most general expression for any set of currents that satisfies Eq. ( 18): Here Φ γ (t) ≡ u γ (λ t ) + v γ (t).Because the v γ (t) are arbitrary, the functions Φ γ (t) are also arbitrary, and we still have the same ∆ degrees of freedom to span the solution space.Combining Eq. ( 31) with Eq. ( 23), we can write the general solution for the CD currents as where C is the E × ∆ matrix whose columns are the fundamental cycle vectors c (γ) , γ = 1, . . ., ∆, and Φ(t) is a ∆-dimensional vector whose components are the arbitrary functions Φ γ (t).

E. Solving for the CD transition rates
Once we know J (t) from Eq. ( 32), the final step to find the CD protocol is solving for the CD transition rates via Eq.(17).Using the fact the diag(x)y = diag(y)x for any vectors x and y, Eq. ( 17) can be rewritten as where . Since every row of ∇ ± T has exactly one nonzero element equal to 1, and we focus on systems where the elements of ρ(λ t ) are all positive, the matrix M ± (t) has positive elements on the diagonal and hence is invertible.We can thus solve Eq. ( 33) for k + (t), Any vectors k + (t) and k − (t) with positive elements that satisfy Eq. ( 34) constitute valid forward/backward CD transition rates.If one can control both k + (t) and k − (t), valid solutions to Eq. ( 34) always exist.
A special case of Eq. ( 34) often arises in biological contexts: for each edge α, we can only externally control one rate, which we take to be the forward rate without loss of generality.The backward rates do not vary in the CD protocol, k − (t) = k − .In this scenario, when we use Eq. ( 34) to solve for k + (t), we have to choose the arbitrary function vector Φ(t) in Eq. (32) to ensure that k + (t) has positive elements.If this is not possible then CD driving along the given target trajectory cannot be achieved by only changing the forward rates.

F. Thermodynamic costs of CD driving
To quantify the thermodynamic costs of driving via the CD protocol, one can calculate the total entropy production rate Ṡtot (t) at time t [56], where the E-dimensional edge affinity vector χ(t) is given by Throughout this section we will use the convention that a function like ln v or exp v for a vector v is a vector with components (ln The structure of Eq. ( 35), where every term in the inner product is of the form (x−y)(ln x−ln y) for non-negative quantities x and y, gives Ṡtot (t) ≥ 0, in accordance with the second law of thermodynamics.Plugging Eq. ( 34) into Eq.( 36) we can express χ(t) as The form of Eq. ( 37) has an interesting consequence.Imagine a hypothetical scenario where one can control both the forward k + (t) and backward k − (t) rates, satisfying Eq. (34).For a given CD current solution J (t), the limit where the backward rates become large, k − α (t) → ∞ for all α, would correspond to χ(t) in Eq. ( 37) becoming arbitrarily small, | χ α (t)| → 0. Note that due to Eq. (34) this limit also means the forward rates k + α (t) → ∞.With fixed J (t) and vanishing χ(t), we see from Eq. ( 35) that the instantaneous entropy production Ṡtot (t) can approach zero if both backward and forward rates can be made arbitrarily large.
The fact that we can drive the system over a trajectory in finite time with negligible thermodynamic cost is only possible because increasing the rates corresponds to making the local "diffusivity" in the system large (if we imagine dynamics on the network as a discrete diffusion process).In other words we are reducing the effective friction to zero in order to eliminate dissipation.In practice this extreme limit is not realistic, particularly for biological systems.There are likely to be physical constraints that prevent us from simultaneously tuning each pair of rates in the network over an arbitrary range, so the CD implementation with Ṡtot (t) → 0 is not realizable.
If there is some constraint on at least one of the rates in each pair, the minimum Ṡtot (t) among all CD protocols for a given target trajectory will have a finite value.For example, let us take the case mentioned above where the backward rates are fixed, k − (t) = k − .Let us also assume that ∆ > 0, so that varying the functions Φ(t) in Eq. ( 32) gives all possible protocols for a particular choice of ρ(λ t ).Among these protocols, the condition for minimizing Ṡtot (t) is given by a gradient with respect to Φ(t) of Eq. ( 35), where 0 and 1 denote vectors of zeros and ones respectively (of size ∆ and E in this context).To derive this, we have used the fact that ∂ J α (t)/∂Φ γ = C αγ from Eq. ( 32), and that J (t) = Γ −1 (t) e χ(t) − 1 from Eq. ( 37), where Note that after finding a set of Φ(t) that satisfies Eq. ( 38), one must also check that the corresponding forward CD rates k + (t) given by Eq. ( 34) all have positive elements.
There is one case where Eq. ( 38) can be solved analytically.If the driving is slow enough that the magnitudes of the components of J (t) are small (and hence also those of χ(t) via Eq.( 37)), then χ(t) ≈ Γ(t) J (t).In this limit Eq. ( 38) becomes From Eq. ( 39) and (32) we can now solve for Φ(t), where The corresponding minimum instantaneous entropy production for this slow driving regime is: where L(t) = I E + CB(t), and I E is the E × E identity matrix.
The protocol that satisfies Eq. ( 39) has a particular interpretation: let us define a ∆-dimensional cycle affinity vector χ cyc (t) = C T χ(t), where each component is the sum of the edge affinities over a fundamental cycle [61].Equation (39) thus states that χ cyc (t) = 0. Satisfying these ∆ conditions is equivalent to specifying that the CD protocol rates obey local detailed balance, in other words that there exists a right singular vector of G(t) with singular value zero.Hence in the slow driving regime, out of all the possible CD protocols for a given target, it is the one with local detailed balance that minimizes Ṡtot (t).This is not generally true outside of the slow driving regime, where Eq. ( 38) does not coincide with χ cyc (t) = 0.This observation is compatible with the results of Ref. [62], where they showed among all driving protocols between fixed initial and end distributions the one with the lowest total entropy production violates detailed balance.Equation (38) shows that this is also true for any specific target trajectory.Note that when ∆ = 0 for a system with a tree-like graph, there is no degeneracy in J (t) and the CD protocol automatically satisfies local detailed balance.

II. GENERALIZING COUNTERDIABATIC DRIVING: NON-STATIONARY TARGETS AND LOCAL CONTROL
Up to now the theory has been framed in terms of driving all N states in the system along a target ρ(λ t ) that is an instantaneous stationary distribution (Eq.( 3)) with respect to some original transition matrix Ω(λ t ).However in biological contexts the control problem may be more general: perhaps one is interested in having only a subset of states follow a target trajectory, and the target trajectory does not necessarily have to be a stationary one.In this section we generalize our theoretical framework to accommodate both these aspects.Doing so allows us to define a set of graphical rules for controllability in discrete-state Markov models, which we will apply to biological examples in Sec.III.

A. Driving along non-stationary target trajectories
Let us imagine an arbitrary target trajectory ρ(t) where Eq. ( 3) is not satisfied, and hence Ω(λ t )ρ(t) = 0 for at least some t during the driving time interval 0 to τ .A special case of this is known as the fast-forward problem, where we start and end in a stationary distribution, Ω(λ 0 )ρ(0) = Ω(λ τ )ρ(τ ) = 0, but allow a nonstationary trajectory for 0 < t < τ .We seek a modified transition matrix Ω(t) that satisfies the master equation ∂ t ρ(t) = Ω(t)ρ(t) during driving.
It turns out that our framework generalizes to arbitrary ρ(t) in a straightforward way: the main difference is that we no longer enforce Eq. ( 15), and we seek a solution to Eq. ( 18) instead of Eq. (19).The general form of J (t) that satisfies ∂ t ρ(t) = ∇ J (t) is given by: Here [ ∇ (1) ] −1 S is the stretched inverse corresponding to one of the spanning trees of the graph, and Φ (t) is a ∆-dimensional vector whose components are arbitrary functions Φ γ (t).From the fact that ∇[ ∇ (1) ] −1 S = I N −1 and ∇C = 0 (since the columns of C are cycle vectors) we see that Eq. ( 42) does indeed solve Eq. ( 18).Once J (t) is known, we find associated transition rates k ± (t) that satisfy Eq. (34).In fact, because the structure of Eq. ( 42) is identical to our earlier CD current expression in Eq. ( 32), the final form of the control protocol solution is the same whether or not the target is stationary.
In order to elucidate the perturbation δJ (t) to our original currents J (t) necessary achieve the driving, we can write it in the form where the new vector Φ(t) is defined via The right-hand side of Eq. ( 44) vanishes when acting on it with ∇ from the left, and hence it is in the null space of ∇.Thus it can be expressed as a linear combination of the fundamental cycle vectors, and hence there must exist a Φ(t) that satisfies Eq. ( 44).Note that [ ∇ (1) ] −1 S ∇ is an E × E matrix that does not equal the identity in general, since [ ∇ (1) ] −1 S is a right, not left, pseudoinverse of ∇.The one case where [ ∇ (1) ] −1 S ∇ = I E is when the original graph is a tree and hence E = N − 1.
For the CD driving scenario, where the target trajectory is stationary and ∇J (t) = 0, Eq. ( 43) reduces to our earlier CD solution in Eq. ( 29).More generally, the structure of Eq. ( 43) allows us to see that only a subset of the currents in the original system need to be modified in order to achieve an arbitrary target.Since Φ(t) is arbitrary, we can set Φ γ (t) = 0 for all γ, and hence from Eq. ( 43) we get that δJ α (t) = 0 only for those N − 1 edges α present in the spanning tree (because the rows of [ ∇ (1) ] −1 S associated with edges not in the tree are all zero).If k ± α (t) are the original forward/backward transition rates associated with edge α, we have to be able to modify one or both of them to new rates k ± α (t) in order to satisfy the δJ α (t) condition, Eq. (43).We call an edge α where it is possible to modify the transition rates via external parameters a controllable edge.No solution exists with less than N −1 controllable edges, and if we set Φ γ (t) = 0 we generally get solutions that require more than N − 1 controllable edges, since the cycle vectors involve currents on edges not in the tree.Since we can use any spanning tree in Eq. ( 43), we can formulate a general rule for global control, the ability to drive every state in the network along an arbitrary target: Global control condition: in order to drive an N -state network along an arbitrary target trajectory ρ(λ t ), the set of controllable edges must span the entire network graph.One consequence is that global control is impossible with less than N −1 controllable edges.
The minimal condition for global control (N − 1 controllable edges forming a spanning tree) is the same whether or not the trajectory is an instantaneous stationary one.Depending on the physical details of a specific system, there may be additional conditions necessary to achieve global control, but the above one must always be fulfilled.For example, if the backward rates cannot be modified, the forward rates k + (t) given by Eq. ( 34) must all be non-negative.A similar story applies to the case where different controllable edges cannot be independently varied (i.e. because they depend on a single external parameter).In this situation it may not be possible to simultaneously satisfy the δJ α (t) conditions at all controllable edges.

B. Local control
The local control problem means we are only interested in having N T < N − 1 of the states in the system follow a target trajectory.If we label the states such that the first N T are target states, we need to find Ω(t) such that the master equation ∂ t p(t) = Ω(t)p(t) is satisfied with a solution whose probability vector has the form where we denote the N T -dimensional target trajectory vector ρ(t) and the (N − N T )-dimensional vector of remaining non-target states as π(t).As discussed below, the case where N T = N − 1 corresponds to the global control scenario, since the probability of the N th state is constrained by the normalization condition The target trajectory ρ(t) is specified beforehand, and we are looking for all possible solutions compatible with a given ρ(t).
Imagine there are E c ≤ E edges in the system that are controllable, and we label the edges such that α = 1, . . ., E c correspond to the controllable ones.We can then partition the current vector J (t) into controllable (supserscript c) and not-controllable (superscript n) components as follows: Analogously, we can partition the forward/backward rates on the edges: Note that by definition the rates k ±n on the noncontrollable edges cannot be externally modified.We additionally assume here that they are time-independent, the typical case for biological systems.This allows us to find the analytical solution for π(t) shown below.If k ±n (t) are fixed, time-dependent functions, Appendix A shows how the general solution requires numerically solving a differential equation for π(t).Finally, we can partition the graph matrices ∇ and G(t) each into four submatrices as shown in Fig. 4A.At the outset of the local control problem, the following quantities are known: the target trajectory ρ(t), the non-controllable edge rates k ±n , the four submatrices of ∇, and the bottom two G(t) submatrices, G nρ and G nπ (which depend only on the non-controllable rates, and hence are time-independent).The goal is to figure out controllable edge rates k ±c (t) that force to system to follow the trajectory on the target states.A detailed derivation of the solution is presented in Appendix A.
Here we summarize the resulting procedure, which consists of three steps.i) Solve for π(t): where Here π(0) is a set of arbitrary initial probabilities for the non-target states, with the constraint that the com-ponents of π(0) and ρ(0) must add up to 1. C c is a matrix with dimensions E c × ∆ c whose columns form a basis for the null space of ∇ ρc .As we will argue below, ∆ c = E c −N T when a local control solution exists.Φ c (t) is a ∆ c -dimensional vector of arbitrary functions.The freedom to choose π(0) and Φ c (t) means that in general the solution for π(t) is non-unique.
A local control solution is not always possible, since there are two criteria that need to be satisfied.The first is that there must exist an The second is that the components π(t) from Eq. ( 48) need to be valid probabilities, π i (t) ≥ 0 for all t during the driving protocol.Note that normalization, where the components of π(t) and ρ(t) sum to 1 at all t, is guaranteed by the structure of the solution, but π i (t) ≥ 0 is not automatically enforced.Given the ability to choose π(0) and Φ c (t), it is often feasible to satisfy this second criterion.For the first criterion to work, ∇ ρc must have rank N T , and there is a simple graphical method to check for this.Let us define a target subgraph as follows: starting from any state in the target subset, this is the connected subgraph of all states reachable via only controllable edges.Examples of target subgraphs are highlighted with dashed curves in Fig. 4B,C.These two panels depict the same network and same set of E c = 6 controllable edges, but with two different sets of target states with N T = 4.There may be multiple target subgraphs in a network, involving disjoint subsets of the controllable edges.As shown in the figure, a target subgraph must include at least one target state, but can also include non-target states.We can now formulate the general rule for local control, encapsulating both criteria: Local control condition: in order to drive N T < N − 1 states from a network of size N through an arbitrary trajectory ρ(t), every target subgraph must include at least one non-target state.One consequence is that local control is impossible with less than N T controllable edges.An additional criterion is that there must be a solution for π(t) with non-negative components during the time interval of driving.
The reason that the subgraph condition is sufficient for [∇ ρc ] −1 S to exist is that we can choose one of the nontarget states in each target subgraph as the analogue of a reference state for that subgraph.This makes the row of ∇ ρc for each target state equal to the row of the reduced incidence matrix for the target subgraph to which the state belongs.We know that rows of a reduced incidence matrix are linearly independent from one another for the same subgraph, and for different subgraphs they are linearly independent because they involve different subsets of edges.Thus overall if the local control condition is satisfied, all rows of ∇ ρc are linearly independent.The condition also implies a lower bound on the number of controllable edges, E c ≥ N T .To see this, let us imagine there are K target subgraphs in the network, κ = 1, . . ., K each containing n κ target states and at least one non-target state.For each subgraph to be connected, it must involve at least n κ controllable edges.We thus get that E c ≥ K κ=1 n κ = N T .The linear independence of the N T rows of ∇ ρc , plus the fact that the number of columns E c is at least N T , guarantees that ∇ ρc has rank N T .Since the rank and nullity of ∇ ρc must sum to E c , this also implies that its nullity ∆ c = E c − N T .
When the local control condition is satisfied, we can write down the stretched inverse [∇ ρc ] −1 S using a graphical procedure analogous to the global control case.First, choose a controllable edge spanning tree for each target subgraph.To find the ith column of [∇ ρc ] −1 S , follow the tree path from the non-target reference state to state i in the corresponding subgraph, and put a ±1 at each row α where the corresponding edge is parallel / anti-parallel to the path.All other entries in the column are zero.

III. DRIVING IN BIOLOGICAL NETWORKS
To illustrate the general theory in specific biological contexts, we consider two examples of driving in biochemical networks, corresponding to global and local control respectively.The first example is a simple genetic regulatory switch involving a repressor protein and corepressor ligand binding to an operator site on DNA, turning off the expression of a set of genes.Here it turns out there are enough control knobsconcentrations of repressors, corepressors, and repressorcorepressor complexes-to implement a whole family of exact global control solutions.Among this family we can then examine which ones satisfy certain physical constraints, or minimize thermodynamic costs.The second example involves a chaperone protein that binds to a misfolded substrate, catalyzing the unfolding of this misfolded protein and giving it another chance to fold into the correct ("native") state.The available control knobs-chaperone and ATP concentrations-are insufficient for global control, but do allow the system to locally control the probability of being in the misfolded state.This local control turns out to be of crucial importance, since rapidly decreasing the misfolded probability is a way to ameliorate the damage due to heat shock.In fact the local control protocols from our theory qualitatively mimic experimental results from yeast and E. coli.

A. Repressor-corepressor model
The first system we consider is a common form of gene regulation in bacteria, illustrated schematically in Fig. 5A: a repressor protein has the ability to bind to an operator site on DNA.When bound, it interferes with the ability of RNA polymerase to attach to the nearby promoter site, preventing the transcription of the genes associated with the promoter.The system acts as a genetic switch, with the empty operator site the "on" state for gene expression, and the occupied operator site the "off" state.In many cases, additional regulatory molecules-inducers or corepressors-influence the binding affinity of repressor proteins [63].In the present model, binding of the bare repressor to the operator site is weak (it unbinds easily), but the binding strength is enhanced in the presence of a particular small molecule-the corepressor.Hence the corepressor concentration acts like an input signal, with sufficiently high levels leading to the promoter site being occupied with high probability, and the associated genes being turned off.
Such genetic switches are basic building blocks of natural and synthetic biological circuits.From the control standpoint, can we drive the switch through a prescribed trajectory, turning it on or off in a finite time, and at what cost?
There are several reasons this system provides a convenient testing ground for our theory.As described below, it can be modeled with three discrete states connected via Markovian transitions, forming a three-state loop (Fig. 5A).This is the simplest graph structure where there exists a whole family of global control protocols for any given target trajectory.Though each of these protocols achieves the same target, they are chemically and thermodynamically distinct, allowing us to explore interesting facets of degeneracy in the driving theory.Since all the forward and reverse rates of the system are known experimentally, taken from in vitro measurements of the purine repressor (PurR) system of E. coli [63,64], the system also provides a simple platform to directly test theoretically predicted control protocols in the future (for example using a time-resolved version of the in bound to a represssor-corepressor complex.Transition rates between the states are shown in green.The binding reaction rates depend on three concentrations of molecules in solution: R(t) for bare repressors, C(t) for corepressors, and X(t) for the complexes.B) One of the spanning trees for the associated network graph, with the edge deleted to form the tree shown in faint red.We take this to be the reference spanning tree for the tree basis.C) The other tree in the basis, with the corresponding fundamental cycle in green.
vitro fluorescence spectroscopy already successfully applied to PurR in Ref. [64]).Finally, the general structure of the network, with edges where either the forward or backward transition depends on the concentration of a regulatory molecule or enzyme, is quite representative of biochemical systems in general.Hence it serves as a jumping off point for the analysis of more complex biological networks.
As is generally the case with genetic regulation in biology, the processes underlying repressor dynamics are stochastic [65].The three discrete states in our Markov model are: 1) free operator; 2) bare repressor bound to the operator; 3) repressor-corepressor complex bound to the operator.Transitions in both directions (clockwise and counterclockwise in Fig. 5A) are possible.In each pair of transition rates between neighboring states there is a binding reaction proportional to the concentration of a chemical species in solution.The relevant concentrations are those of bare repressors R(t), corepressors C(t), and repressor-corepressor complexes X(t).If we label the binding reactions as the forward rates, then k + (t) = (k r R(t), k c C(t), k x X(t)), with associated binding constants k µ for each species µ.The concentration dependence of the forward rates means we have three controllable edges, thus satisfying the global control cri- tocols that all drive the system along the target trajectory ρ(λt).The four protocols are: spanning tree 1 (violet, corresponding to Fig. 5B); spanning tree 2 (teal, corresponding to Fig. 5C); the solution satisfying local detailed balance, ∆µ(t) = 0 at all t (yellow); and the optimal solution that minimizes Ṡtot (t) at all t (thick black).For each case we depict: B) the total entropy production rate Ṡtot (t), with the inset showing the difference δ Ṡtot (t) ≡ Ṡtot (t) − Ṡtot,opt (t) between non-optimal and optimal rates on a log scale [units: kB/min]; C) the instantaneous chemical potential ∆µ(t); D) the concentrations C(t), R(t), X(t).
terion (at least 2 controllable edges for a 3-state network).The reverse rates, describing the unbinding reactions, cannot be externally tuned: To be concrete, we choose parameters based on the purine repressor (PurR) system of E. coli [63,64].The PurR protein turns off genes responsible for the de novo production of purines, a class of molecules including guanine and adenine that are essential ingredients in DNA/RNA and energy transducing molecules like ATP.If the cell has an excess of purines (for example from environmental sources), this is signaled by an abundance of the corepressors guanine or hypoxanthine (a purine derivative) that form complexes with PurR, enabling it to bind strongly with its operator site.This way, the cell can switch off the energetically expensive de novo production of purines when it is not needed.The parameter values, as well as the calculation of the control protocols, are shown in Appendix B.
For a given target trajectory ρ(t), the goal is to find forward rates k + (t) = (k r R(t), k c C(t), k x X(t)) that force the system to be on-target.These define concentration protocols R(t), C(t), and X(t) for the three species.Fig. 6A shows our chosen target ρ(t), mimicking a biological scenario where the genetic switch is rapidly turned off over the course of a couple of minutes: the free operator (state 1) probability is decreased, with a corresponding increase in the repressor-bound states 2 and 3.This particular ρ(t) (details in Appendix B) consists of instantaneous stationary distributions for the system, so the driving is CD, but as described in Sec.II A any other ρ(t) could have been chosen.
Since ∆ = 1 for the oriented current graph, there will be many possible concentration protocols that drive the system along exactly the same target, via different choices of Φ 1 (t) in Eq. ( 32) for the CD currents.Fig. 6D shows four different examples of concentration protocols.The violet and teal curves are for the tree 1 and 2 solutions (Fig. 5B,C) where J 1 (t) = 0 and J 2 (t) = 0 respectively.The other two protocols are described below.Despite leading to the same system behavior, the protocols have quite distinct physical characteristics, with concentrations varying up to an order of magnitude among the four examples shown.They also differ in the cycle affinity χ cyc (t), which is a scalar since ∆ = 1.This affinity has a more direct physical interpretation as the chemical potential ∆µ(t) = k B T χ cyc (t) for the repressor-corepressor binding reaction, and is plotted in Fig. 6C for the four protocols.One of the protocols (yellow curve) has rates chosen so that ∆µ(t) = 0, satisfying the local detailed balance condition.As described in Sec.I F, we know that this protocol should be the one with the smallest entropy production rate in the slow driving limit.Here we are away from that limit, but the protocol still does well in economizing thermodynamic costs, as seen in the plot of Ṡtot (t) in Fig. 6.The protocol that satisfies Eq. ( 38), and hence optimizes Ṡtot (t), is shown as a thick black curve for comparison.It is close, but not exactly equal to, the detailed balance protocol, exhibiting slightly negative ∆µ(t) at intermediate times (Fig. 6C).The inset of Fig. 6B shows the difference δ Ṡtot (t) ≡ Ṡtot (t) − Ṡtot,opt (t) between each non-optimal solution and the optimal one.The detailed balance solution is significantly closer to optimal entropy production than the two tree solutions.
The repressor-corepressor model illustrates the variety of physically realizable control solutions that can exist in certain cases.This gives nature (or an experimentalist engineering a synthetic system) a rich array of options to achieve a specific probabilistic target.When we observe, for example, a genetic switch within a biological circuit being rapidly turned off by changing concentrations of external species, there will typically be a variety of alternatives that would have led to same state distribution at each instant of time.An interesting question for future studies would be to ask whether certain options would be evolutionarily favored over others because of selection pressures due to energetic costs [66].

B. Chaperone model
Many newly synthesized proteins, susceptible to misfolding, become trapped in long-lived metastable states that are prone to aggregation.Since aggregates present a danger to the survival of the cell, there exists an elaborate rescue machinery of molecular chaperone proteins that facilitate unfolding or disaggregating misfolded proteins [7][8][9][10].In the case of E .coli, which has the most extensively studied chaperone network, certain components like the GroEL-GroES system are obligatory for survival [67].Environmental stresses further exacerbate the problem, and an increase of ambient temperature by just a few degrees can significantly enhance protein misfolding and consequently aggregation [11].Responding to a heat shock requires creating extra capacity, since even under normal conditions the majority of chaperones are occupied by misfolded proteins [9] (i.e.occupancy for GroEL can approach 100% for fast-growing E. coli [10]).This is accomplished by rapidly upregulating the number of chaperones to cope with additional misfolded proteins [11,12].
Most chaperones require constant power input in the form of ATP hydrolysis.As a result the stationary probability distribution of conformational states for a protein interacting with a chaperone will generally be out of equilibrium (non-Boltzmannian) [68,69].When the chaperone concentration increases after a heat shock (for example following a sudden rise to a new temperature [70]), the protein is driven away from the previous stationary distribution, and eventually relaxes to a new stationary distribution once the chaperone concentrations reach steady-state values at the new temperature.Chaperone upregulation during heat shock therefore serves as a natural example of nonequilibrium driving in a biological system.Fig. 7 shows experimental results for the heat shock response of two representative organisms.In Fig. 7A relative mRNA levels for six different chaperone genes in S. cerevisiae yeast are plotted as function of time [71].Since higher mRNA expression generally leads higher concentrations of the proteins coded for by the mRNA, the mRNA levels can be seen as a proxy for chaperone concentration.At the start of the experiment, the temperature is raised from 30 • to 39 • C, and then held constant.Chaperone gene expression rises sharply in the first half-hour (in some cases by more than an order of magnitude), then peaks and levels off at a value roughly half that of the peak.E. coli shows similar behavior (Fig. 7) for mRNA levels of the dnaK chaperone gene [70], in this case following a heat shock from 30 • to 42 • C. In both organisms the chaperone levels overshoot and then remain elevated for a long duration after the shock, a characteristic feature of the heat shock response [11,72].Interestingly, E. coli shows another, less common, behavior: ATP concentration transiently increases by about a factor of two in the first minutes after the shock, an observation additionally supported by metabolic evidence [73].How do such changes in chaperones and ATP affect the state distribution of a protein targeted by chaperones?In the analysis below, we will see that these two control knobs enable local control of states involving the misfolded protein.
Our starting point is a four-state Markov model for chaperone-assisted protein unfolding, inspired by earlier models like those of Refs.[68,69].We focus on a network of four states for a particular substrate ("client") protein, and one type of chaperone, depicted in Fig. 8A: 1) a misfolded protein state, prone to aggregation; 2) the misfolded protein bound to chaperone; 3) an intermediate conformational state of the protein, along the folding pathway between the unfolded and native states; 4) the correctly folded "native" state.These four states can interconvert with transition rates denoted in the figure (further details below).The model is a small biochemical module within a broader set of processes, some of which are depicted schematically with dashed arrows in the figure: protein synthesis and the initial folding to the intermediate state, and aggregation of the misfolded proteins.Our focus will be on a single protein once it enters the intermediate state, and then transitions among the four states.Similarly we ignore the aggregation process, occurring over much larger timescales than the transitions in the network.We model the dynamics in the aftermath of a heat shock [11,72]: a sudden jump to some high temperature T , which then remains fixed as the system adapts.The conditions favor misfolding over the native folding pathway.In the absence of chaperones, state 2 (misfolded) would be most likely, and over longer timescales this would eventually result in a build-up of aggregates.
To understand how this system can be controlled, let us summarize the various transitions in the network (Fig. 8A).The protein can interconvert between states 3 and 1 with rates k m and k −m .A chaperone can bind to the misfolded protein at rate k c C(t), where C(t) is the concentration of unoccupied chaperones and k c is the binding constant.Once bound, the chaperone cat- FIG. 7. Experimental results for the heat shock response of yeast and E. coli.A) Relative mRNA expression of six chaperone genes in S. cerevisiae yeast versus time after a the start of a 30 • → 39 • C heat shock.The genes, listed on the right, all have the Gene Ontology database annotation 0051082 [74,75], indicating that their products exhibit chaperone activity (binding to unfolded proteins).B) The black curve shows the relative mRNA expression of the E. coli chaperone gene dnaK versus time after the start of a 30 • → 42 • C heat shock [70].The blue curve shows ATP concentration for the same system.
alyzes the partial unfolding of the misfolded state to the intermediate state at rate k a .This conversion may involve several substeps and hydrolysis of multiple ATP molecules, but we simplify the process to a single reaction step hydrolyzing one ATP molecule, with some rate function k a .Though typically negligible compared to the forward rate k a , the reverse rate k −a C(t), proportional to chaperone concentration, must be formally defined in order to have a thermodynamically complete description of the system.Transitions from the intermediate to native state occur with rate k n , and from the native to misfolded state with rate k u .The full details of the model, including parameter values estimated from experimental data on the chaperone GroEL assisting the folding of the substrate protein MDH [68], are described in Appendix C. The model allows us to determine what kinds of external control are possible in the system.We will consider two different scenarios: the first is the typical case with one control knob, as seen in yeast, where chaperone concentration C(t) can vary as part of the heat shock response, but ATP levels (and thus k a ) are fixed.t) is negligible, there is effectively only one controllable edge in the network, the one between state 1 and 2. We can thus choose a target subgraph to enable one-state local control, with state 1 (misfolded) as our target.C-D) Example of a one-state local control solution, with a chaperone concentration protocol C(t) in panel C forcing the state 1 probability p1(t) to follow a target sigmoidal decrease in the misfolded probability ρ1(t) in panel D. The state probabilities pi(t), i = 1, . . ., 4 are calculated by numerical solution of the master equation, with parameters described in Appendix C. E-G) Same as panels B-D except with an additional controllable edge due to varying ATP levels, allowing ka(t) to be time-dependent.We can now have two-state local control, targeting states 1 and 2, and find a protocol of C(t) and ka(t) to make the system follow a chosen set of targets ρ1(t) and ρ2(t).
In principle changes in C(t) affect two edges in the network, via the rates k c C(t) and k −a C(t).However k −a is usually so small that changes in C(t) make no noticeable difference to the current between states 2 and 3, which is dominated by the rate k a .For example in our parameter set k −a /k a = 0.0952 M −1 .Since chaperone concentrations are usually on the µM scale or smaller, k −a C(t) at least seven orders of magnitude smaller than k a .If k a is fixed, we can treat the edge between states 2 and 3 as effectively non-controllable.Hence in this scenario we have a network with a single controllable edge (edge 1 between states 1 and 2) and thus can choose a target subgraph as shown in Fig. 8B.The target state in the subgraph is state 1, and state 2 is included as the non-target state to fulfill the local control condition.We will dub this scenario one-state local control.Given the danger of having a high probability of misfolded proteins, state 1 is a natural target for control.One could imagine of course other biologically plausible targets, for example trying to control the probability of the protein being in its functional native state 4. However given the available control knob C(t), no target subgraph including state 4 would satisfy the local control condition, since there are no controllable edges involving state 4.This highlights the usefulness of the condition to help us rationalize the influence of external factors on the system.
To illustrate how one-state local control works, we choose a target trajectory for ρ 1 (t) that is a rapid sigmoidal decrease in the misfolded probability over the course of a few minutes, shown as the dotted curve in Fig. 8D.This would be a desirable heat shock response for the system.Details of the target and the local control solution are described in Appendix C. Fig. 8C shows the necessary chaperone concentration protocol C(t) needed to drive state 1 along the target, derived from Eq. (51).Note how C(t) closely resembles the qualitative features of actual chaperone expression in yeast experiments (Fig. 7A): a large initial increase, a peak, and then a more gradual leveling off.This kind of protocol is the typical solution if your target is a rapid suppression of the misfolded probability.If we plug in the C(t) result and numerically solve the master equation for the system, we find p 1 (t) (solid blue curve in Fig. 8D) agreeing exactly with ρ 1 (t), as expected.The remaining non-target state probabilities are compatible with Eq. ( 48) for π(t) = (p 2 (t), p 3 (t), p 4 (t)).As explained in Appendix C, the second clause in the local control condition, that all components of π(t) must be nonnegative during driving, actually gives us information about the kinds of target functions ρ 1 (t) that we can implement.Virtually any ρ 1 (t) that is monotonically decreasing leads to non-negative π(t) solutions, regardless of how steep the decrease.On the other hand, if one chose an increasing ρ 1 (t) one could violate the nonnegative π(t) criterion.This distinction agrees with our biological intuition: changing chaperone concentration gives us fine-grained control in decreasing the misfolded probability, via the outgoing transition k c C(t) from state 1.An arbitrary increasing ρ 1 (t) target would be biologically detrimental, and thus it is not surprising the system has not evolved the proper control knobs to achieve it.
In the second scenario we have two control knobs: in addition to C(t), we imagine ATP levels can change, as seen in E. coli, allowing k a (t) to be a time-varying function.This gives two independently controllable edges (again ignoring the negligible role of k −a ), and allows us to choose a target subgraph as seen in Fig. 8E.We fulfill the local control condition for two target states (1 and 2), and thus dub this scenario two-state local control.Despite the additional control knob, we still cannot include state 4 as a target.Having both ρ 1 (t) and ρ 2 (t) as targets allows one to control the total probability of observing a misfolded protein, both free (state 1) and bound to chaperone (state 2).For example one can avoid the peak for state 2 seen in the one-state local control results of Fig. 8D.Such transient accumulation of misfolded proteins bound to chaperones might not be ideal if the chaperones need to be turned over quickly to handle multiple different substrate proteins.Fig. 8G shows the same ρ 1 (t) target as in the one-state control case, but with ρ 2 (t) chosen to be a small sigmoidal step, avoiding transient accumulation (red dotted curve).The solutions for C(t) and k a (t) are shown in Fig. 8F.The chaperone concentration protocol is nearly identical, and ka (t) exhibits a transient peak, necessary to suppress the build-up of misfolded proteins on the chaperones.Comparing with the E. coli experimental results in Fig. 7B, we again have qualitative similarities in both chaperone expression and ATP concentration behavior.At least for our representative set of model parameters, the influence of the additional k a (t) control knob is not as large as that of C(t).This may in part explain why varying ATP levels is a fairly atypical heat shock response, while chaperone upregulation is universal.However the two-state local control case does clearly illustrate the increasing precision of influence available with additional control knobs.

IV. CONCLUDING REMARKS
Our theory of classical stochastic driving and its biological applications open up a variety of questions for future work.We have focused here on discrete state Markov models, but taking the continuum limit of such models allows one to connect to diffusive dynamics described by Fokker-Planck equations.In Appendix D we show the simplest such connection, using the continuum limit of a 1D lattice to recover the Fokker-Planck CD driving theory of Refs.[50,76].However there are still open questions, like what a continuum version of local control would look like.
The fact that there can exist many CD protocols for the same target trajectory, with distinct thermodynamic properties, means that one can search among these protocols for those that optimize certain quantities-like minimizing dissipated work under given physical constraints.Optimal control of nonequilibrium and finitetime processes is an active research area [77][78][79][80], with connections to techniques like Monge-Kantorovich transport theory [81] and trajectory-observable biasing within the framework of large deviations [82,83].Situating our driving theory within the broader context of these earlier optimal control approaches is an interesting topic for further study, both generically and in specific biological implementations in areas like ecology and evolution [84,85].
Driving a system between long-lived states is also subject to universal bounds or "speed limits" [86][87][88] that constrain the speed of driving in terms of dissipated work.Does CD or non-CD global driving saturate these bounds in certain circumstances?If so, are there biological implications, for example cases where natural selection has pushed a control process close to the theoretical limit?
Finally, are there analogous bounds when we only specify local control targets?
In summary, stochastic processes and their biological realizations are an ideal laboratory for investigating nonequilibrium control ideas.The driving framework we have developed is a particularly useful starting point, because the control protocols can be expressed analytically in terms of easy-to-calculate graph properties of the underlying Markov model.We can thus in principle explore a wide swath of control solutions, and identify generic features of control in diverse biological systems sharing similar graph topologies.The practicality of our formulation makes it well suited for deriving driving prescriptions in specific experimental contexts, like evolving cell populations [54] or the example systems in the current work.Because the protocols involve accessible control knobs-like varying drug/protein concentrationswe believe near-term experimental validation is within reach.Thus our approach may help with implementing control of biological systems in the lab, and also understanding how that control operates in nature.
The matrix can be decomposed as ∇ = ∇ − − ∇ + using Eqs.( 8)-( 9), where The reduced incidence matrix ∇ is given by the first two rows of Eq. (B1).Because ∆ = E − N + 1 = 1, we have ∆ + 1 = 2 trees in a tree basis.Taking the tree with edge 1 missing as the reference (tree 1 in Fig. 5B), we choose the other tree in the basis to be the one with edge 2 missing (tree 2 in Fig. 5C).Using the graphical algorithm, we can easily write down stretched inverse reduced incidence matrices for these trees: (B3) One can readily check that ∇[ ∇ (γ) ] −1 S for γ = 1, 2 is the 2×2 identity matrix.There is a single fundamental cycle vector for the graph, shown as a dashed line in Fig. 5C, given by c (1) = (1, 1, −1).
Using Eq. ( 32) we can write the currents J (t) to achieve a certain target ρ(λ t ) as where Φ 1 (t) is an arbitrary function.The backward rates k − (t) = k − = (k −r , k −c , k −x ) are fixed, and we can solve for the forward rates k + (t) = (k r R(t), k c C(t), k x X(t)) using Eq.(34).
The diagonal matrices M ± (t) = diag(∇ ± T ρ(λ t )) in Eq. (34) are given by Substituting the expressions from Eqs. (B3)-(B5) into Eq.( 34), we can solve for the concentration protocols R(t), C(t), X(t) that determine the forward rates: Different choices of Φ 1 (t) correspond to different control protocols that drive the system through the same trajectory ρ(λ t ).For example Φ 1 (t) = 0 gives the protocol associated with tree 1 (Fig. 5B), and Φ 1 (t) = ∂ t ρ 2 (λ t ) gives the protocol associated with tree 2 (Fig. 5C).The one additional constraint is that only Φ 1 (t) functions that lead to non-negative concentrations in Eq. (B6) at all t during driving are physically allowable.
To illustrate a family of control protocols, we need to choose a specific target trajectory ρ(λ t ).Based on the discussion in Sec.II A, Eq. (B6) describes the control protocol regardless of whether the target trajectory is an instantaneous stationary one or not.However to be concrete, we will choose an instantaneous stationary trajectory, making our control protocols CD.To mimic a rapid switch in gene expression from on to off, we will select ρ(λ t ) to be the stationary distribution associated with a set of concentration functions that serve as control parameters in the original system, λ t = (C(t), R(t), X(t)).We choose C(t) to sharply increase in a sigmoidal fashion, with R(t) kept at a constant level and X(t) in detailed balance with C(t) and R(t): where R 0 = 20 nM, C 0 = 0.2 µM, C f = 20 µM, k = 3 min −1 , t 0 = 5 min, and the remaining parameters are described above.Eq. (B7) determines the forward rates k + (t) that enter into the transition matrix Ω(λ t ), and hence allows us to use Eq.(3) to solve for the instantaneous stationary state ρ(λ t ).The components of this stationary target distribution are shown in Fig. 6A.They represent a transition from a system dominated by state 1 at the beginning of the protocol to one dominated by state 3 at the end (the gene turning mostly off).
The local detailed balance and optimal control protocols shown in Fig. 6 can be found numerically in a straightforward way.For the detailed balance case, we need to satisfy the condition χ cyc (t) = C T χ(t) = 0, where χ cyc (t) is a scalar because ∆ = 1.Plugging in the expressions from Eq. (B6), the χ cyc (t) = 0 condition becomes a nonlinear equation for Φ 1 (t).We can then use numerical root finding to determine the Φ 1 (t) at each t that satisfies the condition, thus defining the local detailed balance protocol.In a similar way, the optimal protocol (minimizing entropy production) should satisfy Eq. (38).When we plug Eq. (B6) into Eq.( 38), we get another nonlinear equation for Φ 1 (t), which can be numerically solved at each t.We verified that the resulting Φ 1 (t) is exactly the same as what we would get by direct numerical minimization of Ṡtot (t) in Eq. (35).We also checked that the solution for Φ 1 (t) always gives non-negative rates k + (t), or equivalently non-negative concentrations R(t), C(t), and X(t).
where m > 0 is the free energy difference between the intermediate and misfolded states.The rates k −n and k −u are related to their counterparts k n and k u through Here n and u are the free energy differences between the intermediate and native, and between the native and misfolded states respectively.Since going from states 1 → 4 → 3 should yield the same cumulative free energy difference as going directly from 1 → 3, we know that Since a full traversal of the left loop clockwise (states 3 → 1 → 2 → 3) involves hydrolysis of an ATP molecule, the product of clockwise/counterclockwise rate ratios over the entire cycle is related to the chemical potential difference ∆µ of ATP hydrolysis: We base the parameter values in our model on those associated with the chaperone GroEL assisting the folding of the substrate protein MDH, estimated from fitting to experimental data [68]: In cases where only upper or lower bounds on the parameters could be determined, we used the values at the bound.Using Eq. (C2) and the values of k u and k −u yield an estimate of u = 1.17 k B T .We do not know the precise value of m from the experimental fitting, but we assume a typical value of m = 3 k B T , which then gives k a , and hence we can neglect the effect of the k −a transition.

One-state local control
Let us first consider the one-state local control solution corresponding to the target subgraph of Fig. 8B.Here the number of target states N T = 1, and hence we want p 1 (t) = ρ 1 (t) for some chosen target trajectory ρ 1 (t), while the non-target states are given by π(t) = (p 2 (t), p 3 (t), p 4 (t)).The number of controllable edges is E c = 1, with only edge 1 in Fig. 8B amenable to external control (assuming fixed k a and negligible k −a ).The goal is to solve for the controllable edge rate vector k +c (t) = (k c C(t)) from Eq. ( 51), using the method outlined in Sec.II B. The various quantities needed to construct the solution are as follows.The submatrices of ∇ and G(t) are given by where we have set k −a ≈ 0 in G nπ because its effect is negligible.Because the target subgraph associated with ∇ ρc is tree-like, the stretched inverse [∇ ρc ] −1 S is just the ordinary inverse: [∇ ρc ] −1 S = (−1).Given these submatrices, we can calculate the matrix B G nπ and vector a(t) needed to evaluate the expression for π(t) in Eq. ( 48): Note that there is no dependence on arbitrary functions Φ c (t) in a(t), since ∇ ρc has no null space (∆ c = E c − N T = 0, and hence C c does not exist).The integral in Eq. ( 48) can then be carried out numerically to find π(t).Knowing π(t) allows us to evaluate the currents at the controllable edges J c (t) from Eq. ( 50).Finally we plug these currents into Eq.( 51) to find k +c (t) and hence our control protocol C(t).The matrices that appear in this equation are M +c (t) = (ρ 1 (t)), M −c (t) = (π 1 (t)), and the reverse rate vector is k −c = (k −c ), which cannot be externally varied.
To get the control results shown in Fig. 8C,D we chose a sigmoidal decreasing target function for ρ 1 (t) of the form neighboring lattice points.If we imagine the states as actual positions in a d-dimensional space, and allow the lattice spacing to become infinitesimal as the number of states N → ∞, then the behavior of such models should approach continuum diffusive dynamics described by Fokker-Planck equations.Thus, taking appropriate limits, we should be able to use our formalism to derive control solutions for Fokker-Planck systems.Here we describe how to do this for a d = 1 lattice in the global CD control case, rederiving the Fokker-Planck CD driving results of Refs.[50,76].Beyond this validation, we demonstrate how CD driving works for systems exhibiting position-dependent diffusivity, not considered in Refs.[50,76].
To connect our formalism to Fokker-Planck dynamics, let us first describe a one-dimensional Fokker-Planck equation for the time evolution of a probability density p(x, t), where x is our position variable, A(x, t) is the drift function, and D(x) is the position-dependent local diffusivity.Though D(x) is often taken to be a constant, D(x) = D, we allow it to be position-dependent for generality.We focus on the case where the drift A(x, t) = −D(x)∂ x E(x, λ t ), and hence arises from forces due to a potential energy E(x, λ t ) that may be dependent on time-varying control parameters λ t .Eq. (D1) can then be rewritten as where From the structure of Eq. (D2) it is clear that ρ(x, λ t ) is the instantaneous stationary distribution that makes the right-hand side vanish.We assume the energy function E(x, λ t ) → ∞ as x → x L and x → x R , defining a domain of x of width ∆x = x R −x L .Thus the partition function Z(λ t ) = x R x L dx exp(−βE(x, λ t )) is well-defined.An infinite domain would correspond to the special case where ∆x → ∞.The second line of Eq. (D2) defines a probability current density J(x, t), in terms of which the Fokker-Planck equation takes the form of a continuity equation for probability.
To apply our global CD control approach for discrete Markov systems, let us construct a one-dimensional lattice graph Markov model with N states, shown in Fig. 9A, that approximates the Fokker-Planck equation as N → ∞.State i corresponds to position x i = x L + ia, where a = ∆x/N is the lattice spacing, which becomes infinitesimal for large N .In this limit the probability p i (t) of being in state i is related to the probability density p(x, t) through a −1 p i (t) → p(x i , t).

(D4)
Here D i ≡ D(x i ) and E i (λ t ) ≡ E(x i , λ t ) are the discrete versions of the local diffusivity and potential energy.The ratio of the forward and backward transitions satisfies the local detailed balance relationship = e −β(Ei+1(λt)−Ei(λt)) . (D5) As a result, the instantaneous stationary distribution for this system assumes a form analogous to Eq. (D3), where Z(λ t ) = N i=1 exp(−βE i (λ t )).To check whether the transition rates of Eq. (D4) give the correct Fokker-Planck equation in the continuum limit, we note that the master equation for the discrete system can be written as: ∂ t p i (t) = j Ω ij (λ t )p j (t) = −J i+1 (t) + J i (t), (D7) where the current from state i to i + 1 is given by Eq. (D7) is the discrete analogue of the second line in Eq. (D2), with the conversion J i (t) → J(x i , t), a −1 p i (t) → p(x i , t).Plugging Eq. (D4) into Eq.(D8), we can rewrite the current J i (t) as (D9) Eq. (D9) goes to the correct limit in the continuum case, becoming the current density in the square brackets in Eq. (D2).To see this, note that a −1 ρ i+1 (λ t )ρ i (λ t ) → [ρ(x i + a, λ t )ρ(x i , λ t )] 1/2 ≈ [(ρ(x i , λ t ) + a∂ x ρ(x i , λ t )) ρ(x i , λ t )]   shown in the SI) but any valid discretization should lead to the same CD results in the continuum limit.
With the discretization validated, we can now proceed to applying the general solution procedure.The oriented current graph (N states, E = N − 1 edges) is tree-like, so the graph itself is the only spanning tree.Using the graphical algorithm we can write down the (N − 1) × (N − 1) dimensional stretched inverse reduced incidence matrix for this tree, [ ∇ (1) Because the graph is tree-like, the stretched inverse is also the ordinary inverse of the reduced incidence matrix, [ ∇ (1) ] −1 S = ∇ −1 .From Eqs. (D4) and (D6) we can deduce that the stationary currents have zero magnitude: Hence we know that J (t) = δJ (t).Moreover, since there are no cycles in the graph, Eq. ( 23) gives us the full CD current solution: Let us assume CD rates k + i (t) and k − i (t) of a form analogous to Eq. (D4), where D i (t) represents a modified, potentially timedependent, local diffusivity which we allow for generality, and E i (t) is the energy associated with state i in the CD protocol.In many cases it may not be possible to control the local diffusivity via external parameters, and hence it remains unchanged, D i (t) = D i .However as will be seen from the structure of the CD solution described below, we have in principle the freedom to choose D i (t) to be any non-negative function.The energy perturbation at each site due to the CD protocol is U i (t) = E i (t)−E i (λ t ).
To solve for these CD perturbations U i (t), the first step is to rewrite Eq. (D13) using Eq.(D11) and the expression for J (t) in terms of the CD transition rates: ∂ t ρ j (λ t ).(D15) After plugging in Eq. (D14) for the CD rates, and Eq.(D6) for the stationary distribution, Eq. (D15) can be written as: We can invert this to find a recursion relation for the .
(D17) Given an arbitrary choice of function U 1 (t) (which corresponds to the freedom of redefining the zero level for energies), we can use consecutive applications of Eq. (D17) to solve for U i (t), i = 2, . . ., N .
The final step is to transform the CD results back to the continuum, where the CD energies can be expressed as E(x, t) = E(x, λ t )+U (x, t).The perturbations U (x, t) can be found from the continuum analogue of Eq. (D17), ∂U (x, t) ∂x = 1 β D(x, t)ρ(x, t) To derive this we have expanded in small a and used the fact that sinh −1 ( ) ≈ to lowest order in .In the continuum limit i j=1 a ∂ t ρ j (λ t ) → x x L dx ∂ t ρ(x , λ t ) and a −1 ρ i (λ t ) ρ i+1 (λ t ) → ρ(x, λ t ), to leading order.This follows from the same argument as Eq.(D10), setting x = x i .
From Eq. (D18) we can directly solve for U (x, t), U (x, t) = U 0 (t) where x 0 is an arbitrary reference position and U 0 (t) is an arbitrary energy offset function (which does not affect the driving).
In practice, a particular CD protocol means simultaneously implementing the diffusivity D(x, t) and perturbing the energy landscape by U (x, t).As mentioned earlier, in many experimental scenarios control of diffusivity will not be possible, so the only available CD protocols will involve keeping the diffusivity equal to the value in the original system, D(x, t) = D(x).One special case of this is a position-independent diffusivity D(x) = D that is not varied during the CD protocol.This was solved by Li et al. [76] and Patra & Jarzynski [50] using alternative approaches, and their expressions for the CD perturbation are equivalent to our Eq.(D18) with the substitution D(x, t) = D.
From the perspective of thermodynamic costs, Eq. ( 35) for our discrete-state system takes the form where we have used the CD rates from Eq. (D14) and J i (t) = − i j=1 ∂ t ρ j (λ t ) from Eqs. (D11)-(D13).The functional form for Ṡtot (t) is always non-negative, since y sinh −1 (cy) ≥ 0 for any y when c ≥ 0. In the limit of adiabatically slow driving, ∂ t ρ j (λ t ) → 0, we see that J i (t) → 0 and hence the entropy production rate Ṡtot (t) → 0. As noted in Sec.I F, under the (unlikely) scenario that one can control the local diffusivity D i (t) and make it large during the CD protocol, then Ṡtot (t) can be made small even for fast driving.
In the continuum limit, Eq. (D20) becomes where J i (t) → J (x, t) is the continuum CD current.

5 FIG. 3 .
FIG.3.A two-loop discrete state Markov model, with N = 4 states and E = 5 edges.A) The black arrows correspond to entries in the transition matrix Ω(λt): transition rates k ± i (λt) that depend on an external protocol λt.B) The red arrows labeled α correspond to the oriented stationary currents Jα(λt), defined in Eq. (13).C) On the left, one of the spanning trees of the graph, chosen to be a reference for constructing the tree basis.Edges deleted to form the tree are shown in faint red.On the right, two trees in this set derived from the reference one.Each such derived tree has a one-to-one correspondence with a fundamental cycle of the graph (highlighted in green).
FIG.4.A) Partition of the graph matrices ∇ and G(t) into submatrices to facilitate solving the local control problem.B) A network with NT = 4 target states and EC = 6 controllable edges, depicted as indicated by the legend.The three target subgraphs, each containing at least one target state and all other states connected to it via controllable edges, are outlined in dashed curves.Because each target subgraph includes at least one non-target state, the local control condition is satisfied.C) Same network as B, except that the set of NT = 4 target states is different.Here neither of the two target subgraphs contain any non-target states, and hence local control is impossible.

1 FIG. 5 .
FIG.5.A) Biochemical network of a repressor-corepressor model, showing an operator site on DNA in three different states: 1) free; 2) bound to a bare repressor protein; 3) bound to a represssor-corepressor complex.Transition rates between the states are shown in green.The binding reaction rates depend on three concentrations of molecules in solution: R(t) for bare repressors, C(t) for corepressors, and X(t) for the complexes.B) One of the spanning trees for the associated network graph, with the edge deleted to form the tree shown in faint red.We take this to be the reference spanning tree for the tree basis.C) The other tree in the basis, with the corresponding fundamental cycle in green.

FIG. 6 .
FIG.6.A) Components of the target stationary distribution trajectory ρ(λt) (solid curves) for the repressor-corepressor system.B-D) Characteristics of four different control protocols that all drive the system along the target trajectory ρ(λt).The four protocols are: spanning tree 1 (violet, corresponding to Fig.5B); spanning tree 2 (teal, corresponding to Fig.5C); the solution satisfying local detailed balance, ∆µ(t) = 0 at all t (yellow); and the optimal solution that minimizes Ṡtot (t) at all t (thick black).For each case we depict: B) the total entropy production rate Ṡtot (t), with the inset showing the difference δ Ṡtot (t) ≡ Ṡtot (t) − Ṡtot,opt (t) between non-optimal and optimal rates on a log scale [units: kB/min]; C) the instantaneous chemical potential ∆µ(t); D) the concentrations C(t), R(t), X(t).

FIG. 8 .
FIG.8.A) Conformational states of a protein interacting with a chaperone.Transition rates in our kinetic network model are indicated by solid green arrows.Related transitions outside the scope of the model are shown as dashed arrows.B) If ka is fixed and k−aC(t) is negligible, there is effectively only one controllable edge in the network, the one between state 1 and 2. We can thus choose a target subgraph to enable one-state local control, with state 1 (misfolded) as our target.C-D) Example of a one-state local control solution, with a chaperone concentration protocol C(t) in panel C forcing the state 1 probability p1(t) to follow a target sigmoidal decrease in the misfolded probability ρ1(t) in panel D. The state probabilities pi(t), i = 1, . . ., 4 are calculated by numerical solution of the master equation, with parameters described in Appendix C. E-G) Same as panels B-D except with an additional controllable edge due to varying ATP levels, allowing ka(t) to be time-dependent.We can now have two-state local control, targeting states 1 and 2, and find a protocol of C(t) and ka(t) to make the system follow a chosen set of targets ρ1(t) and ρ2(t).

Appendix C: Details of the chaperone model calculations 1 .
Transition rates and model parametersThe transition rates of the chaperone model in Fig.8Asatisfy certain local detailed balance relationships.The rates k m and k −m , describing interconversion between states 3 and 1, obey

n
= m − u = 1.83 k B T .Similarly, we set k −c = 0.1 min −1 as the unbinding rate of the chaperone, a typical scale assuming strong binding affinity between the chaperone and substrate.The remaining unknown parameters can now be determined using Eqs.(C1)-(C3) (setting the ATP hydrolysis potential difference ∆µ = 22 k B T [89]): k −m = 0.0184 min −1 , k −n = 0.0585 min −1 , k −a = 0.381 M −1 min −1 .As mentioned in the Sec.III B, given typical chaperone concentrations C(t) ∼ O(1 µM) or smaller, we get k −a C(t)

1 / 2 =
ρ(x i , λ t ) + O(a), (D10)where O(a) denotes corrections of order a.Thus Eq. (D4) is a valid discretization of the Fokker-Planck system.It is not unique (other choices are possible, as