Riemannian Geometry of Optimal Driving and Thermodynamic Length and its Application to Chemical Reaction Networks

It is known that the trajectory of an endoreversibly driven system with minimal dissipation is a geodesic on the equilibrium state space. Thereby, the state space is equipped with the Riemannian metric given by the Hessian of the free energy function, known as Fisher information metric. However, the derivations given until now require both the system and the driving reservoir to be in local equilibrium. In the present work, we rederive the framework for chemical reaction networks and thereby enhance its scope of applicability to the nonequilibrium situation. Moreover, because our results are derived without restrictive assumptions, we are able to discuss phenomena that could not been seen previously. We introduce a suitable weighted Fisher information metric on the space of chemical concentrations and show that it characterizes the dissipation caused by diffusive driving, with arbitrary diffusion rate constants. This allows us to consider driving far from equilibrium. As the main result, we show that the isometric embedding of a steady state manifold into the concentration space yields a lower bound for the dissipation when the system is driven along the manifold. We give an analytic expression for this bound and for the corresponding geodesic, and thereby are able to dissect the contributions from the driving kinetics and from thermodynamics. Finally, we discuss in detail the application to quasi-thermostatic steady states.


I. INTRODUCTION
In general thermodynamics, the importance of the Hessian of the free energy function as a Riemannian metric has been recognized as early as 1975 by Weinhold [1][2][3] and by Ruppeiner [4]. In 1983 it was discovered by Salamon and Berry that this metric controls the dissipation when the system is driven endoreversibly between equilibrium states [5]. For this reason, the concept to determine curves of minimal dissipation from the Riemannian geometry of the parameter space was termed thermodynamic length. The derivation required that the driving is given by the explicit form dη i /dt = k(η i e − η) for the extensive variables η i characterizing the equilibrium state of the system and η i e being the respective variables of the reservoir, with a universal rate constant k. More recently, a derivation using a statistical approach was given by Crooks [6,7] based on a cycling driving protocol. In 2012, Zulkowski, Sivak, Crooks and DeWeese have presented an exact and explicit computation of geodesics for a driven particle diffusing in a one-dimensional harmonic potential [8].
These previous approaches were valid for endoreversible driving, i.e. for such processes that both the system and the driving reservoir are in local chemical equilibrium and dissipation takes place only at the boundary between the system and the reservoir. This theory does not directly apply to systems in nonequilibrium or those with algebraic constraints. Chemical reaction networks (CRN) are an important class of such systems and we extend the concept of thermodynamic to chemical reaction networks in equilibrium and nonequilibrium steady states.
The importance of chemical reaction networks derives from the fact that they form the basis to model and understand complex chemical processes on a molecular level [9,10]. They are central to the understanding of biological function as well as in many industrial applications [11][12][13][14][15][16], where tight control and high efficiency of the reactions are mandatory.
For CRN, a physically meaningful driving dynamics should be given by dx i /dt = k i (x i e − x), where x i is the concentration the i-th chemical, x i e the concentration of the respective chemical in the reservoir and k i the diffusion rate constant across the boundary. This dynamics generalizes the one given by Salamon and Berry and, more importantly, it can drive the CRN away from chemical equilibrium or from a steady state. Furthermore, models of CRN found in biochemical applications are often open CRN, and as such might not exhibit equilibrium states. For such networks, the dissipation caused by the driving between steady states is of interest, yet it cannot be determined within the framework available thus far.
The mathematical foundation of CRN theory was laid by Horn, Jackson and Feinberg [17][18][19][20][21], and the theory has been refined and developed by employing various mathematical disciplines such as graph theory [22][23][24], homological algebra [25], and differential geometry [26]. In recent years, it has become clear that CRN are most naturally analyzed using geometrical methods: Algebraic geometry has led to various important foundational results [27][28][29][30] for CRN with mass action kinetics as well as applications to biologically important networks and phenomena [31][32][33]. Meanwhile, information geometry has been successfully used to analyze the geometry of thermodynamics within a stoichiometric compatibility class [34,35]. Finally, both approaches could be merged within the framework of Hessian geometry [36][37][38][39], generalizing the global structure of the concentration space as studied in algebraic geometry and equipping it with the thermodynamic structure from information geometry as well as including the more classical results [40][41][42][43][44][45][46][47] on the thermodynamics of CRN. In [34,35,[37][38][39], the information geometric aspects of Hessian geometry, based on the Bregman divergence, have been thoroughly explored and linked to thermodynamics. In this work, we focus on the differential geometric aspects of this Hessian geometry.
We show that the dissipation rate caused by the driving of a CRN is determined by the square of the norm j D 2 g K X of the driving flux j D with respect to the Riemannian metric diag(x 1 k −1 1 , . . . , x n k −1 n ) on the concentration space. This holds arbitrarily far from chemical equilibrium and even for CRN without equilibrium states. The Riemannian distance L X := d X (x 0 , x 1 ) between two concentration vectors x 0 and x 1 on the concentration space is called thermodynamic length. Its square is proportional to the minimal dissipation due to the driving in the case that the driving flux j D and the reaction flux j R obey the time scale separation j D j R , which allow to ignore the reaction effects. The geodesic equations can be solved explicitly to yield the optimal trajectory of the system x(t) and the optimal driving protocol x e (t). This is Result 1, which is derived in Sec. II. The total minimal dissipation Σ min X of this process can be explicitly computed as This is Result 2 of this work. Without any assumptions on the time scales of j D and j R , the minimal total dissipation Σ min V ss caused by the driving of the system along the steady state manifold V ss between x 0 anf x 1 on V ss is proportional to the square length of the Riemannian distance between x 0 and x 1 on V ss (Result 3). The optimal system trajectory is thus a geodesic and the optimal driving protocol can be determined from the respective geodesic equations. The main result (Result 4) follows from the isometric embedding of V ss into concentration space and states that Σ min V ss is bounded from below by Σ min X , i.e.
This lower bound is universal in the sense that it holds for all possible steady state manifolds irrespective of the reaction kinetics of the CRN. The Results 1-4 are derived in Sec. II. In Sec. III, the general results from Sec. II are applied to the important class of quasi-thermostatic CRN. By definition, these are all CRN whose steady state manifold is given by where x ss is any steady state and S is the stoichiometric matrix. This class includes all CRN with equilibrium states that obey the thermodynamics of an ideal solution and all complex balanced nonequilibrium steady states. An explicit parametrization akin to the exponential family in information geometry is used to derive an explicit expression for the Riemannian metric on the parameter space, which yields numerically solvable geodesic equations for any optimal driving problem. This parametrization is complemented by a parametrization via conserved quantities, and its physical importance is discussed. In particular, the parametrization via conserved quantities is an important tool to treat the problem of arbitrary driving protocols with the time scale separation j R j D , which we aim to investigate in the future. Finally, in Sec. IV, we discuss possible extensions of our approach to systems with non-diagonal Onsager relations.

Setup and notation
Notation. The logarithm and exponential functions of vectors are taken componentwise and the resulting vectors are elements of the linear dual of the original space, i.e.
where {e i } is a basis for the vector space and {e i } the respective dual basis. The reason for this convention is that the exp and log maps appearing here are Legendre transformations.
By the symbol • we denote the Hadamard product between vectors, i.e.
For the differential geometric formalism, we refer to the textbook [48].
Chemical reaction networks (CRN). In this work, we consider a chemical reaction network consisting of n chemicals X 1 , . . . , X n and r reactions R 1 , . . . , R r . The jth reaction is given by with nonnegative integer coefficients (S + ) i j and (S − ) i j . These coefficients determine the reactants and products of the reaction [49]. The structure of the network is thus encoded in the n × r stoichiometric matrix S = (S i j ) with matrix elements The state of the reaction network is characterized by a vector of nonnegative concentration values x = (x 1 , . . . , x n ) ∈ R n >0 , where x i represents the concentration of the chemical X i . We denote the concentration space by X := R n >0 . Finally, the dynamics of the CRN is governed by the equation where j R = (j 1 R , . . . , j r R ) is the vector of reaction fluxes and j D = (j 1 D , . . . , j n D ) is the vector of the external driving fluxes. Choosing a kinetic model is tantamount to specifying j R as a function of x. For any given kinetic model j R (x), we denote the respective steady state manifold by The results in this work have general validity as they are not based on the choice of a particular kinetic model for j R .
Chemical thermodynamics. In chemical thermodynamics, the state of the system is characterized by concentration vectors x ∈ X ⊂ R n or, equivalently, by chemical potential vectors µ, which live in the dual space Y := (R n ) * . We denote the bilinear pairing between the dual vector spaces by ., . . The concentration and potential vectors are connected by Legendre duality via the strongly convex free energy function ϕ(x) and its Legendre dual ϕ * (µ) as ϕ * (µ) = max x { x, µ − ϕ(x)} and µ(x) = argmax µ { x, µ − ϕ * (µ)} and analogous variational characterizations of ϕ(x) and x(µ). Additionally, the equality is satisfied for the Legendre dual pair of x and µ. Hereby, the potential ϕ(x) is the Gibbs free energy [50], which takes the form for an ideal dilute solution (or, equivalently, an ideal gas). The vector µ 0 ∈ (R n ) * is the vector of standard chemical potentials and we choose the energy scale such that k B T = 1. This explicit form of ϕ(x) will be used in calculations throughout the text. This yields the xdependence of µ as and, using Eq. (2), the explicit form of the Legendre dual potential function [51] as The strictly convex function ϕ(x) on the concentration space X gives rise to a Riemannian metric via its Hessian This is known as the Fisher information metric in information geometry [52] and Weinhold or Ruppeiner metric in thermodynamics [1,4]. Because the spaces we work with have global coordinate systems, the metric can be globally represented by a matrix. For the specific form of the convex function ϕ(x) given in Eq. (3), the metric g X is represented by the diagonal matrix To account for the driving kinetics in the next section, we introduce a diagonal weight matrix K = diag(k 1 , . . . , k n ) with k i ∈ R >0 and the weighted Fisher information metric on X as which relates to the Hessian metric g X as g K The Legendre duality between X and Y , together with the Hessian metric g X and the Bregman divergence is the setup leading up the Hessian geometry of CRN, which was thoroughly explored in [37] and [38]. Here, we focus on the Riemannian geometry based on the metric g K X .

II. DISSIPATION VIA THERMODYNAMIC LENGTH
Throughout this section, we analyze the case of an externally driven CRN between two states x 0 , x 1 ∈ X in finite time T . The full dynamics of the driven CRN is given by Eq. (1). We let x(t) t∈[0,T ] denote the resulting integral curve and also have x 0 = x(0) and x 1 = x(T ). The driving is due to the flux j D , which is controlled by the concentrations of chemicals in the reservoir, and the resulting dissipation is evaluated in Sec. II A. If the driving is faster than the chemical reactions, i.e. j R j D , then the dynamics reduces to and the optimal trajectory x(t) t∈[0,T ] is a geodesic in the concentration space X, endowed with the weighted Fisher information metric g K X . In this case, the optimal trajectory and the minimal dissipation can be calculated analytically, as shown in Sec. II B.
If such a time scale separation does not exist, then the problem of treating an arbitrary driving field j D is rather intricate as the contributions to the dissipation and to the dynamics generated by j D and by j R mix.
Even if the time scale separation j R j D holds, although this constrains the trajectory x(t) t∈[0,T ] to the steady state manifold, the contribution of the dissipation due to chemical reactions cannot be neglected. Because of the dependence on the details of the particular CRN and of the driving protocol, we do not treat this general case here.
A possibility to circumvent the mixed contributions of j D and j R consists in restricting the driving field j D to the tangent space T x(t) V ss of the steady state manifold and to impose the initial condition x 0 ∈ V ss . This implies that the dynamics is again described by Eq. (8) and that the trajectory x(t) t∈[0,T ] is restricted to the steady state manifold. In this case, as shown in Sec. II C, the optimal trajectory is given by a geodesic on V ss and can be explicitly computed whenever a parameter manifold for V ss is known.

A. Dissipation
We assume that the system is in contact with an external reservoir and can exchange all its chemicals X 1 , . . . , X n via diffusion with the reservoir, with diffusion rate constants k 1 , . . . , k n . The reservoir is thereby characterized by the time-dependent concentration vector x e (t) = (x 1 e (t), . . . , x n e (t)) and the corresponding chemical potential vector µ e (t), which satisfy µ e (t) = µ 0 + log x e (t). The diffusion flux is given by Fick's law, i.e. the flux vector field on X is where K = diag(k 1 , . . . , k n ) is the matrix of diffusion rate coefficients. In the following, we suppress the timedependence of the variables for the sake of clarity of exposition. In addition, we assume that the external reservoir potential vector is close to the system potential vector, i.e. µ e ≈ µ. Then the relation holds [53], where g X is the Fisher information metric introduced in Eq. (6). This leads, using Eqs. (9) and (10), to the dissipation rate at the boundary as where g K X is the weighted Fisher information metric introduced in Eq. (7). This expression is the reason for introducing the weights to the standard Fisher information metric: The dissipation rate at the boundary between the system and the reservoir is given by the squared norm of the driving vector j D (x, t). Together with Eq. (8), this allows to relate the dissipation rate to the squared speed of the integral curve x(t) and the total dissipation to its energy (here, energy is used in the context of differential geometry, which has no relation to the thermodynamical energy).

B. Optimal fast driving
If the time scale separation j R j D holds, which we call the fast driving regime, the optimal system trajectory x(t), the corresponding driving protocol x e (t) and the minimal dissipation can be determined analytically. In the fast driving regime, the dynamics of the driven system obeys Eq. (8) and therefore the dissipation rate is given by Eq. (11): Applying the Cauchy-Schwarz inequality to the integral yields where L X is the length of the curve x(t). For the trivial weight matrix K = I, this is known as thermodynamic length [5] and we also use this term for the lengths of curves on spaces endowed with the weighted Fisher metric. The equality in Eq. (13) holds if and only if the speed of the curve, dx dt g K X , is constant. In this case, x(t) is a geodesic. If it is a minimal length geodesic, then both the length L X of the curve and its energy Σ X are minimized. In this case, L X is equal to the distance between x 0 and x 1 on the space X, i.e.
For the driving field given by diffusion, i.e.
, and the weighted Fisher information metric g K X , the geodesic equations are for i = 1, . . . , n, cf. Appendix A. They have the explicit solution for t ∈ [0, T ]. The second term of the expression describes a linear trajectory in concentration space, which would yield the shortest path in the standard Euclidean metric, and the first term shows the correction due to the non-Euclidean geometry of the concentration space. Interestingly, the optimal curve x(t) does not depend on the details of the driving kinetics, i.e. it is identical for all possible values of the diffusion rate constants k i . In other words, the optimal curve is determined purely by the thermodynamics of the system via the Hessian g X . This is the first main result of this section: Result 1. The optimally driven curve x(t) in the fast driving limit is given by Eq. (14). It does not depend on the diffusion rate constants but only on the thermodynamics of the system. The optimal driving protocol x e (t) is given by which follows from Eqs. (8) and (9).
With the explicit solution at hand, Eq. (11) yields the dissipation rate and Eq. (12) the total dissipation in the driving process and thus the second main result: Result 2. The optimally driven system under fast driving has the constant dissipation rate and the minimal total dissipation is The results are not only significant in the fast driving regime but they provide lower bounds for dissipation when the trajectory x(t) is restricted to the steady state manifold, as is shown in the following section.

C. Driving between arbitrary steady states
In this section we consider the case that the driving field j D is constrained to the tangent space T x(t) V ss of the steady state manifold with the initial condition x 0 ∈ V ss . This implies that the integral curve x(t) is restricted to V ss and therefore the term Sj R (x(t)) in Eq. (1) vanishes for all points on the curve by definition of V ss . Thus, the dynamics is given by dx dt = j D and the dissipation rate is given by σ D = dx dt g K X as before. Now, the integral curve is restricted to V ss and by the same argument as before, the minimal length curve on V ss between x 0 and x 1 will minimize the dissipation. This can be made precise (and accessible to direct computation as demonstrated in Sec. III) as follows: Without loss of generality we assume that V ss is connected and thus parametrized by a connected manifold M , i.e. there is an embedding with Im[f ] = V ss . Defining the metric on M as the pullback makes the embedding f isometric and thus the dissipation rate can be written as where f (m) = x. The total entropy production is bounded from below by where L V ss is the length of the curve m(t) on M , which is bounded from below by the distance d V ss (m 0 , m 1 ) be- This leads to the following result: Result 3. In the case that the system is driven on the steady state manifold, the dissipation Σ V ss , which is caused by the driving, is minimized if and only if m(t) is a geodesic on M . In this case, the minimal dissipation reaches its lower bound, which is given by The integral curve on the concentration space is x(t) = f (m(t)) and the corresponding driving protocol is given explicitly by which follows from Eq. (15) and where J f is the Jacobian of f at m(t). This result holds independently of any time scale separation between j D and j R .
Finally, the isometry of the embedding f implies that and thus Σ min X ≤ Σ min V ss . This is illustrated in Fig. 1. We are led to the main result of this work: Result 4. The total dissipation Σ V ss for a driving process on the steady state manifold is bounded from below by This lower bound is achieved if and only of the system is driven on the optimal curve for the fast driving regime, given by Eq. (14).
x 0 1. Illustration of the approach to the lower bound for the total dissipation Σ min V ss . The isometric embedding of V ss (represented by the green surface) into the concentration space X allows to compare the lengths of the geodesics dVss and dX on V ss and X between the two points x0, x1 ∈ V ss as dVss ≥ dX . The relations Σ min The latter can be explicitly computed, cf. Eq. (16), and the former is a lower bound for the dissipation when the system is driven on V ss with any protocol, thus yielding Result 4. Note that nonequilibrium steady states have a nonzero dissipation rate caused by continuously ongoing chemical reactions. Our approach does not take this dissipation into account, but only gives a recipe for minimizing the dissipation caused through the driving.
For an arbitrary CRN, the parametrization of the steady state manifold f : M → X is, in general, not available. In fact, even when j R is given by mass action kinetics, the determination of steady states involves highly nontrivial algebraic geometry and must be carried out for each CRN individually [54]. Yet, the Result 2 holds irrespectively of any of the complicated details of V ss .
For the class of quasi-thermostatic CRN, an explicit parametrization of V ss is available and the Result 3 can be made precise and accessible for direct calculation. This is presented in the next section.

III. DRIVING IN QUASI-THERMOSTATIC CRN
Quasi-thermostatic CRN, as defined by Horn and Jackson in [17], are characterized by the following particular form of the steady state manifold where x ss is a particular solution of the equation Sj R (x ss ) = 0. The quasi-thermostatic property derives its importance from the fact that all equilibrium states of ideal solutions and complex balanced steady states of CRN with mass action kinetics are quasi-thermostatic. The class of quasi-thermostatic CRN is, however, even broader because there exist quasi-thermostatic CRN that belong to neither of the two classes mentioned above [55].
In the following, we give two possible parametrizations of quasi-thermostatic steady states and illustrate each with an explicit calculation of the geodesic and minimal dissipation based on Result 3.

A. Exponential Parametrization of Quasi-thermostatic CRN
By choosing a basis u * 1 , . . . , u * q of Ker[S T ] and writing U * = (u * 1 , . . . , u * q ) for the n × q matrix of basis vectors, the characterization in Eq. (17) can be rewritten as Thus, we obtain the parametrization of V ss via the embedding This structure has been thoroughly studied in algebraic geometry [27,30], where V ss was shown to be a toric variety and, using the matrix-tree theorem explicit formulae for x ss were obtained. In the particular case that V ss is the equilibrium manifold of a closed CRN, x ss is given by exp(−µ 0 ). In this work, the analytification V ss of the toric variety from algebraic geometry is used. We call V ss a toric manifold.
We now treat the Riemannian geometry of the quasithermostatic steady state manifold V ss . The pullback of the metric g X to R q can be computed as With the Jacobian explicitly given by J f = DU * , where D := diag(x 1 , . . . , x n ). Note that D can be explicitly expressed in the η * coordinates by using Eq. (19). We obtain the desired metric on the parameter space R q as Although we could not solve the geodesic equations on R q in an analytically closed form, they can be solved numerically for any given CRN. The following example illustrates this. Example. Consider the following nonlinear chemical reaction network A 4B, with stoichiometric matrix S = (−1 4) T . The space Ker[S T ] is spanned by the vector u * = (4 1), the space of η * variables is one-dimensional, i.e. it is given by R 1 , and the metric g K R 1 is given by according to Eq. (20). The geodesic equation for η * is ss . Fig. 2A shows the geodesic on the steady state manifold determined by Eq. (21) and the geodesic on the space X, which is determined by Eq. (14), with the numerical parameter values given in the figure caption. The curves are color coded by the speed of the driving, which is given by the Euclidean norms of dx dt and df (η * ) dt , respectively. The speed of the curves is also shown in Fig. 2B and as a function of time. The speed increases with increasing norm of the coordinate x for both curves, whereby the system driven on V ss has the greater acceleration. As both curves are geodesics, their speed with respect to the metric g K X must be constant, as is confirmed by the numerical results shown in Fig. 2(B). The total dissipation for the optimal driving along V ss 17. For linear CRN, i.e. CRN where all the reaction complexes consist of exactly one chemical, the geodesic on X, given by Eq. (14), and the geodesic on V ss coincide (see Appendix C). In this case the equality Σ min V ss = Σ min X holds and Σ min V ss can be computed by Eq. (16). In order to illustrate the difference between Σ min V ss and Σ min X we have thus chosen an example with a strong nonlinearity. But even in the nonlinear case, Σ min X provides a good estimate for Σ min V ss . This is not only the case for the given example but for various other CRN that we have analyzed numerically. In future work, we aim to analytically quantify the error in estimating Σ min V ss by Σ min X by relating it to the curvature of V ss .

B. Parametrization via conserved quantities
In addition to the parametrization of the steady state manifold V ss of quasi-thermostatic CRN as a toric manifold, given by Eq. (19), the manifold V ss can be with x0 = f (η * 0 ) and x1 = f (η * 1 ). A. Trajectories of the system being driven from x0 to x1. The lower curve x(t) is the optimal curve on the whole concentration space, given by Eq. (14) and the curve xe(t) (yellow) is the corresponding optimal driving protocol, given by Eq. (15). The curve f (η * (t)) is the optimal curve on the steady state manifold V ss (light blue). The curves x(t) and f (η * (t)) are color coded according to the speed of the driving, given by the Euclidean norm of the tangent vector dx dt . In B, the Euclidean norms, which give the speed of the driving, and the . g K X norms (their 10-fold values are plotted here for better visibility), which govern the dissipation, of the tangent vectors x(t) and f (η * (t)) are plotted for the course of the driving.
parametrized by the vector of conserved quantities η = (U * ) T x, with the matrix U * introduced in the previous section and x ∈ V ss [56]. This is based on Birch's theorem [27,57], which states that the intersection of V ss with the stoichiometric polytope P (η) := x ∈ X|(U * ) T x = η exists and is unique, i.e. the map (U * ) T : V ss → R q , defined by x → (U * ) T x, has a unique inverse given by where E := (U * ) T V ss is the image of V ss . This is the parametrization of quasi-thermostatic steady states by the vector of conserved quantities. The geometrical reason for the uniqueness of this intersection is the generalized orthogonality between V ss and P (η), as discussed in [37,38]. Fig. 3 illustrates this. 3. Illustration of the parametrization of points on a manifold of quasi-thermostatic steady states V ss via the vector of conserved quantities η. The codimension of V ss in X is the dimension of the stiochiometric polytope P (η) and therefore a zero dimensional intersection set is expected in general. The dual orthogonality from Hessian geometry [37,38] or, alternatively, Birch's theorem from algebraic geometry [27], ensure that the intersection consists of exactly one point for each η ∈ E. As such, this structure is a folitation of X, with the base manifold being V ss and the leaves being P (η). In other words, V ss is the moduli space of stiochiometric compatibility classes of the CRN. In the figure, it is illustrated how the P (η) and P (η ) are parallel affine spaces in X, that cover all of X and are indexed by points on V ss .
Although the map h cannot be given in analytically closed form, its Jacobian can be computed. First note that the map h −1 f is the coordinate change from η * to η variables and its Jacobian is given by where f is the map given in Eq. (19). Using J f = DU * leads to Thus, defining the weighted Fisher information metric g K E on E by makes the parametrization h : E → V ss isometric [58]. This is summarized in the following diagram With the respective metrics g K E , g K R q and g K X V ss , the maps in the upper row of the diagram are isometries and thus the dissipation through finite time driving can be computed on either of the three spaces. The embedding ι : V ss → X is isometric and allows for the explicit computation of lower bounds for the total dissipation due to the driving as stated in Result 4.
Finally, we remark that the coordinates η and η * are Legendre dual, see [38] for details.

C. Applicability of the two parametrizations
The metric g K R q does have an analytically closed expression in the η * -coordinates via Eqs. (19) and (20). Therefore, the space R q of η * -parameters, equipped with the metric g K R q , is suitable for explicit computations. The example in Sec. III A can be generalized to arbitrary quasi-thermostatic CRN, i.e. the geodesic equations in η * -coordinates on R q can be explicitly written and numerically solved for a given set of parameters to obtain the optimal system trajectory, the optimal driving protocol and the minimal total dissipation.
In contrast, there is no analytically closed expression for the metric g K E on the space E of η-coordinates. However, the vector of conserved quantities η has an intuitive physical meaning and is a natural tool to analyze the case of the time scale separation j R j D . In the case of arbitrary driving fields j D , the time scale separation j R j D results in a slow time scale dynamics on V ss , which is determined by the pushforward (U * ) T x(t) = η(t) of the trajectory x(t) to the space of conserved quantities. While the detailed analysis of this scenario is dependent on the details of the CRN, the slow dynamics on V ss can be reproduced by the trajectory h((U * ) T x(t)). The corresponding driving protocol is given by However, this formula requires the evaluation of the function h. This can be circumvented by transforming the original driving field j D into a driving field j * D on the tangent space T V ss of the steady state manifold. For each x ∈ V ss and t ∈ [0, T ], the change in steady state, which is the state that the system will eventually relax to after the disturbance caused by j D (x, t), is described by the change of the vector of conserved quantities. The latter quantity is given by the pushforward J (U * ) T j D (x, t) to the tangent space T (U * ) T x E. There is one and only one vector field j * D (x, t) on T x V ss , which causes the same change in the vector of conserved quantities as j D (x, t). It must be given by the pushforward of Using the explicit formula for the Jacobian J h given in Eq. (23), one obtains the desired vector field as This approach is illustrated in Fig. 4. One verifies by a direct calculation that J (U * ) T j * D = J (U * ) T j D and therefore j D and j * D indeed induce the same slow dynamics on V ss . Moreover, in the case that j D is already a vector field on the tangent space of V ss , the equality j * D = j D holds. Therefore, when considering arbitrary driving fields and the time scale separation j R j D , this construction plays an important role. In future work, it has to be supplemented by an analysis of the dissipation caused by the relaxation of the system to V ss caused by j R .

IV. DISCUSSION AND OUTLOOK
In this work, we have analyzed the application of the Riemannian geometry with a weighted Fisher information metric to the dissipation in driven chemical reaction networks. This geometry should be thought of as an infinitesimal and weighted version of the Hessian geometry established in [37,38], where information geometric aspects based on the Bregman divergence were analyzed. Inspired by the classical results on thermodynamic length, we have shown that the dissipation rate is given by the speed of the integral curve if either the diffusion is fast or the driving generates curves that lie on the steady state manifold. Thereby, we have presented a mathematically and physically precise derivation of this formula based on the physical model of diffusion of chemicals between the system and the reservoir. Moreover, we were able to explicitly include the diffusion rate constants into the calculation. This was the reason to introduce weights in the various Fisher information metrics.
The careful mathematical setup allowed us, almost effortlessly, to obtain several new results for the dissipation in driven chemical reaction networks. If the diffusion of 4. The construction of a driving vector field j * D tangent to V ss from an arbitrary driving field jD such that both fields generate the same slow dynamics on V ss . At each point xss ∈ V ss , the vector field jD will drive the system away from V ss in general. However, if the time scale separation jR jD holds, the steady state reached after such a perturbation will be determined by the change in the vector of conserved quantities ∆η = η − η, i.e. by the pushforward J (U * ) T jD to the tangent space T (U * ) T xss E. All tangent vectors at xss, which lie in P (η ) have the same pushforward. Among them, there is exactly one, which is tangent to V ss . This is the desired driving field j * D . Note that in this construction, both stoichiometric polytopes P (η) and P (η ) lie in Tx ss X and not in X.
the driving is faster than the chemical reactions, a complete analytical solution was given: The optimal curve in concentration space, the optimal driving protocol and the minimal total dissipation were determined. Interestingly, it was found that the optimal curve does not depend on the diffusion rate constants, i.e. it is determined purely by the thermodynamics of the system, whereas the optimal driving protocol depends on the kinetics.
Then, without any assumption on the time scales, if the driving proceeds along the steady state manifold, the geodesic equations were formulated on the respective parameter space. For the case of quasi-thermostatic steady states, the equations were derived explicitly and thus they can be numerically solved for any given CRN. More importantly, the isometric embedding of the steady state manifold into the concentration space was used to obtain a lower bound for the total dissipation. This bound is universal as it is independent of the particular kinetic model j R (x) and is valid for any steady state manifold. In all numerical examples we investigated thus far, this bound turned out to provide a good estimate for the actual total dissipation. Thereby the estimation error seems to grow with the curvature of the steady state manifold and it will be a future challenge to make this precise.
Seen through the lens of general thermodynamics, this work shows that the concept of thermodynamic length can be extended to the nonequilibrium situation. This is necessitated by the physically meaningful driving of the concentration variables x instead of the extensive parameters η, which characterize the equilibrium system. In regard to the original derivation in [5], where only a specific driving field on the space of η parameters is considered, we have shown how the isometric embedding into a larger state space can yield a more detailed understanding of the driving kinetics and enable the calculation of explicit bounds.
The explicit analytical results obtained in this work shed new light on time-reversed driving: Results 1 and 3 show that the optimal system trajectory for the driving from x 1 to x 0 will be the time-reversed trajectory caused by the optimal driving from x 0 to x 1 . However, the driving protocols will follows different paths (see Appendix C for a discussion of time reversal illustrated with an example). This leads to the conclusion that the approach taken in [6] to compute the total dissipation for the driving between x 0 and x 1 is strictly valid only in the T → ∞ limit for the trivial weight matrix K = kI with a uniform rate constant k ∈ R >0 .
In future work, it will be interesting to address the general case of an arbitrary driving protocol and without any assumptions on time scale separation or restrictions of the direction of the driving field. Thereby, the results from [37], in particular the Pythagorean theorems, will be useful to evaluate the contribution to dissipation through chemical reactions.
However, there are more conceptual questions that have been raised by the presented analysis. The weighted metric g K X should be thought of as an Onsager matrix for the fluxes j D and the corresponding forces µ e −µ. This is in line with the interpretation of the metric as the generalized friction tensor given in [8]. It will be a rewarding future challenge to work out the geometry of optimal driving for arbitrary Onsager matrices. Moreover, the fact that the calculated optimal curve did not depend on the particular details of the driving kinetics and was of thermodynamical origin requires a thorough physical explanation. It might be necessary to employ a mathematically more advanced framework to disentangle the contributions from kinetics and thermodynamics, which are mixed within the matrix g K X in our current approach. Apart from the theoretical insights and questions raised in this work, we have provided another tool to better understand biochemical reaction networks and to aid in the design and optimization of operational protocols in industrial chemical applications.
where η ss = u * , x ss = n i=1 x i ss . Thus, the metric on E = R >0 is given by Eq. (24) as where η 0 = η(0) and η 1 = η(T ). The ideal driving on V ss between the points x 0 = η0 ηss x ss and x 1 = η1 ηss x ss thus follows the trajectory The total dissipation is given by This dissipation Σ min V ss is equal to the lower bound Σ min X stated in Result 2 and thus the driving on the equilibrium manifold coincides with the optimal driving for the fast driving regime given by Eq. (14). This can also be verified by directly substituting Eq. (B1) into Eq. (B2).

Appendix C: Reversal of driving
The Results 1 and 3 imply that a system which is optimally driven from x 0 to x 1 will follow the same trajectory as the optimally driven system from x 1 to x 0 with time-reversal: Denoting the former trajectory by x f (t) and the latter by x b (t), the relation x b (t) = x f (T − t) for all t ∈ [0, T ] holds. However, the same results show that this is not true for the driving protocols x f e (t) and x b e (t). They lie on different curves and cannot be identified under time reversal, even if one would introduce a lag time. This is illustrated in Fig. 5. Therein, the example from Sec. III A is considered with the additional assumption that there are no chemical reactions between A and B and thus that all of X is the equilibrium manifold. The driving protocols x f e (t) and x b e (t) are different non-intersecting curves.
FIG. 5. Reversal of driving: an example. The example from Fig. 2 with the assumption that no reactions take place between A and B is used here. The optimal protocol x f e (t) for the driving from x0 to x1 is shown in yellow and the optimal protocol x b e (t) for the driving from x1 to x0 is shown in blue.
For x f e (t) and x b e (t) to become the same curves, it is necessary and sufficient that the matrix K is proportional to the identity matrix and that the time limit T → ∞ is taken. In this case, they are identified with the geodesic x(t).