Inferring metabolic phenotypes from the exometabolome through a thermodynamic variational principle

Networks of biochemical reactions, like cellular metabolic networks, are kept in non-equilibrium steady states by the exchange fluxes connecting them to the environment. In most cases, feasible flux configurations can be derived from minimal mass-balance assumptions upon prescribing in- and out-take fluxes. Here we consider the problem of inferring intracellular flux patterns from extracellular metabolite levels. Resorting to a thermodynamic out of equilibrium variational principle to describe the network at steady state, we show that the switch from fermentative to oxidative phenotypes in cells can be characterized in terms of the glucose, lactate, oxygen and carbon dioxide concentrations. Results obtained for an exactly solvable toy model are fully recovered for a large scale reconstruction of human catabolism. Finally we argue that, in spite of the many approximations involved in the theory, available data for several human cell types are well described by the predicted phenotypic map of the problem.


I. INTRODUCTION
The wealth of biological data acquired in recent years via high-throughput techniques has lead to increasingly refined descriptions of a cell's content in terms of macromolecules and metabolites, as well as of the complex interactions between them that determine a cell's physiology. Such interaction patterns represent the networks underlying cellular organisation. Of course, various such networks (metabolic, signaling, protein-protein interaction, regulatory, etc.) operate on separate time-and space-scales, continuously cross-talking and reciprocally feeding into each other in single cells. In addition, because cells do not live in isolation, inter-cellular interaction networks can also be considered as representations of multicellular organization (e.g. tissues, organs, etc.). Physiologic functions emerge through the integration of collective biological interactions across different scales, from cellularlevel, to tissue and organ level. While this concept has always been at the center of biology, a formal and mathematical description of physiology in terms of the dynamics of network of networks has become possible only recently [1,2]. Obtaining a comprehensive picture of cellular activity from large-scale data is however still a major theoretical and computational challenge [3].
Mathematical network-based models represent a natural way to encode the inherent complexity of the interaction structure revealed by integrated biological data. Different approaches have been developed for a variety of cellular networks [4]. In many cases, the basic idea on which such schemes rely is that of optimization: assuming biological networks are optimized to perform a specific, context-dependent biological function (e.g. biomass or energy production in metabolic networks) one focuses on the network states that maximize the specific functionality. While such methods have met predictive and/or explanatory success [5], optimal states normally depend on a multitude of interaction parameters, which makes robustness an issue. Furthermore, identifying objective functions is not always straightforward, as testified e.g. by the variety of optimality criteria that can be used to describe, to different degrees, the same system [6,7]. In some cases, it is even hard to argue that something is being optimized at all. A different approach consists in placing the emphasis not on functional aspects but on the physical constraints under which such networks operate, and in trying to identify generic variational principles (similar to those derived for physical systems) by which the space of possible network states can be reduced to the 'physically relevant' ones. Clearly, states selected by physical variational rules will in general be unable to pin down a specific biological functionality. The picture they provide will however be inherently robust and, in certain cases (specifically whenever physical constraints set the relevant limits to the network's operation), one may hope to obtain biologically relevant insight.
One such case is possibly that of cellular metabolic networks. A cell's metabolism is, in essence, the complex network of enzyme-catalyzed reactions that processes nutrients to derive energy and molecular building blocks, while harvesting free energy from the environment and allocating it in the multiple tasks a cell has to accomplish. Metabolic network models have been extended to the scale of the whole genome and optimization-based schemes are routinely employed for their mathematical analysis [8,9]. Stationarity can be argued to be an appropriate assumption in order to investigate the productive capabilities of the network in terms of output metabolites or maximal flow of specific reactions (including biomass production). Therefore, upon defining suitable objective functions, one may resort to frameworks such as Flux Balance Analysis (FBA) for their study [5]. At the computational level, FBA is normally a linear programming (LP) problem, and it has been particularly successful in predicting the growth rate of microorganisms in batch cultures. In cases in which a clear objective function is lacking, as is typical in multicellular organisms, an unbiased sampling of the steady states can provide useful information on the organization of reaction pathways or on the design of experiments [10]. The practical feasibility of sampling algorithms on genome-scale systems is however still a concern.
It is therefore important to understand how physics, and arXiv:1410.6364v1 [q-bio.MN] 23 Oct 2014 thermodynamics in particular, constrains the solution space of metabolic networks. Recently, a simple and intuitive physical description of steady states based on a thermodynamic variational principle has been proposed, according to which reaction networks occurring in a given volume and exchanging compounds with a stationary environment (a 'bath') tend to minimize the rate of decay of entropy production [11]. The latter quantity can be written explicitly in terms of the stoichiometry of the reaction network, of intracellular reaction fluxes and of extracellular metabolite levels. Therefore, the rule allows in principle to explore the physically viable intracellular flux configurations upon changing the composition of the environment.
In this article we shall apply this framework to infer viable steady-state flux patterns in the catabolism of human metabolic networks, with the goal of clarifying the extent to which physics accounts for the switch in cellular energetic strategies from fermentative to oxidative phenotypes that is observed in many cell types [12][13][14][15]. We will show that, perhaps surprisingly, thermodynamic principles alone predict the crossover between different metabolic phenotypes using a small number of 'environmental' control parameters, such as the external glucose, lactate, oxygen and carbon dioxide levels. Indeed, despite the crude approximations on which the theory is based, experimental observations fall remarkably well within the derived scenario. These results ultimately suggest that a more thorough understanding of the physical and chemical constraints under which biological networks operate may provide us with much conceptual insight and possibly predictive power.

II. METHOD
Let S denote the stoichiometric matrix of a given metabolic network with N reactions and M metabolites, with S µi the stoichiometric coefficient of compound µ in reaction i. Upon neglecting the discrete nature of molecules, noise and spatial gradients, the mass balance equations for the concentration c µ of each compound µ may be written in terms of the reaction fluxes {f i } asċ The Gibbs energy of reaction i may in turn be decomposed in terms of the chemical potentials g µ of the different compounds, i.e.
where, for a well-mixed and diluted system, the chemical potentials at constant pressure and temperature are given by where R is the ideal gas constant and concentrations are assumed to be measured in units of a fixed reference level.
Differentiating (2) with respect to time and applying (3) and (1), one sees that Introducing the shorthands Equation (4) can be re-written aṡ where This elementary derivation suggests that the time evolution of Gibbs energies follows the gradients of the quadratic function H of the fluxes. In turn, steady flux states correspond to the minima of H. (See [11] for a more precise microscopic derivation.) Notice that H ≥ 0 by construction. It can furthermore be seen [11] that H corresponds to the rate of entropy decay, i.e. H = −S/R, with S the internal entropy of the system [16]. Note that, in terms of the concentrations, H takes the form It is convenient to distinguish the levels of intracellular compounds (c µ,int ) from those of extracellular ones (c µ,ext ), so that In turn, the variations of intracellular and extracellular concentrations are linked by the exchange fluxes u µ . If we single out the latter (assuming a positive sign for fluxes entering the cell), we get where is the ratio of intracellular and extracellular volumes: A trivial solution to the H-minimization problem is obtained by taking vanishing fluxes f i and uptakes u µ , leading to H = 0. We shall focus on solutions carrying non-vanishing fluxes, corresponding to non-equilibrium steady states. Since it is reasonable to think that the volume ratio will typically be small, the second term in (10) dominates H. In metabolic networks, however, the constraints provided with models are normally compatible with internal homeostasis, implying i S µi f i + u µ = 0. Therefore the H is minimized by minimizing the first (exchange) term alone. In summary, for given extracellular levels c µ,ext , the variational principle takes the form of the optimization problem .
(11) For sakes of simplicity, we shall henceforth re-define In short, the above problem amounts to the minimization of a positive definite quadratic convex function in a convex polytope defined by the conditions (11). Finding a solution requires in general polynomial time [17]. Sampling all solutions uniformly, on the other hand, can either be done exactly upon knowing the vertices of the polytope, or by stochastic methods (e.g. Monte Carlo). The latter procedure is also polynomial when a Hit-and-Run Markov chain is employed [18][19][20][21][22].
In the following, we shall explicitly solve the above problem for specific metabolic networks and compare its solutions both with empirical data and with solutions derived from another variational principle that is widely employed to describe out-of-equilibrium systems.

A. Minimal model for ATP production
As a first example, we consider a minimal model for the production of ATP from either glucose (through oxidative or fermentative metabolic pathways) or lactate (via oxidative pathways only), with the goal of inferring how a cell employs the different energy-producing pathways as a function of the extracellular levels of nutrients and waste products. Fig. 1 displays schematically the reactions of our minimal model. Each glucose molecule enters the cell and is transformed to 2 pyruvate molecules thereby converting two ADPs to ATPs (note that the stoichiometry is not depicted in the figure). Each pyruvate can either be used to produce further 15 ATPs, with the concomitant consumption of 3O 2 and production of 3CO 2 (oxidative pathway), or be transformed to lactate (fermentative pathway). The latter reaction is reversible so that lactate can be both expelled and intaken by the cell; in contrast, O 2 can only be intaken and CO 2 can only be expelled. We set the stoichiometry of this simple model to match the one of the complete network presented below; in general, the precise stoichiometry of ATP production via oxidation varies across species, with a yield in the range of 9 to 18 ATPs obtained per pyruvate molecule [23]. Denoting by u µ the exchange reaction of metabolite µ (keeping in mind that u µ > 0 for intakes and u µ < 0 for outtakes), by f ox the flux trough the oxidative pathway (labeled 5 in Fig. 1), and by f atp the ATP production flux (the sum of fluxes through reactions 2 and 5), the homeostatic intracellular steady state is defined by the equations where the subscripts glc, lac, co2, and o2 stand for glucose, lactate, carbon dioxide, and molecular oxygen, respectively, and where we have used the fact that, by homeostasis, the flux through glycolysis equals u glc . Since we have four extracellular species, H as defined in (12) is given by Fixing the ATP production flux conventionally to f atp = 1 (this serves no specific purpose except fixing a scale for fluxes) and using the steady-state conditions (13), (14) can be re-cast as In turn, the value of f ox that minimizes H can be easily obtained by differentiating the above equation. The minimum of H is obtained when Let us analyze how the emerging scenario changes with the extracellular levels. From the steady state conditions (13), if there is no exchange of lactate, i.e. if u lac = 0, then glucose, which in this case is the only carbon source, is completely oxidized and f ox = 1/16. Substituting this value in (16), we obtain that the glucose concentration corresponding to zero lactate exchange, which we denote as c glc , verifies Note that c glc is independent of the external lactate concentration c lac . For the physiological levels in the blood plasma (c co2 30 mmol, c o2 5 mmol), one finds c glc 1.8 mmol.
On the other hand, combining (13) with (16), and using expression (17), one finds that u lac ∝ (c glc − c glc ). This suggests that (17) defines a threshold separating, for any given levels of oxygen and carbon dioxide, different metabolic phenotypes. In particular, if c glc > c glc one has u lac < 0 (there is a net lactate secretion), while u lac > 0 for c glc < c glc (corresponding to a net lactate intake). We start by considering in detail the case c glc > c glc . Here, ATP is produced using glucose exclusively, which can be channeled to both the oxidative and the fermentative pathway. By recalling that one glucose produces two pyruvates, the fraction of oxidized glucose can be written as O = f ox /(2u glc ) ≤ 1, which, in the steady state described by (13), becomes simply Substituting expression (16) for f ox , one obtains the percentage of oxidized glucose as a function of the external metabolite levels, namely where In the blood plasma, one has typically c glc 5 mmol and c lac 1 mmol, so that O 0.92. In other words, oxidation is the dominant energy-producing strategy for cells in contact with the blood.
If instead the glucose level is below threshold (c glc < c glc ), then there is a net lactate intake in addition to glucose and no fermentation is possible, as it would imply lactate secretion. In this case, the relevant quantity to consider is the fraction of carbons that are intaken as glucose, or (recalling that one glucose molecule is converted to two pyruvate molecules) G = 2u glc /(2u glc + u lac ) ≤ 1, which at steady state is To summarize: if c glc > c glc (or c glc /c glc > 1), the only source of carbon is glucose (G = 1) and the fraction of carbons that are oxidatively processed is given by O. If instead c glc < c glc (or c glc /c glc < 1), oxidative carbon processing is the only possibility and O = 1. The fraction of carbons intaken as glucose is concomitantly given by G. Noticing that G = 1/O, one can put different regimes together by defining a new function L, representing the normalized number of carbon atoms exchanged as lactate, being positive for outtakes and negative for intakes, which is given by Note that L ∈ [−1, 1]. This allows us to organize possible metabolic phenotypes in a single diagram. The contour plot of L in the plane (c glc /c glc , c lac /c glc ) is displayed in Fig. 2, where four regions can be distinguished: a mainly fermentative regime with lactate outtake (red), a mainly oxidative regime with glucose as the exclusive carbon source (white, c glc /c glc > 1), a purely oxidative regime with glucose as the main carbon source (white, c glc /c glc < 1), and a purely oxidative regime with lactate as the main carbon source (green).
Within the H-minimization framework, the diagram maps out the internal metabolic states of a cell as a function of the extracellular concentrations of glucose and lactate. The levels of oxygen and carbon dioxide are implicitly included via the parameter c glc and serve as scaling factors. In essence, Fig. 2 is the complete solution to the inference problem of determining the main carbon source and the pattern of pathway usage by a cell when the extracellular concentration of exchanged metabolites is known. We shall see that, despite its crudeness, the scheme just described captures the key features of cellular energetic strategies. Indeed, the same picture will emerge from a much more detailed model of cell metabolism.

B. H-minimization versus Minimum entropy production
It is instructive to compare these results with those obtained by minimizing entropy production, a well-known variational principle for biochemical systems, described e.g. in [16]. The general expression for the entropy production is given by [16] where T is the temperature and g µ is the chemical potential of metabolite µ.  (17), separates the plane in two zones depending on whether cells intake (c glc /c * glc < 1, L < 0) or secrete (c glc /c * glc > 1, L > 0) lactate. The intensity of the colors represents the relative value of the lactate flux. For small values of c lac /c * glc the lactate flux is always negligible, while for c lac /c * glc 10 the amount of lactate exchange is large and, by crossing c glc /c * glc = 1, switches rapidly from a large lactate intake to a large lactate secretion. L is defined in (22) and represents the fraction of carbon atoms exchanged as lactate. section, it takes the form In addition, the constraints u glc ≥ 0 (glucose is entering the cell) and f ox ≥ 0 lead, via (13), to 0 ≤ f ox ≤ 1/15. Entropy production is easily seen to be minimized by two states only, depending on the sign of the coefficient α ox . In particular, it is minimized by taking f ox = 0 (complete fermentation) if α ox > 0, or f ox = 1/15 (corresponding to no glucose uptake, and complete oxidation of lactate) if α ox < 0. Notice that, within this simple model, the minimization of entropy production predicts no intermediate (or mixed) ATP producing strategy, and in particular no glucose oxidation. Therefore, according to this variational principle, fluxes will saturate their lower or upper bounds depending on the sign of certain linear functions of the chemical potentials. Such a scenario inevitably leads to "extreme" regimes separated by sharp switches between them, which seems unlikely to agree with biological reality.

C. Large-scale model of ATP production
A realistic model of energy production by cells can be obtained by including the backbone of four ubiquitous pathways leading to ATP production, namely Glycolysis, Pentose Phosphate Pathway, Citric Acid Cycle, and Oxidative Phosphorylation. We have built such a network, extracting relevant reactions from the human reactome Recon-1 [24] (the complete list of reactions is reported in Table I). Altogether, the network comprises 49 chemical species and 45 reactions (18 of which are irreversible). The reactions include the shuttling of six metabolites: molecular oxygen, carbon dioxide, water, hydrogen, lactate, and glucose. Our goal is to analyze its feasible steady states according to the H-minimization principle, along the lines followed for the reduced toy model discussed above.
To apply the variational principle in this case, we again impose steady state conditions for intracellular fluxes and single out the exchange reactions. Intracellular homeostasis is described by These equations generate a large number of linear dependencies among fluxes, which can be resolved explicitly by transforming the stoichiometric matrix to its Reduced Row Echelon Form (RREF) through Gaussian elimination [25]. It turns out to be possible to represent internal fluxes in terms of five degrees of freedom only, which we label as (x 1 , x 2 , x 3 , x, u).
The RREF directly provides an expression of the different fluxes in terms of these parameters. Such expressions, emphasizing the biological meaning of the five independent degrees of freedom, are shown in the last column of Table I. In short, x 1 describes the so-called Rapoport-Luebering shunt (reactions catalyzed by DPGM and DPGase); x 2 represents the ATP consumption, which will be fixed to 1 to match the corresponding production flux (see Sec. III A); the flux through the Citric Acid Cycle is described by x 3 ; u corresponds to the glucose uptake; and, finally, x represents the flux through the superoxyde dismutation that reduces O − 2 An expression for H can be obtained directly through the expressions of exchange reactions, which are also defined by the steady state condition, as functions of the five independent degrees of freedom. The latter can be read off from Table I. In particular, the expressions for the exchange fluxes of lactate, glucose, carbon dioxide, and oxygen are given by Note that they depend on the two parameters (x, u) exclusively. We neglect the exchanges of water and hydrogen ions assuming to work in biochemical standard conditions, where the water level is taken to be large and hydrogen is buffered: in other words, c h2o 1 andċ h 0. As a consequence, the corresponding terms in H are negligible. Notice also that corresponding to the mass balance of carbon atoms. To characterize the domain of values of (x 1 , x 2 , x 3 , x, u) where the minimum of H should be sought, one has to consider how the reversibility assignments encoded in the reaction network transfer to the independent variables. Indeed, while homeostasis defines linear dependencies among fluxes (i.e., linear equalities), the 18 irreversibility constraints (see Table I, with a positive (resp. negative) flux conventionally taken for the forward (resp. reverse) direction) define linear inequalities for the five remaining degrees of freedom. From Table I it can be easily recognized that all five degrees of freedom must be non-negative. A direct analysis of the redundancies of the constraints on the irreversible fluxes in Table I shows that the non-redundant constraints are the ones determined by the enzymes ACYP, ICIDHy, PFK, PGK, and PGL. In particular, one finds that the independent variables are cross-linked by the conditions where we used the fact that x 2 = 1. Inequalities (28) together with the non-negativity constraints determine the convex domain where the minimum of H must be found. Now, using (26), H is given by Note that H is a quadratic function with absolute minimum in the origin. The latter point however lies outside the convex domain defined by (28). Therefore the feasible solution of the H-minimization problem will lie on an edge of the domain.
Simple algebraic calculations reveal that H is minimized on the edge defined by The condition 250x + 2u = 1 reduces H to a function of a single variable that, for convenience, we redefine to be s = 250x with s ∈ [0, 1]. With this substitution, we are left with the simple problem of minimizing the function where a = 1 + 101 1500 . The minimum of H turns out to be attained when which is a single point on the boundary of the convex domain. As for the minimal, we are interested in describing the emerging metabolic phenotypes in terms of (a) the pattern of pathway utilization, and (b) the substrates from which ATP is preferentially produced. The emerging scenario is indeed very similar to that derived in the simpler case. In specific, • s → 0 for c o2 → 0 or c co2 → 0, corresponding to glucose intake and complete fermentation with lactate outtake; • s → 1/a for c lac → 0, corresponding to glucose intake and complete oxidation; • s → 1 for c glc → 0, corresponding to lactate intake and complete oxidation.  In other words, as s is changed in [0, 1/a] one passes continuously from fermentative (s = 0) to oxidative (s = 1/a) phenotypes, with lactate outtake, whereas for s ∈ [1/a, 1] one has complete oxidation without fermentation, and the preferred fuel switches continuously from glucose (s = 1/a) to lactate (s = 1). The value s = 1/a defines the curve 240 101 through which a threshold value c glc for the glucose level c glc is obtained that, in parallel with the minimal model, can be conveniently used to separate different metabolic phenotypes.
In physiological conditions for the blood (c o2 5 mmol and c co2 30 mmol), one has c glc 1.8 mmol. Below threshold, one finds partial lactate intake, whereas above threshold partial lactate outtake (fermentation) is observed. In this realistic model, the fraction of oxidized glucose is given by For the average value of glucose and lactate in blood (c glc 5 mmol and c lac 1 mmol), (34) predicts a small lactate outtake, the phenotype being almost completely oxidative Below the line the phenotype is oxidative with partial intake of lactate, above the line the phenotype is partially fermentative with lactate outtake. (b) Fraction of lactate uptake, with respect to glucose (L < 0, green), or with respect to carbon dioxide (L > 0, red) in the plane (c glc /c * glc , c lac /c * glc ). L is defined in (35) and represents the fraction of carbon atoms exchanged as lactate.
(the fraction of glucose oxidized is O 0.92). In Fig. 3  (top panel) we display the critical curve c glc in the plane (c glc /c co2 , c glc /c o2 ).
To make contact with empirical results, we have considered the experiments discussed in [26][27][28], where data are given for glucose and lactate levels for the Blood Plasma (BP) and for the Extra Cellular Fluid (ECF) in brain tissues, both in human and in rat. By plotting these data on the (c glc /c co2 , c glc /c o2 ) plane, we can in principle assess whether the metabolisms of brain cells in such environments is fermentative or oxidative and what is their preferred energy source. However, lacking data for carbon dioxide and oxygen levels, in order to carry out a detailed comparison we have employed the physiological values c o2 5 mmol and c co2 30 mmol. Strikingly, experimental points distribute so that H-minimization predicts a lactate intake for cells in the ECF and a lactate secretion for cells in the BP, in agreement with the conclusions drawn in [26]. Since Fig. 3 (a) does not describe the dependence on lactate levels, it cannot provide the relative value of lactate uptake, with respect to glucose intake or carbon dioxide outtake. This information can however be retrieved by studying how experimental points distribute in the (c glc /c glc , c lac /c glc ) plane with respect to the contour plot of the quantity L (the normalized number of carbon atoms exchanged as lactate), which can be defined in analogy to (22) as The resulting phenotypic map is shown in Fig. 3 (b), where  [27] and Rats [28]. The first column reports the source of the experimental data. The second and third columns report the measured levels of glucose and lactate, respectively. Finally, the last two columns show our estimate for the direction and of lactate and for the absolute value of |L|, respectively. L is defined in (35) and represents the fraction of carbon atoms exchanged as lactate.
the emerging scenario is that of a partial fermentation in the BP and partial lactate intake in the ECF. Values for the experimental data measured for both ECF and BP are reported in Table II.

IV. DISCUSSION
Metabolic flux analysis for cell-autonomous systems and for populations [29] is a rich, fruitful and expanding field of research, displaying multiple connections between theory and experiments. For the greatest part, modeling schemes rely on the possibility to identify objective functions by which the complexity of the space of possible solutions can be reduced, focusing on optimal flux patterns exclusively. In this way, given a set of uptake fluxes describing a cell's exchanges with the environment, the physiologically relevant states of the cell can be mapped onto a small set of flux configurations (possibly reduced to a single point), which can be studied in detail by different computational techniques. Sampling methods allow in principle to explore the solution space beyond optimality, to characterize, e.g., sub-optimal states, correlations, etc, and improving such techniques is one of the current frontiers of the field. At the other end of the modeling spectrum, one may be interested in characterizing how physics constrains the solution space by means of variational principles that characterize steady states in terms of the minima of specific functionals that are physically-(rather than biologically) motivated. The thermodynamics of reaction networks, a subject that goes back to at least [30] and it has recently gained attention from a stochastic perspective [31,32] suggests that steady states should be characterized in terms of their entropy production. Since such principles should hold for cellular systems as well [16,[33][34][35], it is interesting to study what type of information can be obtained about a cell's metabolism by applying these ideas.
In this work, we have studied the problem of inferring the metabolic phenotype from the level of metabolites in the extracellular space (i.e. from the exometabolome) using a thermodynamic variational principle for the steady states of a chemical reaction network. The variational principle, which holds for slowly varying chemical potentials, amounts to the minimization of the rate of decay of entropy production. From a conceptual viewpoint, it merely allows us to select intracellular flux states that are compatible with the observed extracellular concentrations. From a technical viewpoint, it requires the minimization of a semi-positive definite quadratic function of the fluxes in the space of feasible steady states, and can be carried out in polynomial time.
We have applied it to the catabolic core of a detailed genome-scale reconstruction of the human metabolic reactome with the goal of characterizing the conditions under which cells switch from a fermentative to an oxidative phenotype as a function of the external levels of key environmental indicators like glucose, lactate, oxygen and carbon dioxide. Our results indicate that cells transition from one phenotype to the other in a continuous, modulated way, and that mixed phenotypes are possible. This is in line with empirical knowledge [36,37] and at odds with the scenario predicted by another widely used variational principle (that of minimal entropy production), which predicts sharp transitions [16] rather than smooth cross-overs. The scenario we obtain is also recovered in an exactly solvable toy model that only retains the main features of the key energy producing pathways. Quite remarkably given the crudeness of the variational principle we employ, upon inferring from experimental values of glucose and lactate levels in the brain we find very moderate levels of lactate exchange, in specific a small intake in the extra-cellular fluids and a small outtake in the plasma. Experimental evidence compares well with the predicted "phase structure", suggesting that indeed fundamental physical considerations might suffice to explain at least part of the evidence on cellular energy production strategies. One of the limitations of the approach discussed here is that, in principle, knowledge of intracellular substrate levels (besides extracellular concentrations) is required to solve the full-fledged variational problem. In this work, we have circumvented this difficulty by assuming that a mass-balanced flux pattern for the intracellular state. This assumption is justified for the type of systems we consider, but cannot be expected to hold generically. It would be interesting to see how strongly solutions for large-or genome-scale networks depend on the particular internal metabolite pools. In turn, characterizing robustness to fluctuations in metabolite levels may highlight the presence of bottlenecks in the reaction network. On a more abstract level, it would also be important to generalize these considerations to a dynamical setting, e.g., characterizing trajectories (as opposed to steady states) in terms of physical or thermodynamical variational principles.