Reaction networks as systems for resource allocation: a variational principle for their non-equilibrium steady states.

Within a fully microscopic setting, we derive a variational principle for the non-equilibrium steady states of chemical reaction networks, valid for time-scales over which chemical potentials can be taken to be slowly varying: at stationarity the system minimizes a global function of the reaction fluxes with the form of a Hopfield Hamiltonian with hebbian couplings, that is explicitly seen to correspond to the rate of decay of entropy production over time. Guided by this analogy, we show that reaction networks can be formally re-cast as systems of interacting reactions that optimize the use of the available compounds by competing for substrates, akin to agents competing for a limited resource in an optimal allocation problem. As an illustration, we analyze the scenario that emerges in two simple cases: that of toy (random) reaction networks and that of a metabolic network model of the human red blood cell.

Introduction

The dynamics and thermodynamics of chemical reaction networks is a subject that goes back at least to [1]–[6]. Recent years have seen considerable interest in the problem at different levels, from the characterization of their mass-action kinetics [7] and of their stochastic thermodynamics based on the chemical master equation [8], [9], to the analysis of their non-equilibrium steady states (NESS) [10]–[13]. Besides their relevance for fundamental understanding, these approaches provide an important frame for the study of real biochemical systems, like genome-scale reconstructions of cellular metabolic networks [14]–[16]. The most basic information about these systems is usually encoded in the matrix of stoichiometric coefficients, representing in essence the (weighted) topology of the couplings between chemical species and reactions. With uncertainties about kinetic parameters and reaction or transport mechanisms often preventing large-scale kinetic approaches (with exceptions like the the metabolism of human erythrocytes [17]), the challenge at the simplest level is that of building stoichiometry-based predictive models of metabolic activity at genome scale. Much information on the organization of reaction fluxes in NESS can indeed be obtained from constraint-based models that rely on minimal mass-balance [18]–[20] or stability [21], [22] assumptions. Such descriptions revolve around pre-defined (and physically motivated) sets of local constraints, enforcing for instance mass balance at each metabolite node in the network. It would be interesting (and instructive) to derive the relevant local constraints as the result of a mathematical analysis of the dynamics taking place on the network, both to further clarify the assumptions behind the models and to highlight their limitations.

This work revisits the joint dynamics of concentrations and reaction fluxes in chemical networks in a statistical mechanics frame. In specific, we obtain a variational principle that relates fluxes (i.e. the average number of microscopic transitions per time per volume for each process) in NESS to the minima of a global function where stoichiometric coefficients and steady-state concentrations appear as parameters. This function is reminiscent of the Hamiltonian of a Hopfield model with Hebbian couplings [23], with stoichiometry playing the role of the patterns. An analysis of its physical meaning explicitly shows that reaction networks dynamically converge towards states where the use of the available compounds is optimized and the rate of decay of entropy production is minimized. The flux organization problem turns out to have remarkable similarities with that of optimal resource allocation by heterogeneous agents, as described e.g. by Minority Games [24], [25]. Systems of this type generically undergo a transition from an ergodic phase (the NESS is independent of the initial conditions of the dynamics) to a non-ergodic one when the ratio between the number of reactions and that of chemical species is changed. In our case, the two regimes are described by different sets of local constraints. We shall first explore this scenario in toy “random” reaction networks where such a transition can be fully analyzed numerically. Then a simple real system will be considered, namely the reduced metabolic network model of human erythrocytes. Finally, we shall discuss the relevance of these results for the quantitative analysis of cellular metabolism.

Analysis

The Variational Principle

We consider an open system enclosed in a volume V (a reactor or “cell”, for brevity), formed by M distinct chemical species (labelled ) that can be processed by N distinct reactions (labelled ) at fixed temperature T, pressure, ionic strength and pH. The reaction stoichiometry is described by the coefficients , with the convention that negative (resp. positive) coefficients identify substrates (resp. products) in the ‘forward’ direction of reaction i. Individual reaction events occur stochastically with rates (number of events per unit time) proportional to the substrate concentrations, which vary in time. At stationarity and for ideal systems, the Gibbs energy (GE) change per mole associated to each reaction i, i.e. [26] (where R is the gas constant, the intracellular concentration of species , a reference concentration and the GE change in standard conditions at concentration ) characterizes its distance from detailed balance. Specifically, the forward () and reverse () fluxes of i, i.e. the average number of transitions per time per volume, satisfies the relation with [10]. At equilibrium, and for each i.

We focus on the time evolution of the internal composition of such a chemical reactor. The dynamics of this system is driven essentially by two factors: (a) the fact that molecules can stochastically enter or leave the system, and (b) the fact that reactions events occur (stochastically) inside the system. Reasonably, then, the NESS will depend (a) on the rate at which molecules can cross the system’s boundaries per unit volume (intake or outtake fluxes, denoted by ), and (b) on the rates of the internal reactions (more precisely, on the average net number of microscopic transitions per time per volume, denoted by ). In the following we will show specifically that, given the boundary fluxes (characterizing the environment), over time scales for which the chemical potentials can be assumed to vary slowly (so that concentrations can be assumed to be roughly constant) the internal fluxes in a NESS minimize the functionwhere denotes the concentration of species .

Assume that at time the system is characterized by molecular populations (number of molecules) . Consider a time interval of size and let denote the net number of transitions of reaction i that take place between time t and time . The latter is a random number governed by a probability law (which shall depend, e.g., on the concentration of substrates) that we leave unspecified for the moment. The corresponding variation of ’s is given bywhere stands for the net (random) number of molecules of species taken in (if ) or out (if ) between time t and time . Taking V to be fixed, the time-evolution of concentrations is simply

We now focus on the quantitywhose change between times t and is given bywith initial conditions . Equation (3) contains sources of stochasticity in the ’s and the ’s. As a consequence, the “macroscopic” variables and will fluctuate stochastically as well. Our aim is to characterize the steady state(s) of (6). We make the following simplifying assumptions:

A1. Molecular populations are large enough to allow us to treat as a continuous variable. This is generically assumed to be the case in real cells, although e.g. in E. coli the number of copies of certain small molecules can be as low as a few tens (corresponding to a concentration of the order of 10 nM [27]). The effects induced by molecular noise can be non trivial [9] and accounting for it might alter the emerging picture [28].

A2. The quantity is small (i.e. the chemical potential , with the standard chemical potential, is roughly constant), so that can also be taken to be small.

Under the above assumptions the right-hand side of (6) is easily linearized to yield

This equation highlights the way in which concentrations affect the time evolution of the system and, in principle, one would now need to analyze the coupled system formed by (4) and (7). However, for simplicity, we replace with some time-independent limit value . This approximation can only be justified as long as one considers evolution over time scales shorter than . If however concentration changes are sufficiently slow (in agreement with homeostasis) it is reasonable to expect that it will hold over time scales much longer than those required to reach a NESS. With this, (7) takes the formwhere the “couplings” are defined as

The steady state of (7) can now be obtained straightforwardly by dividing both sides by and averaging over time. This gives

Note that is the net flux of reaction i (the average net number of microscopic transitions per time per volume) while represents the uptake of species (the average number of molecules of species per time per volume entering or leaving the cell). Definingwe therefore havewith

This implies that, on average, the stochastic dynamics of reactions collectively minimizes the function H so that the NESS of (12) correspond to the solutions of the minimization problemconstrained by the uptake values and by the fact that the ’s are assumed to be bounded by enzyme availability, so that .

Following [24], [25], we can characterize the behavior of fluxes in NESS by noting that the solutions of (12) for (which are expected to be linked to the Gibbs energy change of reaction i in a NESS) are generically of the form with , since is constant in NESS. In practice, the stochastic dynamics of for different reactions will be characterized by different values of , depending on the asymptotic behavior. If tends to a finite value as then . Recalling that the Gibbs energy change is related to the forward-to-reverse flux ratio through the detailed balance condition, this describes the case of a reaction whose microscopic transitions occur bidirectionally even for , i.e. such that the forward and reverse fluxes are both non-zero in the steady state. Specifically, its flux is determined by the condition . When , instead, increases or decreases linearly in time (after a transient) and diverges for so that, in the corresponding NESS, reaction i will be characterized by unidirectional microscopic transitions for . Note that if then H is indeed minimized by solving , whereas for the minimum of H is achieved by taking as large (if ) or as small (if ) as possible, i.e. (resp. ) if (resp. ).

Physical Meaning of H

For a start, notice that the change in is expressed in equation (8) as the sum of two terms. The first accounts for the coupling of i to other reactions in the network via shared compounds (the coupling coefficient is non-zero only if i and j have a metabolite with finite concentration in common). Consider a species and two reactions i and j such that ( j produces ) and ( is consumed by i). Then and a positive net advancement of reaction j will contribute to the increase of , see (8). Now looking at (6) it is reasonable to expect that the probability of observing a forward transition for reaction i between time t and time will be larger the larger (and positive) is (in agreement with the fluctuation theorem [8]; we shall see an explicit example of this later on). Therefore in this case a positive net advancement of reaction j will ultimately increase the probability of a concomitant advancement of reaction i. In other words, this situation favors the emergence of a positive correlation between reactions i and j. Likewise, if both i and j are either producing (, ) or consuming (, ) species , their coupling will tend to anti-correlate and . In this case the dynamics discourages over-production or over-consumption of a chemical species. A similar picture holds for the coefficients , which are related to the presence of sources (like nutrients) and sinks (e.g. outtakes) in the network. A non-zero acts as a force of magnitude proportional to that tends to polarize in a particular direction a reaction i that is stoichiometrically connected to a compound with . The favoured direction depends on whether is a source or a sink and on the sign of . This effect can propagate to other nodes connected to i if is sufficiently large. The role of the concentrations appearing in and is in essence that of modulating the strength of the couplings and of the forcing fields with the availability of the intermediate metabolites. Indeed, ’s get stronger when the intermediate compounds are present in smaller amounts, stressing the emergence of positive correlations between processes (if ) or the limits imposed by competition for a limited resource (if ). When the concentration of the intermediate is large, instead, the coupling gets weaker and i and j may become effectively independent. The impact of is understood along similar lines. In this way, the original bipartite network of reactions and metabolites can be re-cast as a system of interacting reactions with Hebbian couplings, as shown in Figure 1. This is strongly reminiscent of Hopfield models of neural networks. In such a scenario, finding the steady state(s) of (7) is equivalent to finding the ground states of a system of reactions interacting with “energy” H.

Figure 1

Mapping of a reaction network to a system of interacting reactions. Panel

(a) Toy reaction network with reactions represented as circles and chemical species as squares. Continuous, dashed, incoming drawn and outgoing drawn arrows denote stoichiometric coefficients and uptakes, respectively , , and . Panel (b) Reduced reaction network with couplings and “fields” given by (9). Continuous, dotted, incoming grey and outgoing grey arrows denote respectively , , and . For instance, , . Grey arrows are double-headed when the sign of h depends on the precise values of stoichiometric coefficients and uptake fluxes. For instance, the value of depends on the choice of the ’s and ’s, since the first term in the sum is negative while the second is positive.

From a physical standpoint, H quantifies the resource mis-usage by the network so that, by minimizing H, the system strives to reach states in which compounds are used as optimally as possible, given the initial conditions (that also account for the standard GEs and hence, to some degree, for the a priori reversibility), the stoichiometry, the available nutrients, the production goals, etc. Whether for a given network the minimum of H is zero or not then depends on several factors, including the bounds on fluxes and the specific form of the uptakes. Note that would imply that at stationarityi.e. that in the steady state a Kirchhoff-law type of scenario holds in which strict mass-balance conditions are satisfied for each chemical species. The network in this case organizes the fluxes so that consumption and production exactly match for each species and meet the nutrient availability and outtake requirements described by the vector . On the other hand, the physically relevant states with can more generally be thought to haveotherwise the system could e.g. consume a nutrient in excess of its availability. Both sets of conditions have been employed for the modeling of cellular metabolic networks. In particular equations of the type of (15) are the basis of the highly successful flux-balance-analysis (FBA) [19], [20], where biological functionality is included as an additional ad hoc constraint usually represented by the maximization of a specific score function (e.g. biomass production for E. coli in an optimal environment). The conditions (16) are instead reminiscent of Von Neumann’s model of reaction networks [29], where self-consistent flux states with a net positive production of intracellular metabolites can be allowed. This type of approach can be helpful in analyzing the metabolic capabilities of an organism and, if statistically robust production profiles emerge, in inferring (rather than postulating) cellular objective functions, with the idea that chemical species that are globally produced (e.g. amino acids) are employed in macromolecular processes (like protein formation) that are not encoded by the reaction stoichiometry. We remark however that within the above setting finding the NESS of the system means minimizing H, which obviously is a priori different from solving (15) or (16).

A broader physical insight can be obtained by studying how H relates to standard quantities used to characterize non-equilibrium behavior in reaction networks, such as the entropy production. Indeed the entropy production per volume of the system enclosed in the “cell” of volume V is given bywhere is the concentration of species and is its chemical potential. Taking time derivatives keeping in mind that at stationarity is constant, one easily finds that . Now since , . In addition, H is constant in NESS. It follows that for time scales shorter than (i.e. for time scales over which chemical potentials are roughly constant) one hasH is thus seen to play the role of the rate at which the entropy production changes in a NESS. If (i.e. if the NESS is flux-balanced), then must vanish as well (since ), leading to a state with constant entropy. If , instead, the entropy production is non-zero and decreases at the smallest allowed rate. Correspondingly, the entropy “slows down” quadratically over the time scales for which the theory holds. (Note that for sufficiently long times the entropy production can become negative if , implying that this limit may be unphysical.) In summary, the NESS that can be obtained are either characterized by , zero entropy production and hence constant entropy, or by , positive entropy production decreasing in time as slowly as possible and entropy increasing in time accordingly.

It is noteworthy that this scenario essentially characterizes the Lyapunov condition described in [30] for the stability of stationary states (Chapter 18). The Lyapunov function in our case can be computed explicitly, takes the form (13), its minimization bears a further physical meaning in terms of optimal resource allocation, and provides an equivalent description of the reactor as a system of processes interacting via Hebbian couplings. Indeed the quantity H appears in the Gibbs theory of thermodynamic stability for ideal systems. Let us consider a system at equilibrium with vanishing net fluxes and evaluate the effect of a perturbation that drives the system away from equilibrium. Such a perturbation corresponds to a (small) change in the turnover of each reaction i, which can be achieved e.g. by forcing a (small) non-zero flux through each reaction i for a time , so that . The free energy per volume associated to the perturbed state is easily found to be given byin agreement with the second law of thermodynamics. After turning off the perturbation, the system will relax back to equilibrium by minimizing G (i.e. H).

Results

A Toy Model: Random Reaction Networks

A simple algebraic argument allows to understand that, generically, a qualitative change in the solutions of (14) is expected to take place as one varies the ratio between the network parameters N and M. Let denote the number of reactions that remain asymptotically bidirectional, i.e. such that . The conditions that they must satisfy form a set of equations, M of which at most are independent (strictly speaking, the rank of the matrix equals that of the stoichiometric matrix , see (9); it is however always possible to eliminate dependencies in the latter so as to attain full rows rank). Neglecting the fact that variables are bounded, we can say that when the number of variables () exceeds that of equations the system is underconstrained and multiple solutions occur. So multiple steady states exist when . Note that is determined dynamically by (7) with initial conditions . This suggests that different choices of will lead to different steady states, i.e. that when ergodicity (i.e. independence of the steady state on initial conditions) won’t hold. A transition is thus expected to occur when .

Such a transition can be fully investigated within the following, simple model: the stoichiometric coefficients are independent, identically distributed random numbers (such that each reaction uses and produces a finite fraction of the possible compounds), and the dynamics of the network advances in discrete time steps of size . In specific, at each time t the ’s take on values stochastically in withi.e. with . Note that the above probability ratio is proportional to the ratio between the concentration of substrates and that of products of the reaction. Clearly, (20) does not allow to capture the temporal structure induced by the Arrhenius law, because by forcing each reaction to take place at every time step it neglects the fact that activation energies (and hence characteristic timescales) can differ significantly across reactions. It is however reasonable to think that the steady state will be unaffected by these transients. On the other hand (20) makes the theory considerably easier from a mathematical viewpoint and the final result for the steady state (which is the focus of the present work) more transparent. Setting for simplicity, so that for each i, one can follow [31] to derive the continuous-time limit of (7) for large N and M. The result is the Langevin processwhere is a Gaussian noise with zero mean. Clearly, time-averaging leads back to (12) with . Recalling that and noting that this implies where, by (20), , this in turn suggests the relationwhich explicitly links the thermodynamic driving force of a reaction to the (stochastic) dynamics of the quantity . The fact that the final result differs from the naïve intuition is a direct consequence of the stochastic fluctuations encoded in (20). Note that the dynamics (20) thus converges to NESS (minima of H) that are thermodynamically feasible, in agreement with the second law of thermodynamics. We can study the linear stability of (21) by settingwhere is a NESS and is a zero-average noise representing (small) perturbations to the trajectory. Reactions with , for which diverges, will be insensitive to . Hence it suffices to focus our attention on the response of reactions for which . To first order in , fluctuations are easily found to obey the condition

where . Now if all eigenvalues of the matrix are positive the dynamical system will be linearly stable as small perturbations occurred along the trajectories will die out in time. The term is clearly positive, hence it suffices to check that the smallest eigenvalue of the matrix is positive. Assuming that are independently and identically distributed, the spectrum of can be computed, in the limit with finite, using the results of [32]. For the smallest eigenvalue one finds , where and is a constant. Hence stability (and ergodicity) requires . When , instead, the dynamics will be sensitive to small perturbations and, reasonably, its steady state will be selected by the initial conditions . The marginal stability condition coincides with the rough algebraic estimate of the transition point made above, i.e. .

To illustrate the scenario underpinned by the above theory, we have simulated the dynamics defined byfor an ensemble of artificial reactors formed by M species and N reactions in order to analyze the dependence of its steady state(s) on . Stoichiometric coefficients were chosen randomly, so that with probability p and with probabilities , independently on i and . Similarly, for the boundary metabolites we took with probability q and with probability , independently on . For sakes of simplicity, we set for all (it is clear that the choice of does not affect the transition point where H vanishes; it only changes the value of in the ergodic phase). No prior assumption on reversibility was made, i.e. all microscopic transitions can initially occur in both directions for each process. The coefficients and are defined as in (9) and (11) except for the fact that is re-scaled by to ensure that the system is well-behaved when , while evolves according to (20). For consistency, H for this case is defined as

Results for are shown in Figure 2. One sees that initializing the dynamics from “equilibrium” conditions for each i one ends up in stationary states with for whereas for . Similarly, for while for , so that the critical point is indeed marked by the condition . When the dynamics starts from for each i, instead, the system ends up in a different steady state for , as signaled by the different value of the fraction of asymptotically unidirectional processes at stationarity, . In the phase characterized by , the steady state is unchanged. These results fully confirm the theoretical predictions derived above. We also note that (data not shown), if for each , i.e. if there is no exchange with the surroundings, the corresponding networks converge to a steady state in which for each i, and for each n, i.e. to chemical equilibrium. In other words, expectedly, boundary fluxes induce NESS.

Figure 2

Results for random reaction networks.

Average stationary values of H, fraction of asymptotically unidirectional reactions () and relative number of asymptotically bidirectional reactions versus obtained from (26) (with no prior assumption on reaction reversibility for an ensemble of random reaction networks with constructed as described in the text. Averages are taken over 200 realizations for each value of n. Unbiased i.c. (initial conditions) refers to steady states of (26) for ; biased i.c. instead correspond to for each i. Note the two phases with ( and () as predicted. The critical point coincides within numerical error with the point where . Finally, the phase with is non-ergodic: different initial conditions lead to different NESS, characterized by different values of .

Red Cell Metabolism

We now turn to a somewhat more realistic case in which nevertheless a full analytic study of the scenario underlying the minimization of H is possible, namely the standard reduced model of the human red blood cell (hRBC) metabolism, which includes the glycolytic and the pentose phosphate pathway only (21 reactions among 30 metabolites plus two ionic pumps, namely ATPase and NADPHase; see Tables S1 and S2 for the numerical details and the network structure used here). They operate by consuming glucose (GLC) in order to, respectively, maintain the osmotic balance through the sodium-potassium ionic pump (ATPase), and reduce the amount of free radicals through glutathione reductase (NADPHase). Moreover, at a simplified level, one can think that the final products of these pathways are, respectively, lactate (LAC) plus an exchange of K+ and Na+ ions, and CO2. The partitioning of GLC between the two pathways depends on the level of oxidative stress faced by the cell. Experimental estimates based on enzyme activities range from 70% or more in favor of glycolysis in unstressed conditions to 70% or more in favor of the pentose-phosphate pathway under oxidative stress [27]. We want to use the theory described here to estimate the range of variability of the fraction of glucose consumed by each pathway as a function of the oxygen concentration in the environment.

We start by assuming (reasonably) that the steady state of the hRBC metabolism is compatible with (i.e. flux balance). A straightforward analysis of the emerging equations reveal that only three of the 23 reactions are linearly independent: we choose the glucose uptake , the flux through the Rapoport-Leubering shunt (or through the enzyme 2,3-DPG mutase, DPGM, a key step that regulates the haemoglobin’s affinity with oxygen) and the flux through the pentose phosphate pathway (or through the enzyme glucose-6-phosphate dehydrogenase, G6PDH) . All fluxes can be written in terms of these. In particular, one finds

Now the variation of the extracellular concentrations of GLC, LAC, K+, Na+ and CO2 due to the operation of a single hRBC is easily seen to be given by

In turn, for the extracellular medium one has

Minimizing this at fixed and one finds thatwhere a and b are defined respectively as

One sees that if the concentration of CO2 is much larger than the others (implying ) and then so that the pentose phosphate pathway consumes roughly all of the glucose (the factor 6 is in agreement with the stoichiometry of carbon atoms in GLC and CO2). We can therefore define the fraction of GLC consumption through the PPP by re-scaling by :

From this it is also immediately clear that, in general conditions for CO2, F varies between (corresponding to ) and (corresponding to maximal ). Using the typical concentration values for the above metabolites in the blood, namely M, M, M and M one has and so that . Despite the roughness of the network reconstruction we employed, the bounds just obtained are in remarkable agreement with experimental evidence. On the other hand, using the empirical estimates M/s and M/s we find .

Discussion

In this paper we have derived a variational principle for the NESS of chemical reaction networks, showing that, for time scales over which chemical potentials can be considered constant, stationary non-equilibrium fluxes minimize the function H, see (13), in which stoichiometry, intakes, outtakes and concentrations appear as parameters. The cost function is closely related to Hopfield models of neural networks. This allows to rephrase such networks as systems of reactions interacting via Hebb-like rules and subject to an external forcing provided by the boundary fluxes. Furthermore H bears two simple physical interpretations. First, minimizing it amounts to finding the flux organizations that keep the overall “waste” of chemical species to a minimum, given the boundary conditions to be satisfied. Second, it equals (modulo a constant) the time derivative of the entropy production per volume: NESS thus correspond to flux configurations such that the entropy production decreases in time at the smallest allowed rate, in line with the ideas exposed in [30]. We have investigated the implications of such a picture in toy (random) chemical networks – where the dependence of the solutions on the network’s structural parameters can be fully explored, revealing the existence of a phase transition between flux balanced states with and states where unbalances emerge – and in a small real biochemical network, namely the metabolic network of the human red blood cell, a system that is possibly the closest to the theoretical situation described here.

Making contact with stoichiometric models of metabolic networks is relatively straightforward as long as boundary fluxes are taken to be fixed and one does not include additional objective functions that NESS are required to maximize. For instance, the maximization of biomass flux (a frequent optimization criterion for bacterial metabolism [33]) provides a source of entropy production even if one focuses on states with (as in Flux-Balance-Analysis, FBA). On the other hand the theory developed here suggests that the set of local constraints that describes NESS is provided by the minimization of H rather than by taking or a priori.

An important property that NESS should possess is thermodynamic feasibility, i.e. they should not contain infeasible loops [34], [35]. Inspired by the fluctuation theorem [8], we have proposed here a simple dynamical rule (see (20)) that ensures that the NESS obtained as the minima of H are indeed void of cycles. In general (i.e. when straightforward minimization of H is carried out), it is possible to get rid of infeasible cycles by complementing the variational problem described here with the minimization of the square norm of the flux vector. This is a consequence of the Gordan theorem of alternatives [36]: assuming that the matrix has full rank, only one of the following systems has a non-trivial solution: (a) (for real); (b) (for ). For a reaction network with stoichiometric coefficients and fluxes , defining one sees that system (a) expresses the condition of thermodynamic stability , whereas system (b) defines thermodynamically infeasible cycles, so that either a flux configuration is thermodynamically feasible or it contains at least one cycle. Now let us consider a steady-state flux configuration satisfying and let us assume it is thermodynamically infeasible, i.e. that system (b) has a solution . We can then construct a new flux configuration as , with a constant. Evidently, still . However defining , it is easily seen that choosing so that one gets . In other terms, starting from a thermodynamically infeasible steady-state flux configuration one can construct another steady-state flux configuration whose total flux is lower. In turn, has to be minimum when the flux configuration is thermodynamically feasible. It would be interesting to find a dynamical justification of this criterion. (Notice that flux minimization is an optimality principle frequently used for metabolic network modeling, see e.g. [37].).

Perhaps not too surprisingly, the formal aspects of this theory present many common traits with those developed over the past decade for the analysis of large games with heterogeneous interacting agents, specifically with Minority Games [24], [25]. For obvious reasons, such an analogy shouldn’t be stretched and we do not elaborate in detail on it here. In short, however, reactions (or, more properly, enzymes) in chemical networks can be thought to be involved in a competition for the use of a set of possibly limited resources (the different substrates). Whether a reaction can operate or not depends on how much substrate is available to it, i.e. on the overall substrate concentration and on how many other enzymes can bind the same substrates. Each reaction disposes of two ‘strategies’ (or ways to access the set of resources), corresponding to the vectors of its input and output metabolites in the forward and reverse direction, respectively. Such strategies are anti-correlated: substrates in the forward direction are products in the reverse, and vice-versa. Rules like (26) can then be read as ‘learning’ processes through which enzymes try to anticipate at each time step whether the substrates needed for the forward or reverse processes are most likely to be available, in order for it to operate. This parallel provides an elementary quantitative flavor to the idea that enzymes in biochemical reaction networks compete for the substrates. In absence of boundary fluxes ( for each i in (7)), the situation is completely equivalent to the Minority Game with anti-correlated strategies studied in [38]. In this case, the dynamics converges to and for each i, i.e. the system asymptotically reaches chemical equilibrium. As it should be, NESS are induced by non-zero boundary fluxes, i.e. by non-zero fields . This additional term turns out to be the main difference between standard anti-correlated Minority Games and the systems discussed here.

Besides a possible theoretical interest in deepening the analogy just described (e.g. by extending the dynamical approaches employed for the analysis of multi-agent systems [39] to models of chemical reaction networks), it will be interesting to see how well the variational principle (14) describes flux states in real biochemical networks.

Metabolites (abbreviation, full names, estimated concentrations) appearing in the reduced model of hRBC metabolism.

(PDF)

Click here for additional data file.

Reactions (abbreviation, enzyme name, formula) appearing in the reduced model of hRBC metabolism.

(PDF)

Click here for additional data file.