Non-steady state mass action dynamics without rate constants: dynamics of coupled reactions using chemical potentials

Comprehensive and predictive simulation of coupled reaction networks has long been a goal of biology and other fields. Currently, metabolic network models that utilize enzyme mass action kinetics have predictive power but are limited in scope and application by the fact that the determination of enzyme rate constants is laborious and low throughput. We present a statistical thermodynamic formulation of the law of mass action for coupled reactions at both steady states and non-stationary states. The formulation uses chemical potentials instead of rate constants. When used to model deterministic systems, the method corresponds to a rescaling of the time dependent reactions in such a way that steady states can be reached on the same time scale but with significantly fewer computational steps. The relationships between reaction affinities, free energy changes and generalized detailed balance are central to the discussion. The significance for applications in systems biology are discussed as is the concept and assumption of maximum entropy production rate as a biological principle that links thermodynamics to natural selection.


Introduction
One hundred and fifty years ago Peter Waage and Cato Maximillian Guldberg published their first article describing the law of mass action, that the rate of a chemical reaction is proportional to the concentration of the reacting species [1][2][3][4]. For a reversible reaction, labeled by +1 for the forward reaction and −1 for the reverse reaction, in which γ A moles of chemical A are transformed into γ B moles of chemical B, the forward rate due to reaction 1 is simply, The brackets indicate the concentration, and the constant of proportionality k 1 is known as the rate constant. A similar relation exists for the reverse reaction −1. The net rate is given by, All introductory chemistry texts describe the law of mass action in one form or another. Although the relationship is simple and can easily be applied to many reactions, the application to more complex systems such as biological metabolism is challenging because most rate constants are not available and measuring the missing rate constants is very labor intensive. Typically rate constants are measured by first isolating and purifying the individual enzymes and then determining the rate constants spectrophotometrically in vitro. The process is challenging to implement in a high-throughput manner and in principle, measurements would have to be done for each enzyme in the system of interest because even enzymes having high sequence similarity to the orthologous enzyme in another species can still have widely differing rate parameters. For example, for the enzyme dihydrofolate reductase, the turnover rates for the substrate 7,8-dihydrofolate measured in vitro vary five orders of magnitude across species-from 284 s −1 to less than 1 s −1 . Consequently, thousands of rate parameters would be needed for each new genome, and even then the traditional kinetic analysis provides measurements of rates in vitro rather than in vivo. Thus, while a large number of enzyme rate parameters have been determined for a few well-studied model systems such as Escherichia coli and Saccharomyces cerevisiae, this is Non-steady state mass action dynamics without rate constants: dynamics of coupled reactions using chemical potentials not true for most systems of interest. As a result of the difficulty in obtaining rate constants, constraint-based flux models have been the method of choice for largescale modeling of biological processes such as metabolism. Unfortunately, the lack of physics in regard to the law of mass action limits the ability of these approaches to narrow the solution space to only physically likely solutions. A considerable amount of work has focused on extending constraint-based models to include the physics of mass action in terms of thermodynamics constraints [5], hybrid models [6] or dual variables [7].
Thermodynamic [8][9][10][11][12] and other approaches [13][14][15][16][17] have been proposed for the study of elementary reactions using the law of mass action that do not use rate constants, and instead have used the concept of reaction affinities at steady state. In fact, a clear fundamental understanding of reaction affinities is essential for understanding these methods. IUPAC defines the reaction affinity A as the derivative of the total free energy G with respect to a reaction coordinate, or the extent of the reaction, ξ, such that A ξ = .
G d d Nevertheless, there are conflicting definitions of reaction affinity in the literature where it is common to see the reaction affinity defined as A = −∆G rxn [10], where ∆G rxn is the free energy change in the system due to a reaction. Other non-standard definitions in the literature directly employ the standard chemical potential µ i 0 or the concentration-dependent chemical ( ) [19], where the signed stoichiometric coefficients ν γ = − i i for reactants and ν γ = + i i for products are used. In these definitions of reaction affinity, the loose use of calculus (a derivative is not a difference) is used to associate reaction affinity with free energy changes. Moreover, adding to the confusion is the statement that µ ν ∆ = ∑ G i i rxn . While this is true for a closed system characterized by constant total number of particles (in the limit of large numbers), many reactive systems of high interest, such as biological systems, are open systems. In an open system, the chemical potential µ i is defined as a function of the average count or concentration of the respective chemical species [20]; that is, the chemical potential is a parameter-a constant-rather than a variable. In this case, the value given by µ ν ∑ i i is also a constant; it does not vary with reaction conditions. But the change in the reaction free energy does vary with reaction conditions; it is only approximately constant if a reaction doesn't change the abundance of the reactants or products significantly, as in the case of the thermodynamic limit of large concentrations at steady state. If this is the case, then the usual kinetic description of the reaction rate can be formulated to include a free energy term. In the case of reaction scheme A,  exists when ∆G rxn is highly favorable and in that case the approximation is useful [8,9]. But this is an extreme case and is not appropriate for the variety of reaction conditions in many systems. To see that this is the case, one only need to consider the initial phase of a reaction when only reactants for the forward reaction are present. In this case, for the forward reaction free energy ∆G 1 and a reverse reaction free energy ∆ − G 1 , equation (4) results in a finite reverse reaction rate because ∆ = −∆ − G G 1 1 and ∆G 1 has a finite value. One could argue that if the reactant concentrations for the reverse reaction are at zero, then the exponent of the free energy change for the reverse reaction should be zero (i.e. ∆ − G 1 would have an infinitely large positive value). However, the logic of this argument results in the contradiction that ∆ ≠ −∆ − G G 1 1 and therefore free energy is no longer a function of state. Even at steady state, when the number of particles of the reaction intermediates drop below ~10, equation (4) will not produce the correct results. A single firing of a reaction in this case can easily result in changes in the abundance of chemical species of 10-100%. This discussion suggests that a microscopic function related to the free energy change is needed that does not require any assumptions about linear relationships between rates and thermodynamic driving forces.
In this report, we first review the relationship between reaction affinity and free energy for isolated reactions using concepts from statistical thermodynamic integration to explicitly demonstrate the relationship between kinetic formulations of rates that use rate constants and a statistical thermodynamics formulation that uses chemical potentials. (Readers intimately familiar with the relationship between reaction affinities, free energies and multinomial statistical models can read ahead to the discussion of coupled reactions following equation (12).) Subsequently, a theorem for coupled reactions is presented based on chemical potentials which can provide relative rates for all sequentially coupled reactions under all conditions-steady state or not-obviating the need for rate constants in many applications. However, if just one rate constant is known for the system of sequentially coupled reactions, then exact rates can be obtained, as well. The only requirement is that steady state levels of reaction intermediates be measurable and available to formulate chemical potentials. Unlike previous thermodynamic formulations, the approach is capable of modeling equilibrium, near-equilibrium and far from equilibrium conditions. We discuss the current work in the context of generalized detailed balance, and show that the work described herein extends the applicability of generalized detailed balance. The ability to model the dynamics of systems for which accurate rate constants are not available, such as biological systems, would be a considerable advance for predictive modeling in biology. At macroscopic scale, this approach is shown to be the rescaled deterministic rate law based on a rigorous definition of reaction affinities and a coupling constant is introduced that couples the rates of independent reactions. Furthermore, in the case that steady state concentrations are not available, it is shown that a maximum entropy assumption results in both a thermodynamically and kinetically optimal system in which the rate of entropy production is maximized.

Review of reaction affinity and closed systems
The IUPAC Gold Book defines the reaction affinity as, where ξ is the extent of a reaction and G is the free energy, usually taken to be the Gibbs energy. We will use instead G to indicate a generalized free energy function that can be specifically defined by ensemble parameters such as total number of particles (N), volume (V) and temperature (T). During an incremental increase in a reaction by ξ d , the amount n i of a substance i changes by where ν i is the signed stoichiometric coefficient such that γ ν = i i . Over the course of a reaction, an initial amount of substance changes from n i,0 to = +∆ n n n During the course of the reaction, the free energy change is given by, For the extent of the reaction, we use the accepted criteria for thermodynamic integration [21,22] such that ξ varies from ξ = ′ 0 initial (no reaction) to ξ = ′ 1 final representing a complete reaction in the sense that a stoichiometric unit of products is produced from a stoichiometric unit of reactants. When the extent of the reaction is greater than one, then more than a stoichiometric amount of material is processed through a reaction. Since the free energy is a state function, the actual path of the integral does not matter.
The reaction affinity as a function of concentrations or counts can be found by formulating the free energy as the multinomial Boltzmann distribution and explicitly taking the derivative of the free energy with respect to the extent of the reaction or reaction coordinate, ξ. The multinomial Boltzmann distribu-tion describes how a mixed population of molecules is distributed among M chemical species i, each with a count of n i , such that ∑ = n N i T . The species are not equally likely, but are distributed according to the Boltzmann criteria of the molecular partition func- is the free energy of formation of compound i averaged over all solvent configurations and k T B is the ambient thermal energy (k B is Boltzmann's constant and T is temperature). The justification for using the Boltzmann criteria is similar to that used by De Decker et al, that the heat released due to chemical reactions is either moderate or averaged over, such that microscopic processes like velocity distributions of individual atoms and molecules are decoupled from the statistics of the counts [23]. The microscopic analogue G of the free energy is then, where Stirling's approximation is used in equation (6). G is a microscopic analogue of the macroscopic Gibbs energy G because it does not involve averaging the counts n i over all configurations. (The macroscopic Gibbs energy G is a function of total number of particles (N total ), temperature (T) and pressure (P) and is an average over all microstates at fixed N total , T and P. Thus, once initial transients are averaged out, the macroscopic free energy is a constant for a closed system.) Taking the derivative of the microscopic free energy G with respect to the extent of the reaction gives, Here, N T is an ensemble parameter-a constant-but for an open system ξ N d d T / is small but not generally zero. In the,next section we deal with open systems; here, for the sake of simplicity we will assume Next, we substitute the identity which is unfortunate because this makes the chemical potential a variable rather than a parameter which can lead to confusion when discussing an open system where the chemical potential is strictly a parameter (see below). Likewise, one could argue that the open system parameter µ i should be referred to explicitly as an average chemical potential. While these distinctions are clear in discussions of statistical thermodynamics and ensembles, they are less clear in discussions of metabolic modeling. Regardless, the reaction affinity can now be written as, where the approximation is again due to the use of Stirling's approximation in equation (7). To avoid later confusion with the chemical potential in an open system, we will express the reaction affinity in terms of the equilibrium constant and counts (or similarly concentrations). Taking the exponential of the reaction affinity and writing the counts n i as an explicit function of the extent of the reaction ξ at the initial condition ξ = 0, The constant K 1 is the equilibrium constant for reaction 1, defined as the ratio of the molecular partition functions q i of the species i involved in the reaction, each raised to the power of its stoichiometric coefficient in the reaction [24]. It is important to note that the partition functions … q q , , M 1 are the exponent of standard chemical potentials µ µ … , , M 1 0 0 (scaled by the available thermal energy), which in this case are averaged over all lower degrees of freedom-solvation, translational, rotational, vibrational and electronic.
In the case of reaction scheme A, the reaction affinity becomes, In contrast to the reaction affinity, the microscopic analogue of the free energy change is the log odds of the densities of two different states, which for an reactant state at ξ = 0 and a product state at ξ = 1 its exponential is [25], where the approximation is again due to the use of Sterling's approximation. To demonstrate that equation (9) is a microscopic free energy change, the free energy change of a reaction can be obtained from carrying out the integration in equation (5), Recognizing that the approximation in the second step is due to the earlier use of Stirling's approximation in equation (7), we can rewrite the last equality as an exact relationship, Again recognizing that the equilibrium constant K 1 is the ratio of the molecular partition functions of the species involved in the reaction, each raised to the power of its stoichiometric coefficient in the reaction, for reaction scheme A we get, In the last expression the unsigned stoichiometric coefficients γ ν = i i have been substituted for the signed coefficients. It is clear from comparison of equations (8) and (9) that the rigorous interpretation of the ratio of two mass action functions in equation (3) is the reaction affinity and not the microscopic free energy change for the reaction. If the value of the reaction affinity is constant when averaged over many states, as is the case for large systems at steady state and assuming real valued counts or concentrations, then the average value of the reaction affinity and the free energy change are equal in value. More generally, the reaction affinity can completely describe the relative dynamics of single (isolated) reactions even for nonsteady state systems. The reaction affinity, expressed strictly in terms of standard chemical potentials as parameters, is

Open systems and sequentially-coupled reactions
In biology and many other fields reactions rarely occur in isolation. Instead, multiple reactions may be coupled together that transform initial reactants into final products. Consider the simple coupled reaction system, where A, B and C are chemical species and γ i j , is the stoichiometric coefficient for species i of reaction j.
For an open system, the microscopic free energy is given by a multinomial Boltzmann distribution similar to equation (6) but without any restriction on the total number of particles. In this case, the derivative of the total number of particles in equation (7) must be included in the integral of equations (7) and (10). If the total number of particles before the reaction (ξ = 0) is N T 0 and after the reaction , then the respective reaction affinity and free energy change are, Here the symbol Ξ is used to represent the microscopic free energy of an open system and ∆Ξ to indicate a change in the microscopic free energy with dependence on the total number of particles. In the equation for the free energy change, is negligible for systems with a large number of particles undergoing a single reaction. For the rest of the discussion, we will assume that the total number of particles does not change simply for convenience and clarity, and will consequently use formulations analogous to equations (8) and (9) instead. In order to model coupled reactions in a manner analogous to the use of reaction affinities for isolated reactions, a function is needed for the odds of one reaction relative to the next reaction. That is, for reaction scheme B information is needed on the relative rate of the first reaction to the second reaction. Taking the ratio of the mass action kinetic rate laws for reactions 1 and 2 in scheme B, a term analogous to equation (8) is needed that couples sequential reactions, Here, K K and 1 2 are rate constants that are functions of the volume of the system such that for K , and likewise for K 2 . Here K is used instead of k simply for convenience so that rate equations, thermodynamic functions and probability densities discussed below are all based on counts. The challenge in developing the function expressed in equation (13) is that K 21 cannot be determined from equilibrium measurements or calculated from free energies of formation. Like rate constants, the proposed coupling constants K 21 are system-dependent and have no system-independent standard value. Next, we demonstrate how nonequilibrium steady state measurements can be used to determine coupling constants and how sequentiallycoupled reactions can be modeled using chemical potentials instead of rate constants. In all cases discussed below, the extent of each reaction i is ξ = 0 i and so we will drop the explicit notation for the extent of the reaction.
For reactions shown in scheme B, the deterministic time dependence of each chemical species is governed by the usual set of ordinary differential equations. For species B this time dependence is, The deterministic rate equation can be rescaled in time by dividing by and replacing the ratios of rate constants with equilibrium constants, are equilibrium constants such that the second term in each parenthetical expression corresponds to the reaction affinity for each reverse reaction. The factor K K / , equation (15) is not useful for modeling systems of coupled reactions without rate constants. However, this coupling term can be determined by solving equation (15) at steady state. For the reactions in scheme B at steady state and with steady state concentrations denoted by SS, Equations (15) and (16) The parameters that are needed are the steady state concentrations and equilibrium constants K −1 and K −2 (or equivalently the molecular partition functions q A , q B and q C ). With these parameters in hand, it is easy to determine the coupling constants for the reactions using chemical potentials, The last equality is obtained by substituting for the definition of the chemical potential for an open system. There are many biological systems for which the stoichiometric coefficients are all unity. In this case, the coupling constant has the simple form, The rate equation (14) expressed entirely in terms of chemical potentials and in terms of time relative to τ 1 is then, where again the extent of the reaction is taken to be ξ = 0 for all reactions such that the counts n i are the instantaneous counts for each of the chemical species i. The rate equation can be expressed more succinctly in terms of equilibrium and coupling constants, where again K −1 and K −2 are the equilibrium constants for the reverse reactions in scheme B. Equations (16) and (19) are the key concepts of this report. It will be demonstrated in the Results that for sequentially coupled reactions, using these equations with steady state concentrations or chemical potentials in lieu of rate constants enables one to determine the relative time dependence of the chemical species. In the next section, we generalize these equations to an arbitrary chain of sequential reactions.

General form of coupling
A set of sequentially coupled reactions can be represented by, where Z v is the product of the vth reaction with a stoichiometric coefficient of γ ν ν Z , . For the stoichiometric coefficient γ ν ν Z , , the first subscript ν Z indicates the chemical species and the second subscript ν indicates the reaction. Equation (19) can be generalized using the last two reactions of scheme C as a general case. The coupling term for the vth reaction can be derived in terms of the previous reaction ν − 1 ( ) from the steady state condition, as follows. The relative time dependence of the counts ν− n Z 1 of the reactant of the vth reaction, − Z v 1 , is given by, where ν ν− c , 1 is the coupling term that will be derived. In equation (20)  such that, similarly to equation (16), where SS again indicates that the values of ν− n Z 2 , ν− n Z 1 , and ν n Z were obtained at steady state. The numerical value of ν ν− c , 1 at any state J is easily determined, The sequential nature of the reactions suggests an iterative process, in which the coupling between the first reaction in the sequence and the vth reaction is given by, Consequently, the change in the counts of the product of the ν − 1 th ( ) reaction, ν− n Z 1 , when rescaled by the rate of the first reaction in the sequence is, Using this formulation, all reactions in a reaction network can be rescaled by the transition probability of any specific reaction. In principle, then knowledge of any one rate constant can be used to derive the absolute time-dependence of the entire system.

Simulations of mass action using coupling constants in lieu of rate constants
Like equilibrium constants, coupling constants between reactions are the ratio of the rate constants of the respective elementary reactions. Unlike equilibrium constants, however, coupling constants between reactions such as K 21 are system dependent just as are rate constants and chemical potentials in non-equilibrium steady state systems. For convenience, we refer to equations (18) and (19) as differential equations using coupled reaction theorem. For the reactions of scheme B, figure 1 shows the comparison between a simulation solving the stochastic equivalent of equation (14) and one solving for stochastic equivalent of equation (19). Both simulations used the Gillespie stochastic simulation method [26]. In the kinetic stochastic simulation, reaction/transition probabilities were determined using rate constants and counts, while in the stochastic thermodynamic simulation reaction probabilities were determined using chemical potentials and counts of species. Reactions were selected randomly based on the formulation of transition probabilities shown in table 1. In these simulations, the absolute time is not modeled-all results pertain to the characteristic time of the +1 reaction in scheme B. However, if just one rate constant is available then the coupled reaction theorem can be used to reproduce the exact same timedependent trajectory as the trajectory governed by the (un-scaled) ordinary differential equation (14). The top row in figure 1 compares the steady state trajectories of the reaction intermediate B from stochastic simulations using the coupled reaction theorem (blue) with trajectories from stochastic kinetic simulations (red). Simulations were carried out under a range of conditions, from mildly non-equilibrium conditions with a total free energy change from initial reactants to final products of −2 kcal mol −1 , to strongly non-equilibrium conditions with an overall free energy change of −20 kcal mol −1 . The steady state solutions of the deterministic (real valued) ordinary differential equation are shown as dotted lines for comparison. When using the same set of random numbers, the trajectories from the coupled reaction theorem simulations are exactly the same as that for the stochastic kinetic simulations (the trajectories are offset by +5 counts for clarity), indicating that fluctuations around the steady state are identical. We can also demonstrate that the reaction probabilities are proportional, which is to be expected from the mathematical equivalence of equations (15) and (19).
The use of the coupled reaction theorem is valid away from the steady state, as well. Shown in row 2 in figure 1 are the transient decays from non-equilibrium states to the steady states obtained by setting the initial counts of the intermediate species B at a multiple of the steady state count. The transient is then the decay of the count of the intermediate species B to its steady state level. As expected, both the stochastic kinetic simulations and the coupled reaction theorem simulations produce exactly the same trajectory when using the same set of random numbers. In the case of rows 1 and 2 in figure 1, the forward rate constants for the two reactions differ in scale by four orders of magnitude (K K = − / over a large range. Row 3 in figure 1 shows the average steady state concentrations calculated for both coupled reaction theorem and stochastic kinetic simulations when K K 1 2 / is varied over eight orders of magnitude (10 −4 to 10 4 ) in 100 000 different simulations, each represented as a point in the plot. As the ratio K K Steady state concentrations, however, are less sensitive to variations in the rate constants than the rate of material flow through the system of reactions. The net flow of material through reaction 1 in scheme B during a simulation with n steps timesteps is given by,

Method
Transition probability Coupled reaction theorem of K 1 and K 2 . Once again the coupled reaction theorem and stochastic kinetic simulations produce exactly the same results when the same set of random numbers are used to select which reaction to fire. In fact, the correlation between the trajectories from stochastic kinetic and the coupled reaction theorem simulations across all values of K K 1 2 / is 1.0 within the numerical precision of the software.
As mentioned above, if just one rate constant is available then coupled reaction theorem simulations can reproduce the exact time-dependent trajectory as the trajectory governed by the (un-scaled) ordinary differential equation (14). Measurement of K 1 allows one to formulate the time dependence of the reaction intermediates such as B as (see equation (19)),

Simulations in relative time
More conveniently, simulations can be carried out in relative time using equation (19) and then adjusted post-simulation to absolute time. Furthermore, relative to equations (23) or (14), the use of equation (19) can accelerate the convergence of the non-equilibrium transients according to the time scale of the rescaling reaction. If the dynamic equations are rescaled by the reaction with the fastest dynamics or highest probability, then the values of the respective time derivatives for the set of differential equations governing the dynamics of the system is reduced accordingly. As a consequence, simulations may converge to the steady state in a fewer number of steps. Regardless, the simulated time-to-steady state will be equivalent when time rescaling is used in this manner since reaching the steady state depends on the time scale of the reaction with the lowest probability or slowest dynamics. To illustrate this time dependence, consider a pair of reactions similar to those in scheme B but catalyzed by a pair of enzymes, as shown in figure 2. Enzyme 1 (E 1 ) binds substrate A and produces product B, which is released and in turn bound by enzyme E 2 to produce product C. Starting from a highly non-equilibrium state consisting of only substrate A and enzymes E 1 and E 2 , initially the fastest reaction will be the binding of A to E 1 to form the complex E 1 A. The slowest reaction will be that producing the final product, C. Using a kinetic simulation based on rate constants and coupled reaction theorem simulation based on chemical potentials and rescaling, the dynamics of the unscaled and time-scaled system over a one second window are shown in figure 3(A) for the first enzymatic reaction and in figure 3(B) for the second enzymatic reaction. The differential equations used to model the reactions are analogous to equations (19) and (14), and are provided in the supplementary information (stacks.iop.org/PhysBio/14/055003/mmedia) along with all of the details of the simulation in the form of a browser-based computational notebook. A stiff, adaptive solver using was used to solve the differential equations [27].
Although both simulations start from the same set of concentrations, the simulation using the coupled reaction theorem and rescaling converges to the steady state more rapidly under all reaction conditions that were tested. The kinetic simulation spends a considerable amount of time simulating the fast dynamics that occur below 10 −4 s (figures 3(A) and (B)), while the simulation using the coupled reaction theorem effectively accelerates the convergence of the fast dynamics and spends relatively little time below the millisecond time scale. Regardless, the steady state is reached in the coupled reaction theorem simulation at approximately the same cumulative time point as the kinetic simulation because the dynamics of the slowest reactions are unchanged in the vicinity of the steady state.
That the rescaled ODEs resulting from implementing the coupled reaction theorem require less computational work to reach steady state solution is implied by figure 3(C) which shows that the rescaled ODE's require fewer number of time steps to reach the steady state than the full ODE using rate parameters. Specifically how much computational work is reduced and under what conditions will be a topic of future reports where it will be tested on a sufficiently complex system of equations.

Course graining using summary reactions
In addition to being able to produce the correct dynamics without the need for rate constants, the coupled reaction theorem has an advantage over kinetic formulations of reaction dynamics for multiscale modeling of the steady state: course-graining the dynamics such that the equations are 'telescopic' is relatively easy. One can zoom in or out of the details of the reaction system as needed, which can be a considerable advantage for modeling [13]. If one is only interested in steady state phenomena, it is not necessary to model elementary reactions; summary reactions can be modeled using chemical potentials instead and still produce the correct steady state dynamics. The course-grained dynamics will be the composite dynamics of the collapsed system, however.
To demonstrate the ability to collapse reaction schemes to coarse-grained summary reactions, one only needs to consider a Kirchoff's loop relationship for cycles [28,29]. Consider the first reaction cycle shown in figure 2 from reactant A through reactions 2-4 to product B and then through reaction −1 back to initial reactant A. The free energy change for traversing around a cycle is zero. In the deterministic limit (large number of particles) and at steady state, the relation between rate and free energy is  where R is the gas constant, + r is the forward reaction rate and − r is the reverse reaction rate [11,30,31]. For such a cycle at steady state, For an enzyme catalyzed reaction of a single substrate, + r 2 , + r 3 and + r 4 correspond to the rates of binding of the substrate, the catalytic conversion of the substrate to product, and the rate of release of product, respectively. The ratio + − r r 1 1 / then corresponds to the ratio of rates of the combined steps. The dynamics of the full set of reactions (reactions 2, −2, 3, −3, 4, −4, and 6, −6, 7, −7, 8, −8) and the summary reactions (reactions 1, −1, 5, −5) in figure 2 are compared in figure 3(D), which shows the rate ratio of each reaction as the system approaches steady state, as well as the product of rate ratios in equation (24). In this figure, comparison of the product of the ratio of rates of reaction 2, 3, and 4 (Rxn 234, blue solid line in figure 3(D)) can be compared with the ratio of the rates for reaction 1 (Rxn 1, blue dashed line). While the combined trajectory for reactions 2, 3, and 4 (in the form of the product of the ratio of the rates) start at 10 −8 s, the trajectory doesn't converge until 10 −4 s. In contrast, the trajectory of the ratio of forward and reverse rates for the equivalent reaction 1 takes minimal time to converge in the first few time steps and has converged by 10 −6 s. Likewise, the combined trajectory for reactions 6, 7, and 8 (Rxn 678, blue solid line) start at 10 −8 s and converges by 10 −4 s, while the trajectory for equivalent reaction 5 (Rxn 5, blue dashed line) converges immediately in the first few time steps at 10 −6 s.
Compressing reactions 2, 3, and 4 into summary reaction 1 and reactions 6,7, and 8 into summary reaction 5 are the mass action equivalent of course-graining in atomistic simulations: the degrees of freedom are averaged over and compressed to fewer degrees of freedom, enabling the removal of high-frequency dynamics from the trajectories. Since the coursegrained summary reactions 1 or 5 do not model the dynamics of enzyme binding, the rate ratio converges to the steady state value rapidly.
That the rate ratios can be predicted using coursegrained summary reactions and without the need for rate constants can be a considerable advantage for modeling biological systems. To model the simple two-substrate, two-product enzymatic reaction for the conversion of dihydrofolate to tetrahydrofolate, Fierke et al, demonstrated that the kinetic scheme involves 13 reactions [32]. In the coupled reaction theorem formulation, each pair of coupled reactions has the same canonical form shown in equation (19). One can model either the elementary reactions or a summary reaction that describes the overall phenomena as long as the appropriate steady state concentrations are available.
Using a coupled reaction theorem formulation, measurement of the steady state levels of metabolites and proteins as they exist unbound in the cytosol can be used to derive information on the dynamics of the individual and overall steps of the enzyme-catalyzed reactions, including rate constants.

Branched reactions
For branched reactions, however, steady state concentrations alone are not sufficient to derive the dynamics. Consider the branched reaction where a product of the first reaction A B is then the reactant for two reactions occurring in parallel, When the concentrations of each species is large enough, the deterministic rate law for species B is, Using the approach described above and solving for rate ratios at steady state results in two unknown coupling ratios and one equation, which cannot be solved without additional information. However, if either the flux = − = − − − r r r r r r or 1,net 1 1 3 ,net 3 3 through reactions 1 and 3, respectively, is measured using isotopic labeling assays (metabolic flux analysis [33]), then this additional information allows for determination of the necessary coupling parameters,   (19)) affect the reaction dynamics-in particular, the net rate through a reaction (flux) and the energetic cost to a biological cell. As can be seen in figure 1(D), changing the ratio of the rate constants K 1 and K 2 can dramatically affect the flow of material through the reactions. (Again, there is a many-to-one relationship between the ratio K K 1 2 / and the overall rate because each value of K K 1 2 / can corresponded to multiple combinations of K 1 and K 2 .) In particular, a high or low value of the ratio K K 1 2 / does not lead to a decrease/increase in overall rate; instead the optimum value appears to be closer to 1.
The top row in figure 4 demonstrates the relationship between the coupling term K For each plot from left to right, the overall free energy change for the system A B C changes from −2.0 to −20.0 kcal mol −1 . Within each plot, the net rate of production of product C is a sigmoidal function of the coupling term, specifically the rate parameter k 1 , with the highest rates occuring when ,1 ( ( )/ ) / is equal to or greater than 1.0. As the driving force on the overall reaction is increased from left to right in each of the four plots, the range of values for the coupling term converges to 1.0, as predicted from equation (16) when the reverse reactions −1 and −2 become highly unlikely.
As can be seen in the middle row of figure 4, the average value 1.0 for the coupling term has particular significance-it is the maximum entropy change (or production) steady state out of all the possible steady states. The definition of entropy production used here is the the entropy change in going from initial products to final reactants due to reactions (rxns) [25,34], where C J is the configuration of state J, δS i is a change of state due to reaction i, δ is the time independent state-to-state transition probability given by, The total free energy of the system is also represented in the figures by color-coding. The state of lowest free energy is also the state of maximum entropy production. Thus, the maximum entropy production steady state is the thermodynamically most efficient steady state (least heat dissipation) and is characterized by being the steady state of lowest total free energy and highest thermodynamic likelihood. By comparison of the top and middle rows, and especially for a nonequilibrium driving force of 2 kcal mol −1 , it is apparent that it is possible to maximize the net rate but at the cost of thermodynamic efficiency. That is, while a somewhat higher net rate can be achieved, the system must be maintained at a higher energy. However, at even moderate thermodynamic driving forces of 5 kcal mol −1 for the overall reaction, there is no speed advantage-the net rate is maximized very close to the thermodynamic optimum. As long as both reactions individually have / . (top row) Net rate of formation of final product C from the reaction system shown in scheme B according to equation (25). As the driving force increases from −2.0 to −20.0 kcal mol −1 across the top row, the value of the coupling term at which the rate becomes maximal converges to 1.0. (middle row) Entropy production (equation (26)) as a function of the coupling term. The total free energy of the reaction system, in kcal mol −1 , is represented by colors, as well. The maximal entropy production and total free energy both occur at K = γ γ n n 1.0 / . (bottom row) Entropy production rate (equation (27)) as a function of the coupling term. The entropy production rate is also maximal when the coupling term is 1.0. moderate non-equilibrium driving forces, there is no speed advantage to tweaking rate constants. Rate and thermodynamic efficiency are meaningfully combined in the concept of entropy production rate, hinted at by Lotka [35,36] and proposed by Schrödinger and many others since [37] as the operating principle for biological systems. The definition of the entropy production rate used here is, net r xn (27) As can be seen from the bottom row of figure 4, the entropy production rate is also maximized when K = γ γ n n 1.0 , as are the entropy production and total free energy of the system. This is significant in that it suggests that on the larger scale of biological systems there does not necessarily need to be a tradeoff between growth rate and the thermodynamic efficiency of growth. Furthermore, the bottom row of figure 4 shows that there is a dramatic decrease in steady state entropy production rates as one moves away from the thermodynamically optimal state in which If one assumes that adaptation favors organisms with greater entropy production, then one can create reasonable models without even the need for steady state values of reaction intermediates by setting the coupling term to the value that maximizes the entropy production rate.
The definition of entropy production used above differs from that used recently by De Decker et al, in describing entropy production rates for mass action processes [23], where T is the temperature. Besides the difference in units, one important difference is that the definition of entropy production rate in equation (27) uses the overall rate going from initial reactants to final products of the system rather than the rate of individual processes. One definition of entropy production rates is not necessarily more correct than any other, but rather, like any descriptive statistic, the appropriate entropy production rate formulation depends on the phenomena one is characterizing. For biological systems undergoing natural selection in which one wants to characterize competition among species, it is the overall rate of entropy production via growth that is important and not necessarily the entropy production rates of sub-processes.

Discussion
In the work described above, we show that the state of the system at steady state is sufficient to define the material and energy fluxes through a system of serially coupled reactions. This is particularly important for biology in that both steady state concentrations and fluxes can be determined for many chemical species. Related to this work, Zia and Schmittmann developed a framework for generalized detailed balance which extended the concept of detailed balance to the non-equilibrium domain [38]. The idea is that any steady state is uniquely characterized by the pair * * Pr V , { }, where Pr is the probability of a configuration = … C n n , , M 1 { } and * Pr is the value at the steady state configuration (the most probable state in a system at steady state). * V is the set of probability currents which are proportional to the net steady state fluxes; that is, an element v i of V is such that ∝ − − v r r i i i , the difference of the forward and reverse reaction rates for reaction i. Equilibrium is the special case of * Pr , 0 { }. Zia and Schmittmann demonstrated that the pair * * Pr V , { } uniquely characterize a steady state for both equilibrium and non-equilibrium systems. In contrast, for the special case of sequentially coupled reactions described above, we demonstrated that configuration, = … * * * C n n , , , } the set of counts/concentrations at steady state, alone is sufficient to uniquely characterize a non-equilibrium steady state in terms of relative time, and that the absolute time-dependence can be determined with knowledge of just one rate. The reason for this is that the standard chemical potentials of species in a reaction constrain the dynamics of forward and reverse rates of the same reaction, while the dynamics between reactions are likewise constrained by the observed steady state chemical potentials. The steady state chemical potentials as parameters for an open system determine the probability of any chemical species i, Equation (28) is directly related to the microscopic free energy of the system [25]. It is clear from the work described above that for branched reaction pathways, knowledge of and only a subset of * V is sufficient to characterize a unique non-equilibrium steady state. This begs the question, how many fluxes are needed to uniquely characterize a steady state in addition to the M chemical potentials/ concentrations?
Fleming et al [7] have recently addressed a this question-what conditions are necessary and sufficient for duality between unidirectional fluxes and the concentrations/counts of chemical species in systems of coupled reactions? Their conclusion was that for a stoichiometric matrix of dimensions M × P consisting of M rows corresponding to chemical species and P columns corresponding to chemical reactions, the rank of the stoichiometric matrix determines the number of reactions whose dynamics can be determined based on species concentrations alone. Surprisingly, they found that 26 of 29 genome-scale metabolic models had stoichiometric matrices that were full rank. Of the remaining three stoichiometric matrices of the models, all were less than full rank by only three or less. Since the steady state concentrations are used to set chemical potentials for use as parameters in an open system, the implication for the work presented above is that the relative reaction dynamics could in principle be simulated using chemical potentials in place of rate constants for as many reactions as there were metabolites in the full rank models. This suggests that only M−P fluxes are needed in addition to the M chemical potentials/concentrations to uniquely characterize any steady state system without rate constants.
The demonstration that the maximum entropy production rate solution to the system of equations in scheme B is in many cases also the maximum rate solution prompts discussion on the nature of natural selection. Previously, a correlation between the standard entropy changes of overall reaction pathways and their activities has been demonstrated [39]. If metabolite measurements are available, the hypothesis of maximum entropy for steady states of evolutionary optimized systems is testable [40,41]. It is important to note however, that while entropy production may be maximized at steady state, steady state entropy production rates are constrained by the rates achievable by the respective enzymes [42].
From an evolutionary perspective, conservation of a specific rate constant would be very hard and would require strict error tolerances in replication, which would lead to decreased ability to adapt. More likely, fitness fluctuates around the state of maximum entropy production. However, it must be kept in mind that the coupled reactions of metabolism do not represent the total thermodynamics of the cell and it is even possible for an individual reaction entropy to decrease rather than increase. Feedback regulation of enzyme activities may occur for many reasons, such as regulating the production of metabolites to synchronize the cell's network of coupled reactions with those of the environment [43]. Regardless, predictive models of complex adaptive systems such as those found in biology and elsewhere that do not depend on hardto-measure parameters are urgently needed to accelerate research in systems relevant to medicine, climate and energy challenges. We have outlined a statistical thermodynamic framework that would enable such large-scale simulations. The approach does not rely on kinetic parameters, but rather on standard free energy values and metabolite concentrations for which robust measurement methods are being developed [44].