BACKGROUND
when studying metabolism, it is important to elucidate how fluxes are distributed among the different pathways of the metabolic network with respect to the available quantitative information about the system behavior.
several methods can be used to address this issue.
the first approach consists of building a mechanistic description of transformations and identifying the regulations involved in the system.
continuous dynamical models are often used for this purpose, especially when time-series responses to different treatments are available to infer the dynamics of the network.
static approaches such as petri net can also identify qualitative distributions of fluxes in a metabolic network  <cit> , and their stochastic extensions can even take into account stoichiometric and kinetic information  <cit> .
with a complementary approach, one can study a system at the functional level, based on the study of fluxes at steady states.
this is the purpose of the flux balance analysis  framework, which has evolved considerably over the past decades  <cit> .
with an fba, the identification of external regulations is not necessary because it is assumed that the global behavior of the system can be modeled by optimizing linear combinations of selected fluxes .
roughly, the methods developed in this field aim to explore a convex space of plausible flux distributions and to study the extreme flux distributions obtained when optimizing a linear objective function.
it allows checking whether an extreme flux distribution is consistent with experimental data and to predict new experimental observations  <cit> .
nonetheless, the consistency of the solutions obtained by fba depends on the quality of the constraints integrated in the model.
to overcome this limitation, several extensions have been proposed in the literature.
these extensions can be broken down into two parts.
the first incorporates additional biological knowledge, such as reaction thermodynamics  <cit>  or multioptimization  <cit> .
the second is based on the use of fba to globally analyze large-scale metabolic networks.
for example, in flux variability analysis  <cit> , the minimum and the maximum flux for each reaction in the network are computed under some optimal conditions.
in this paper, our main purpose is to extend this framework to study the variability of a metabolic network at the level of efficiencies instead of fluxes.
indeed, when studying a metabolism at the organ level, a major challenge consists in comparing the efficiency, or yield rates, of two metabolic situations, i.e.
the response to various input patterns.
a typical example of such studies are those concerning animal nutrition, which aim to predict the quality and quantity of animal production  in response to breeding factors.
in this field, the energy and protein conversion efficiencies are derived from the study of the flux distribution of input nutrients between the different branches of the metabolic network  <cit> .
more generally, although its definition depends on the field of application, the concept of efficiency is often linked to energy, mass growth and protein conversion  <cit> .
however, the computing of efficiencies is prone to difficulties and errors, since it requires computing the quantities precursors sets which are required to explain the measured composition of outputs, although some of the internal products may be recycled within cycles  <cit> .
with this goal, we defined the allocation of an input towards an output  to be the proportion of a matter component  in a given input flux that is recovered in the selected output flux.
our definition is based on the choice of a material component, such as carbon or nitrogen, which allows a comparison of the contributions of input metabolites to the composition of output products.
it can be seen as a yield rate, which is uniform among all the outputs of the system, allowing a precise understanding of the precursors used.
as a first methodological contribution, we prove that aio can be uniquely described and computed, even in the case when there are metabolic system cycles, by a matrix whose coefficients are nonlinear functions of the flux variables.
introducing nonlinearity in the definition of aio cannot be avoided because of the presence of these cycles.
studying the variability of aio within a complete space of plausible flux distributions requires the solving of nonlinear optimization problems which are underdetermined in tangible applications.
as a second methodological contribution, we have proposed efficient algorithms to compute lower and upper bounds for aios over the family of flux distributions which are compatible with both the system’s stoichiometry and the experimental datasets, regardless of the choice of a flux distribution for the internal branches of the network.
an important aspect is that when the metabolic network is provided with input-output data, the complete space of plausible distributions appears to have a relatively small size, and can therefore be studied with our method.
our framework is depicted in figure  <dig>
our main example of application is related to milk production.
in this context, several models have been introduced, in relation with the aforementioned classification of models.
one class of small-size dynamical mechanistic models predicts the blood flow and input nutrients of the metabolic system   <cit>  or, alternatively, the nutrients produced by the metabolic system   <cit> .
another class of models predicts the distribution  of fluxes over the pathways from the input and output nutrients of the system  <cit> .
in the latter family, both dynamic and static approaches exist.
indeed, numerical models based on mass-action equations were initially proposed to describe the fluxes to and from individual metabolites, with different levels of description  <cit> .
to that goal, optimizers were used to determine a reasonable set of parameter estimates for the dynamical model of the system.
although such a set of parameter estimates is not unique.
with a totally different approach, in another study, the main biochemical pathways in the different metabolisms were integrated in a generic stoichiometric model called the metabolic spreadsheet <cit> .
its structure was based on a restricted number of intermediary metabolites called carbon-chain pivots, which correspond to important cross-over points between metabolic pathways.
this allowed the computation of the flux rates for all reactions, constrained by a general rule stating that there is no accumulation of intermediary metabolite of cofactors.
a study of manual calculations performed with the metabolic spreadsheet showed that it works, in practice, by maximizing an objective function  in a convex solution space  <cit> .
therefore, this model can be considered as an application of the flux balance analysis framework  <cit> .
nonetheless, in this field of study, there has been little discussion on the impact of the choice of a single model among several  reasonable models  <cit> .
to investigate this issue in a more automatic way, we compared the aforementioned conventional models to the convex space of plausible flux distributions associated with steady states of a model of mammary gland metabolism, as computed in fba.
we checked the consistency of extreme flux distributions of nutrients with experimental data of the mammary gland and milk production, including the contribution of nutrient input-output and isotope balance studies  <cit> .
based on our aio computation framework, our analysis highlighted that the metabolic behavior of the mammary gland cannot be modeled by maximizing atp production or by optimizing a linear combination of flux variables of the model of the mammary gland.
in other words, although an infinite number of flux distributions are compatible with the data, none are extreme within the space of feasible models.
selecting any of these nonoptimal flux distributions to predict the system behavior appears difficult without additional experimentation.
to gain better understanding of the system response regardless of the choice of a flux distribution for the internal branches of the network, we applied our method to estimate the variability of aio coefficients in our model and compared the effects of two different diets on mammary gland metabolism.
our results suggest that the bounds of aio are sufficient to distinguish the effects of different nutritional treatments without selecting a flux distribution for the internal reactions of the metabolic network by any method - optimization of a linear combination of fluxes or a residual score.
overall, the complete study suggests considerable flexibility in mammary gland metabolism.
it provides a view of the functioning of the system although its internal processes still cannot be clarified because of limitations on experimentation on large animals such as ruminants.
RESULTS
we first investigated the set of flux distributions that are compatible with the stoichiometry of our mammary gland model , without taking datasets into account.
the model exhibited a large variability, since thousands of extreme pathways could be identified.
we then successively computed the set of flux distributions compatible with the model and the real datasets of lactation metabolism in dairy cows.
the datasets are given in table  <dig>  they include a control diet , a diet related to an increased protein supply through casein infusion into the duodenum , and a complementary dataset  previously used in a mechanistic model of the mammary metabolism  <cit> .
for the three datasets, the set of plausible flux distributions – solutions to a linear system detailed in eq.
– was an unbounded convex cone of dimension  <dig>  in all cases, the solution spaces shared the same set of five independent variables.
this suggests that there are five independent levels of variability within the system: peptide hydrolysis , nadph oxidation , oaa →pyr , oaa →g3p  and g3p →g6p .
therefore, by including this dataset the model became a fairly small and constrained network.
nevertheless, it was not uniquely determined since there were still several degrees of freedom.
fatty acid output
amino acids balance
fa primer from acetate
data are renormalized in mmol/h/half udder.
control diet  <cit>   higher protein supply by casein infusion in the duodenum  <cit>   generic dataset  <cit> .
this table is used to parameterize input and output vectors v
i
and v
o
together with additional biological linear constraints on some reaction fluxes.
1input i.e.
taken up by the stoichiometric system considered .
2
β
-hydroxybutyrate.
3total triglycerides secreted in milk considering that milk fat was composed of 100% triglycerides and that all the triglycerides were secreted in milk fat  <cit> .
4output i.e.
leaving the system .
5all fatty acids synthesized within the mammary gland i.e.
all c <dig> to c <dig> and 50% of c <dig>  <cit> .
6the balance between amino acid net uptake and amino acid net output in milk protein is calculated with established rules  <cit> .
if the balance is positive it corresponds to an amino acid input that will be catabolized.
if this balance is negative, the values corresponded to an amino acid output .
7serine output corresponded to serine synthesized minus serine utilized in other pathways i.e.
serine required in addition to ser uptake to synthesize milk protein.
8peptide output: number of peptide links required to synthesize the proteins exported out of the system .
9hanigan,  <dig>  <cit> .
10fatty acid primers were synthesized for 50% from acetate and for 50% from bhba except c <dig> fa primer which was supposed to be synthetized only from bhba  <cit> .
11set at zero because their inputs or outputs are set at zero .
investigating the relevance of the optimization strategies for mammary metabolism
the balance between the atp generated by the system and the atp used by the system  was computed for the three datasets ,  and .
the results are detailed in table  <dig>  in this table, the fba approach based on optimization of the atp balance is compared to the natural functioning of the system.
criteria of selection of a solution
three natural assumptions are considered to model the mammary gland behavior: removal of all cycles, optimization of atp production and study of the equilibria of a dedicated ode-based model.
all models exhibit considerable variability in their atp balance , which contradicts the assumption about the behavior of this organ.
moreover, quantitatively, the computed atp balances are much higher than recent measurements.
this suggests that atp maximization cannot be considered as a natural objective function to model cow mammary behavior.
first, we considered the manual computation of fluxes with a tool named “metabolic spreadsheet”  <cit> .
we applied the rules of no accumulation on any intermediary metabolites and cofactors with a utilization rejecting all the cycles.
as expected, since some cycles are atp-consuming, both flux distributions were equivalent in both approaches .
the atp balances were respectively  <dig>   <dig> and  <dig> mmol/h/half udder  in the ,  and  datasets.
according to this model, the atp generated is estimated to be  <dig> mmol/h/half udder  while  <dig> mmol/h/half udder  were estimated to be used for milk component synthesis  <cit> .
therefore, the atp-balance obtained for this model was slightly smaller than the optimal one obtained for our mammary gland model.
independently of the approach considered and the model or dataset at hand, the remaining atp , after use for milk synthesis, appeared to be rather high and nonxconstant.
indeed, as shown in  <cit> , an atp balance of  <dig> mmol/h/half udder  can be expected to be used for other functions not accounted for in the models, such as maintaining membrane potential and synthesizing nucleic acids.
the numerical values obtained here were far from this estimated balance, suggesting that each of these three models have a nonrelevant atp balance.
in addition, both in the ode model and the atp-optimization approach, peptide hydrolysis was obtained at zero, implying an absence of any protein turnover.
this contradicts all the observations about this pathway: considerable use of this pathway has been evidenced in several publications, although the peptide hydrolysis rates differed significantly depending on the technique used for the measurements: peptide hydrolysis  spans from  <dig> , <dig>  to  <dig>  of peptide synthesis  in ,  and , respectively  <cit> .
overall, we concluded with this analysis that the energy-based optimization function may not allow an appropriate simulation of mammary gland metabolism.
this was expected, considering that the system is studied at the complete organ level, involving competing processes which can rarely be modeled with a single linear objective function.
exploring all extreme flux distributions in a refined simplex
in order to study the variability within the space of plausible flux distributions and to identify alternative relevant optimization strategies, an additional constraint was placed on the atp balance of the system.
we considered that an atp-balance of  <dig> mmol/h/half udder  was a relevant measure for this study  <cit> , although we checked that all the results discussed below were still valid when introducing a 10% tolerance on this estimation of atp-balance.
providing a bound for atp yielded a new constraint on several variables of the system.
for the two datasets  and , the flux g3p →g6p  was no longer an independent variable in the system.
the simplex, which was previously unbounded, appeared to be bounded with four independent variables: oaa →pyr , oaa →g3p , nadph oxidation , peptide hydrolysis .
applying flux variability tools  <cit> , it appeared that the minimum of each of the fluxes v <dig> v <dig> v <dig> v <dig> was equal to zero for all treatments.
the maxima of the fluxes  were  for the  treatment and  for the  treatment this suggests that the  treatment generates a more flexible space of plausible flux distributions than the  treatment.
we computed all extreme vertices of the simplex of plausible flux distributions for the two treatments , .
the simplex structure of the polyhedron implies that the optimum of any linear combination of metabolic fluxes involved in the model is either uniquely attained for one of these extreme points or attained by all flux distributions positioned on a face of the simplex.
to gain insight on the pathways involved in the variability of our model, we also computed the linear combination of fluxes optimized for one of the extreme flux distributions.
eight extreme vertices were found for the  and  datasets.
they were obtained when optimizing the same set of linear functionsa.
in table  <dig>  the optimal flux distributions for ,  are classified according to the activation or inactivation of several pathways in the model.
it is worth noting that, according to this table, the two sets of plausible flux distributions associated with  and  have the same topological structure.
indeed, the eight extremal behaviors for  and  have clear combinatorics: first, nadph oxidation is set either at zero or at its maximal value.
then, a single flux within oaa →pyr , oaa →g3p , g3p →g6p , peptide hydrolysis  is strongly activated whereas the three other remaining fluxes are blocked.
r
19
r
15
r
14
r
8
r
64
r
63
r
13
→
the qualitative properties of all vertices are shared in  and  treatments.
both correspond to a simplex with height vertices.
so are six of each of the  and  simplex vertices.
h and g vertices, in the  and  treatments, are plausible with respect to the literature.
as discussed in a previous paragraph, flux distributions with no peptide hydrolysis cannot be considered as relevant  <cit> .
this suggests the extreme flux distributions a-f in  and  are not biologically relevant.
on the contrary, two extreme flux distributions  are consistent with the stoichiometry of the system for the  and  treatments.
these two flux distributions consist in optimizing peptide synthesis and hydrolysis after nadph oxidation.
they corresponded to a ratio between peptide hydrolysis and synthesis v64v <dig> of  <dig>  in the  treatment and of  <dig>  in the  treatment.
these ratios are equal or lower to the maximum ratio of  <dig> v <dig> <cit> .
however, in , treatment peptide synthesis was expected to be higher than in ctrl treatments since in mammals protein synthesis is reported to increase with increasing protein intake   <cit> .
curiously, on the contrary, in both flux distributions g and h, the total mammary protein synthesis decreased for the  treatment when compared to the  treatment: v <dig> equals  <dig> or  <dig> mmol/h/half udder in  and  <dig> or  <dig> mmol/h/half udder in .
study of quantitative contributions of precursors  for plausible extreme flux distributions
as a further investigation to check the relevance of distributions g and h in the  and  treatments, we studied the quantitative contributions of precursors of output nutrients.
this study was inspired by the usual techniques in the field of nutrition - or any domain concerned with organ studies.
these techniques consist in computing yield rates to elucidate how an input nutrient may contribute to the composition of an output product, for instance to clarify what proportion of glucose, acetate or alanine taken up by the mammary gland can be recovered in the milk components  or oxidized and recovered in co <dig> released in blood.
to formalize this issue, we first selected carbon as the component according to which the contributions of precursors were to be computed.
then, in order to determine how much carbon introduced into the system through a given input flux can be recovered in the rate of production of an output metabolite, we introduced a precise model for the allocation of nutrients in measured outputs  ).
as shown in tables  <dig> and  <dig>  the nutrient allocation corresponding to both flux distributions g and h in the  and  treatments provides evidence that glucose is the unique precursor of lactose synthesis.
this contradicts studies of dairy cows suggesting that glycerol and perhaps amino acids contribute to lactose synthesis  <cit> , and that fifteen percent of lactose carbon could not derive from glucose  <cit> .
both models have empty flux through the reactions oaa → pyr , oaa → g3p  and g3p → g6p .
model  shows strong nadph oxidation whereas model  has zero nadph oxidation.
both models have empty flux through the reactions oaa → pyr , oaa → g3p  and g3p → g6p .
model  shows strong nadph oxidation whereas model  has zero nadph oxidation.
this precise analysis of the origin of carbon in lactose synthesis with aio suggests that extreme flux distributions g and h have to be rejected, although both flux distributions are consistent with flux variability criteria.
we concluded that none of the extreme vertices of the set of plausible flux distributions could be considered biologically relevant with respect to the model and data at hand.
notice that the optimization of any linear combination of metabolic fluxes is either reached by these extreme distributions or reached by the infinite number of flux distributions lying in a face of the simplex.
this suggests that the functioning of the mammary gland cannot be uniquely modeled by the optimization of any linear combination of metabolic fluxes involved in the current model.
discriminate treatments despite mammary gland flexibility: maxima of aio on the complete polyhedron of flux distributions
previous studies suggest that the response of the mammary gland cannot be modeled uniquely by the optimization of a linear objective function of fluxes.
however, there are several nonoptimal flux distributions that satisfy the literature-based information that we have used so farb.
more generally, many flux distributions compatible with the additional constraints, including the condition on carbon precursors for lactose can be shown for both  and  treatments.
topological arguments prove that an infinite set of flux distributions exist.
nonetheless, without an idea about the exact shape of the space of feasible fluxes , we cannot select a plausible point within this space.
in other words, the available knowledge appears insufficient to determine uniquely a flux distribution of nutrients among the different branches of the proposed model.
to understand the functioning of the mammary gland and despite this difficulty, we introduced a method to estimate the variability of nutrient allocation among pathways on a carbon basis ) by computing the range  of aio coefficients.
as these coefficients are nonlinear functions of flux variables, computing these min and max over the complete space of plausible flux distribution requires solving nonlinear optimization problems.
our dedicated algorithm detailed in the method section scaled properly to the real case that we are studyingc.
interestingly, with this approach, we did not favor any internal functioning of the system since we parsed all objective functions , in order to have a complete description of the space of plausible flux distributions.
the min-max tables for allocations of nutrients in the different pathways provided a clearer view of nutrient utilization within the mammary gland.
unlike the functioning of the extremal distributions  and  shown in tables  <dig> and  <dig>  the contribution of carbon from glucose to lactose carbon was quite variable within the full set of plausible distribution.
indeed, as shown in tables  <dig> and  <dig>  glucose was the precursor of 85%  to 100% of the carbon lactose in the  and  treatments.
therefore, there exist flux distributions such that glucose is not the only precursor of lactose synthesis  <cit> , and these distributions are consistent with quantitative estimations of the ratio of carbon glucose recovered in lactose  <cit> .
1
1
1
1
1
1
1
1
a local-search algorithm allowed us to compute the minima and maxima of each aio coefficient for the two treatments ,  .
these tables allow discriminating the response of the mammary gland to the two treatments without requiring selection of a flux distribution for reactions in the metabolic network.
treatment  is characterized by a lower proportion of glucose which is oxidized in co <dig> than in .
amino acid input corresponded to positive balances between amino acid net uptake and amino acid and utilization in milk protein .
amino acid output corresponded to negative balance between, amino acid net uptake and utilization in milk protein .
1
1
1
1
1
1
a local-search algorithm allowed us to compute the minima and maxima of each aio coefficient for the two treatments ,  .
these tables allow discriminating the response of the mammary gland to the two treatments without requiring selection of a flux distribution for reactions in the metabolic network.
treatment  is characterized by a lower proportion of glucose which is oxidized in co <dig> than in .
amino acid input corresponded to positive balances between amino acid net uptake and amino acid and utilization in milk protein .
amino acid output corresponded to negative balance between, amino acid net uptake and utilization in milk protein .
more generally, the intervals of distributions of nutrients in different pathways were used to compare the effects of the  and  treatments.
biologically, in comparison to the  treatment, the  treatment was characterized by a lower proportion of glucose  which is oxidized in co <dig>  and a larger ratio used for lactose synthesis.
this hypothesis can be sustained since the intervals of distribution of carbon from glucose in co <dig> and lactose almost did not overlap in the  and  treatments.
more precisely, in the  treatment, from 3961422= <dig> % to 5651422= <dig> % of glucose carbon was oxidized to produce co <dig>  whereas this ratio ranged from  <dig> % to  <dig> % for the  treatment.
in addition, in the  treatment,  <dig> % to  <dig> % of the glucose carbon was required to produce lactose carbon, whereas  <dig> % to 72% of glucose was required during the  treatment.
altogether, this analysis allowed us to discriminate the effects of the different treatments whatever the internal functioning of the system may be: the  treatment  was characterized by a lower proportion of glucose oxidized in co <dig> than in .
it appears to be a suitable strategy to analyze the metabolism flexibility without selecting a precise flux distribution or making any assumption on the internal metabolic fluxes.
as a final study, we used an interior point exploration method to estimate the aio variability over the boundary of the simplex, rather than the complete space that was explored previously.
as both tables appeared to be equal, we concluded that optimized aio are reached on the boundary of the simplex.
however, not all extreme vertices of the simplex are relevant biologically and they do not optimize aio.
this suggests that the flux distributions that optimize an aio coefficient are placed at the interior of the simplex faces.
in the future, it may be interesting to study the biological significance with multi-objective approaches  <cit>  and check whether these can be considered as a “characteristic” point of these faces.
impact of long-chain fatty acids oxidation over the model predictions
as a last study, we studied the robustness of our conclusion with respect to changes in the modeling of input  of long-chain fatty acids , c <dig> and c <dig>  in the tricarboxylic acid  cycle.
in the study of the  and  treatments, such an input of lcfa was not considered in the system since they are not synthesized by the mammary gland and we hypothesized that they are not oxidized within the mammary gland, based on isotope measurements in studies in fed lactating goats and in nonruminants  <cit> .
this hypothesis was sustained by the measured carbon balance in the  and  datasets  <cit> , which did not require increasing the prediction of co <dig> by introducing any lcfa in the tca cycle .
nonetheless, in other contexts and models, lcfa oxidation was either measured in starvation condition  <cit>  or introduced in order balance carbon consumption and production with respect to co <dig> prediction  <cit> .
this lcfa oxidation hypothesis was for instance retained in models studying the  dataset  <cit> .
to study the impact of the variability of fa oxidation, a ratio of long-chain fa  was introduced in the tca cycle.
for datasets ,  and  which were compatible with a given ratio of lcfa oxidation, extreme flux distributions and aio coefficients variability were studied.
our main conclusions are robust to the introduction of lcfa oxidation .
interestingly, as shown in table  <dig>  the structure of the simplex generated by the  diet is more complicated than the  and  treatments, with  <dig> vertices for all hypotheses.
this is due to the fact that r <dig> is highly constrained by measurements, so that this flux cannot vanish when r <dig> is optimized or when r <dig> is maximized.
all the vertices for the -simplex contradict knowledge-based literature.
r
19
r
15
r
14
r
8
r
64
r
63
r
13
→
to study the impact of the variability of fa oxidation, a ratio of long-chain fa  was introduced in the tca cycle.
for datasets ,  and  which were compatible with a given ratio of lcfa oxidation, extreme flux distributions and aio coefficients variability were studied.
our main conclusions are robust to the introduction of lcfa oxidation .
interestingly, as shown in , the structure of the simplex generated by the  diet is more complicated than the  and  treatments, with  <dig> vertices for all hypotheses.
this is due to the fact that r <dig> is highly constrained by measurements, so that this flux cannot vanish when r <dig> is optimized or when r <dig> is maximized.
all the vertices for the -simplex contradict knowledge-based literature.
from this literature study, it appears that lcfa oxidation may depend on both the environmental and experimental contexts, and no single model can be favored yet.
to study the impact of lcfa oxidation, we successively introduced a ratio of lcfa in the tca cycle , and assumed that 50% of c <dig> output are synthesized within the mammary gland  <cit>  .
then we performed a complete study of the ,  and  datasets for these hypotheses .
for the  treatment, by comparing the predicted and measured co <dig> quantities, we concluded that the hypotheses of 20% and 25% of fa oxidized in the tca cycle were not in agreement with our experimental data: the increase of co <dig> prediction was too large both when compared to measured co <dig>  and when confronted to the measured increase of 9% of co <dig> deriving from long lcfa in lactating goats in an extreme starvation condition  <cit> .
similarly, the hypothesis of 0% of lcfa oxidation was not consistent with the  treatment  <cit> .
for all remaining compatible pairs of model and datasets, we first studied the atp maximization hypothesis.
in all cases, our results, shown in table  <dig>  suggest that atp maximization is not biologically relevant: atp balance was even larger than the atp balance of the models where no lcfa was introduced in the tca cycle .
then we enumerated extreme flux distributions and studied their biological relevance.
as shown in table  <dig>  for each ratio of lcfa in the tca cycle, the structure of the simplex of plausible flux distributions for the  and  datasets was similar to that shown in table  <dig> .
using similar arguments, no extreme distribution could be considered as biologically relevant.
when studying the space of solution associated with the  dataset, a more complex structure appeared, based on  <dig> extreme flux distributions instead of  <dig>  explicit  linear functions that can be optimized to obtain these extreme distributions are detailed in table 9d.
more precisely, for the  treatment, we still obtain a division according to the activation of the four fluxes oaa →pyr , oaa →g3p , g3p →g6p , peptide hydrolysis .
nonetheless, an intricate phenomenon appears.
indeed, measurements imply that oaa →g3p  is very constrained.
it cannot vanish at the same time as oaa →pyr  or when nadph oxidation  is at its maximal value.
therefore, for several optimal conditions related to nadph oxidation , oaa →pyr , g3p →g6p , peptide hydrolysis , we need to decide whether oaa →g3p  is slightly non-zero, or whether oaa →g3p is assumed to vanish although the optimal is not exactly reached for the initial flux considered.
despite this difference of structure between the  dataset and the  and  datasets, our analysis suggested that extreme distributions were not biologically relevant, independently from the ratio of lcfa oxidation or the dataset under study.
finally, comparing the variability of the aio coefficients for the  and  treatments when oxidizing 10% of lcfa in the tca cycle still suggested that the  treatment is characterized by a lower proportion of glucose  which is oxidized in co <dig>  and a larger ratio used for lactose synthesis.
altogether, this study suggests that the main characteristics of the  and the  treatments are robust and could be elucidated despite lacking information on the precise internal behavior of lcfa oxidation.
discussions
analyzing the distribution of nutrients in a metabolic network to study the flexibility of a metabolism at the organ level
an important challenge in applicational fields of metabolism studies at the organ level is to understand how the components of inputs are transformed into some expected outputs, under some assumptions about the functioning of a system.
to that end, great use is made of comparisons between yield rates describing the allocation of input nutrients within the set of outputs.
nonetheless, to allow a precise comparison of nutrients, these studies require insights on the distribution of matter components across the output of each reaction involved in the system.
such information may be provided by several experimental marking techniques  that make use of carbon isotopes  <cit> .
nonetheless, their application to mammals is challenging  <cit> .
to elucidate how input nutrients are allocated among the output nutrients of a metabolic system despite experimental limitations, we have introduced novel methods which refine the flux balance analysis of a metabolic system related to an organ of a large animal.
our method can be seen as an extension of flux variability analysis  <cit> , where the emphasize is put on efficiency or yield rate variability rather than on flux variability.
as an example of application, we have studied the mammary metabolism in ruminants   <cit> .
compared to the conventional models available, our stoichiometric model describes very precisely the variations of energy consumptions as atp, allowing us to investigate some optimization hypotheses related to energy variations.
as a methodological innovation, we introduced a method to estimate the allocation of inputs in outputs , that is, the ratio of transformation of each input nutrient into outputs, provided that we are given a flux distribution balancing the production and consumption of intermediary metabolites.
two aio computations are possible: one may fix a linear combination of flux variable to optimize, leading to an extreme flux distribution  whose allocations  are precisely computed with our method.
alternatively, we can reason with regard to the possibly infinite complete convex space of plausible flux distributions, in the spirit of fva analysis, without solving the set of equations.
in this case, our method estimates the minima and maxima of each allocation  coefficient within the complete space.
the latter point is where our main algorithmic innovation lies.
indeed, aio coefficients are not linear with respect to the flux variables v. therefore, they have no reason to be extremal for flux distributions corresponding to vertices, edges or faces of the simplex, unless the gradient of every component of the aio is a nonzero function.
therefore, we have introduced two algorithms allowing us to compute the extremal values of aio coefficients when the flux space is bounded.
the fastest algorithm requires the exhibition of an analytic expression of the aio matrix.
if this is not possible, a local-search algorithm can be used.
these algorithms can be used to check whether the extremal values are reached for flux distributions lying on the boundary of the simplex of plausible flux distributions.
this approach applied to a stoichiometry model taking atp variations into account permitted a better understanding of ruminant mammary gland metabolism in comparison to previous studies based on similar models and datasets  <cit> .
for instance, we were able to characterize the differences in the effects of two treatments, such as the quantitative range of the proportion of glucose which is oxidized in co <dig>  or used for lactose synthesis.
this information was derived from the carbon composition of each metabolite, expert knowledge which is easily accessible and is therefore compatible with the experimental limitations regarding mammals.
optimization strategies within tissue or organs
our study allowed us to revisit optimization-based hypotheses on the functioning of the mammary metabolism.
we have provided evidence that flux distributions corresponding to an optimal production of energy  cannot describe appropriately the metabolism of the mammary gland, as would be the case for bacterial metabolisms, which tend to optimize biomass-related functions  <cit> .
the main reason is that the atp balance is too large and variable among the responses to different treatments, confirming previous observations  <cit> .
as a consequence, the underlying hypotheses used to drive previous studies of mammary metabolism have to be carefully reformulated  <cit> .
as an alternative, we have hypothesized that atp balance remains nonoptimal and almost constant in response to several treatments, and we have introduced in the model a recent estimation for this quantity  <cit> .
our goal was to check whether the observed responses to the system could be explained by the optimization of a linear combination of fluxes, that is, an extremal vertex of the simplex of plausible flux distributions.
this hypothesis was rejected, first because it led to non-realistic orders of magnitude for some fluxes such as peptide hydrolysis, and second because the precursors of some components such as glucose were not biologically relevant  <cit> .
notably, it was necessary to trace back the quantitative origin of the nutrients to complete this analysis and reject all extreme flux distributions, illustrating the advantage of the aio approach.
therefore, no optimization of a linear combination of flux variables could be found to uniquely describe the metabolism of the mammary gland, a multicellular organ, as an extreme flux distribution of our model.
to overcome this limitation of fba-inspired analyses at the mammalian tissue or organ level, we studied the variability of aio coefficients by introducing the computation of min-max ranges of aio.
this led to a general overview of the effects of treatments without precluding any steady state internal behavior.
benefits of studying the variability of the allocation of input in output  and future model refinements
the study of the variability of aio, that is, intervals of allocation of nutrients, made it possible to distinguish the metabolism of the mammary gland in two different nutritional conditions corresponding to an increase in protein supply .
notably, the intervals of allocation of glucose in co <dig> and lactose were different in the two nutritional conditions, reflecting the mammary gland’s metabolic flexibility.
nonetheless, several improvements can be made to the model.
first, the simplex of plausible distributions could be reduced by introducing knowledge currently available on the kinetic bounds for enzymatic activities, introduced in previous numerical models  <cit> .
this may for instance prevent some non-constrained behaviors shown in tables  <dig> and  <dig>  such as the possible nonzero contribution of acetate and β-hydroxybutyrate  carbon to lactose.
the main issue in this area will be to aggregate the enzymatic knowledge to fit with the format of the reactions in the model, such as the reaction r <dig> which transforms g3p to pyr under the regulation of five different enzymes, together with nad+ and atp availability.
therefore, proposing relevant maximal values requires the use of advanced methods for model reductions.
a second improvement of the model is dependent on the production of additional observations to clarify the set of possible behaviors of the system.
particularly, the contribution of acetate and β-hydroxybutyrate  to lactose shown in tables  <dig> and  <dig> remains questionable in ruminants since it suggests that the carbon of acetyl-coa could leave the tca cycle of the mitochondria through the oaa → g3p pathway.
this could be explained by the fact that the carbon allocation in the model did not take account of the carbon positions.
however, one isotope model and one study of gluconeogenesis in man suggested that labeled carbon within glucose came from -acetate   <cit> .
it remains challenging to obtain such precise information for large ruminants such as dairy cows.
nevertheless, the model could be improved by modeling with additional constraints the few items of knowledge we have on the flux distribution of positioned carbons.
the last improvement of the model consists of including constraints that do not correspond to rule-based metabolic equations.
more precisely, in several numerical models, the effects of external fluxes were introduced to predict the response of the mammary gland in a consistent way  <cit> .
among external fluxes, we noticed the regulatory effect of long-chain fatty acids on the activity of enzymes involved in fatty acid synthesis, although the mechanisms are not clearly understood  <cit> .
although some extensions of the fba formalism to dynamic regulation exist  <cit> , the data availability was insufficient to apply the extensions in the model.
therefore, existing regulations of enzyme activities cannot be included directly in the stoichiometric model.
to overcome this limitation, we plan to investigate the effects of external fluxes by introducing additional constraints in the model, to check whether the nutrient output in the milk can be predicted from the nutrient input.
CONCLUSIONS
we have introduced a method in the framework of flux-balance analysis of a metabolic network.
as a main novelty, our approach allows studying the variability of efficiencies  of a metabolic model provided with input-output measurements.
more precisely, our approach allows a quantitative estimation of the minimum and maximum proportions of the carbon quantity of each input nutrient which is recovered in each output component of the system.
the main innovation is to propose a method which does not require determining the quantitative distribution of nutrients between the branches of the system.
to that end, we have performed a parsing of the space of flux distributions which are compatible with both the model stoichiometry and input-output measurements.
this method was applied to study the response of the mammary gland to several treatments.
it allowed us to distinguish two different metabolic responses of the system, corresponding to two nutritional situations and accurately reflecting metabolic flexibility.
overall, our method appears to be configured to study the variability of the yield rates of a metabolic system at the multicellular or organ level without making any hypothesis on the internal behavior of the system.

 when studying metabolism, it is important to elucidate how fluxes are distributed among the different pathways of the metabolic network with respect to the available quantitative information about the system behavior. several methods can be used to address this issue. the first approach consists of building a mechanistic description of transformations and identifying the regulations involved in the system. continuous dynamical models are often used for this purpose, especially when time-series responses to different treatments are available to infer the dynamics of the network. static approaches such as petri net can also identify qualitative distributions of fluxes in a metabolic network, and their stochastic extensions can even take into account stoichiometric and kinetic information. with a complementary approach, one can study a system at the functional level, based on the study of fluxes at steady states. this is the purpose of the flux balance analysis framework, which has evolved considerably over the past decades. with an fba, the identification of external regulations is not necessary because it is assumed that the global behavior of the system can be modeled by optimizing linear combinations of selected fluxes. roughly, the methods developed in this field aim to explore a convex space of plausible flux distributions and to study the extreme flux distributions obtained when optimizing a linear objective function. it allows checking whether an extreme flux distribution is consistent with experimental data and to predict new experimental observations. nonetheless, the consistency of the solutions obtained by fba depends on the quality of the constraints integrated in the model. to overcome this limitation, several extensions have been proposed in the literature. these extensions can be broken down into two parts. the first incorporates additional biological knowledge, such as reaction thermodynamics or multioptimization. the second is based on the use of fba to globally analyze large-scale metabolic networks. for example, in flux variability analysis, the minimum and the maximum flux for each reaction in the network are computed under some optimal conditions. in this paper, our main purpose is to extend this framework to study the variability of a metabolic network at the level of efficiencies instead of fluxes. indeed, when studying a metabolism at the organ level, a major challenge consists in comparing the efficiency, or yield rates, of two metabolic situations, i.e. the response to various input patterns. a typical example of such studies are those concerning animal nutrition, which aim to predict the quality and quantity of animal production in response to breeding factors. in this field, the energy and protein conversion efficiencies are derived from the study of the flux distribution of input nutrients between the different branches of the metabolic network. more generally, although its definition depends on the field of application, the concept of efficiency is often linked to energy, mass growth and protein conversion. however, the computing of efficiencies is prone to difficulties and errors, since it requires computing the quantities precursors sets which are required to explain the measured composition of outputs, although some of the internal products may be recycled within cycles. with this goal, we defined the allocation of an input towards an output to be the proportion of a matter component in a given input flux that is recovered in the selected output flux. our definition is based on the choice of a material component, such as carbon or nitrogen, which allows a comparison of the contributions of input metabolites to the composition of output products. it can be seen as a yield rate, which is uniform among all the outputs of the system, allowing a precise understanding of the precursors used. as a first methodological contribution, we prove that aio can be uniquely described and computed, even in the case when there are metabolic system cycles, by a matrix whose coefficients are nonlinear functions of the flux variables. introducing nonlinearity in the definition of aio cannot be avoided because of the presence of these cycles. studying the variability of aio within a complete space of plausible flux distributions requires the solving of nonlinear optimization problems which are underdetermined in tangible applications. as a second methodological contribution, we have proposed efficient algorithms to compute lower and upper bounds for aios over the family of flux distributions which are compatible with both the system’s stoichiometry and the experimental datasets, regardless of the choice of a flux distribution for the internal branches of the network. an important aspect is that when the metabolic network is provided with input-output data, the complete space of plausible distributions appears to have a relatively small size, and can therefore be studied with our method. our framework is depicted in figure our main example of application is related to milk production. in this context, several models have been introduced, in relation with the aforementioned classification of models. one class of small-size dynamical mechanistic models predicts the blood flow and input nutrients of the metabolic system or, alternatively, the nutrients produced by the metabolic system. another class of models predicts the distribution of fluxes over the pathways from the input and output nutrients of the system. in the latter family, both dynamic and static approaches exist. indeed, numerical models based on mass-action equations were initially proposed to describe the fluxes to and from individual metabolites, with different levels of description. to that goal, optimizers were used to determine a reasonable set of parameter estimates for the dynamical model of the system. although such a set of parameter estimates is not unique. with a totally different approach, in another study, the main biochemical pathways in the different metabolisms were integrated in a generic stoichiometric model called the metabolic spreadsheet. its structure was based on a restricted number of intermediary metabolites called carbon-chain pivots, which correspond to important cross-over points between metabolic pathways. this allowed the computation of the flux rates for all reactions, constrained by a general rule stating that there is no accumulation of intermediary metabolite of cofactors. a study of manual calculations performed with the metabolic spreadsheet showed that it works, in practice, by maximizing an objective function in a convex solution space. therefore, this model can be considered as an application of the flux balance analysis framework. nonetheless, in this field of study, there has been little discussion on the impact of the choice of a single model among several reasonable models. to investigate this issue in a more automatic way, we compared the aforementioned conventional models to the convex space of plausible flux distributions associated with steady states of a model of mammary gland metabolism, as computed in fba. we checked the consistency of extreme flux distributions of nutrients with experimental data of the mammary gland and milk production, including the contribution of nutrient input-output and isotope balance studies. based on our aio computation framework, our analysis highlighted that the metabolic behavior of the mammary gland cannot be modeled by maximizing atp production or by optimizing a linear combination of flux variables of the model of the mammary gland. in other words, although an infinite number of flux distributions are compatible with the data, none are extreme within the space of feasible models. selecting any of these nonoptimal flux distributions to predict the system behavior appears difficult without additional experimentation. to gain better understanding of the system response regardless of the choice of a flux distribution for the internal branches of the network, we applied our method to estimate the variability of aio coefficients in our model and compared the effects of two different diets on mammary gland metabolism. our results suggest that the bounds of aio are sufficient to distinguish the effects of different nutritional treatments without selecting a flux distribution for the internal reactions of the metabolic network by any method - optimization of a linear combination of fluxes or a residual score. overall, the complete study suggests considerable flexibility in mammary gland metabolism. it provides a view of the functioning of the system although its internal processes still cannot be clarified because of limitations on experimentation on large animals such as ruminants.  we first investigated the set of flux distributions that are compatible with the stoichiometry of our mammary gland model, without taking datasets into account. the model exhibited a large variability, since thousands of extreme pathways could be identified. we then successively computed the set of flux distributions compatible with the model and the real datasets of lactation metabolism in dairy cows. the datasets are given in table they include a control diet, a diet related to an increased protein supply through casein infusion into the duodenum, and a complementary dataset previously used in a mechanistic model of the mammary metabolism. for the three datasets, the set of plausible flux distributions – solutions to a linear system detailed in eq. – was an unbounded convex cone of dimension in all cases, the solution spaces shared the same set of five independent variables. this suggests that there are five independent levels of variability within the system: peptide hydrolysis, nadph oxidation, oaa →pyr, oaa →g3p and g3p →g6p. therefore, by including this dataset the model became a fairly small and constrained network. nevertheless, it was not uniquely determined since there were still several degrees of freedom. fatty acid output amino acids balance fa primer from acetate data are renormalized in mmol/h/half udder. control diet higher protein supply by casein infusion in the duodenum generic dataset. this table is used to parameterize input and output vectors v i and v o together with additional biological linear constraints on some reaction fluxes. 1input i.e. taken up by the stoichiometric system considered. 2 β -hydroxybutyrate. 3total triglycerides secreted in milk considering that milk fat was composed of 100% triglycerides and that all the triglycerides were secreted in milk fat. 4output i.e. leaving the system. 5all fatty acids synthesized within the mammary gland i.e. all c to c and 50% of c. 6the balance between amino acid net uptake and amino acid net output in milk protein is calculated with established rules. if the balance is positive it corresponds to an amino acid input that will be catabolized. if this balance is negative, the values corresponded to an amino acid output. 7serine output corresponded to serine synthesized minus serine utilized in other pathways i.e. serine required in addition to ser uptake to synthesize milk protein. 8peptide output: number of peptide links required to synthesize the proteins exported out of the system. 9hanigan,. 10fatty acid primers were synthesized for 50% from acetate and for 50% from bhba except c fa primer which was supposed to be synthetized only from bhba. 11set at zero because their inputs or outputs are set at zero. investigating the relevance of the optimization strategies for mammary metabolism the balance between the atp generated by the system and the atp used by the system was computed for the three datasets, and. the results are detailed in table in this table, the fba approach based on optimization of the atp balance is compared to the natural functioning of the system. criteria of selection of a solution three natural assumptions are considered to model the mammary gland behavior: removal of all cycles, optimization of atp production and study of the equilibria of a dedicated ode-based model. all models exhibit considerable variability in their atp balance, which contradicts the assumption about the behavior of this organ. moreover, quantitatively, the computed atp balances are much higher than recent measurements. this suggests that atp maximization cannot be considered as a natural objective function to model cow mammary behavior. first, we considered the manual computation of fluxes with a tool named “metabolic spreadsheet”. we applied the rules of no accumulation on any intermediary metabolites and cofactors with a utilization rejecting all the cycles. as expected, since some cycles are atp-consuming, both flux distributions were equivalent in both approaches. the atp balances were respectively and mmol/h/half udder in the, and datasets. according to this model, the atp generated is estimated to be mmol/h/half udder while mmol/h/half udder were estimated to be used for milk component synthesis. therefore, the atp-balance obtained for this model was slightly smaller than the optimal one obtained for our mammary gland model. independently of the approach considered and the model or dataset at hand, the remaining atp, after use for milk synthesis, appeared to be rather high and nonxconstant. indeed, as shown in, an atp balance of mmol/h/half udder can be expected to be used for other functions not accounted for in the models, such as maintaining membrane potential and synthesizing nucleic acids. the numerical values obtained here were far from this estimated balance, suggesting that each of these three models have a nonrelevant atp balance. in addition, both in the ode model and the atp-optimization approach, peptide hydrolysis was obtained at zero, implying an absence of any protein turnover. this contradicts all the observations about this pathway: considerable use of this pathway has been evidenced in several publications, although the peptide hydrolysis rates differed significantly depending on the technique used for the measurements: peptide hydrolysis spans from, to of peptide synthesis in, and, respectively. overall, we concluded with this analysis that the energy-based optimization function may not allow an appropriate simulation of mammary gland metabolism. this was expected, considering that the system is studied at the complete organ level, involving competing processes which can rarely be modeled with a single linear objective function. exploring all extreme flux distributions in a refined simplex in order to study the variability within the space of plausible flux distributions and to identify alternative relevant optimization strategies, an additional constraint was placed on the atp balance of the system. we considered that an atp-balance of mmol/h/half udder was a relevant measure for this study, although we checked that all the results discussed below were still valid when introducing a 10% tolerance on this estimation of atp-balance. providing a bound for atp yielded a new constraint on several variables of the system. for the two datasets and, the flux g3p →g6p was no longer an independent variable in the system. the simplex, which was previously unbounded, appeared to be bounded with four independent variables: oaa →pyr, oaa →g3p, nadph oxidation, peptide hydrolysis. applying flux variability tools, it appeared that the minimum of each of the fluxes v v v v was equal to zero for all treatments. the maxima of the fluxes were for the treatment and for the treatment this suggests that the treatment generates a more flexible space of plausible flux distributions than the treatment. we computed all extreme vertices of the simplex of plausible flux distributions for the two treatments,. the simplex structure of the polyhedron implies that the optimum of any linear combination of metabolic fluxes involved in the model is either uniquely attained for one of these extreme points or attained by all flux distributions positioned on a face of the simplex. to gain insight on the pathways involved in the variability of our model, we also computed the linear combination of fluxes optimized for one of the extreme flux distributions. eight extreme vertices were found for the and datasets. they were obtained when optimizing the same set of linear functionsa. in table the optimal flux distributions for, are classified according to the activation or inactivation of several pathways in the model. it is worth noting that, according to this table, the two sets of plausible flux distributions associated with and have the same topological structure. indeed, the eight extremal behaviors for and have clear combinatorics: first, nadph oxidation is set either at zero or at its maximal value. then, a single flux within oaa →pyr, oaa →g3p, g3p →g6p, peptide hydrolysis is strongly activated whereas the three other remaining fluxes are blocked. r 19 r 15 r 14 r 8 r 64 r 63 r 13 → the qualitative properties of all vertices are shared in and treatments. both correspond to a simplex with height vertices. so are six of each of the and simplex vertices. h and g vertices, in the and treatments, are plausible with respect to the literature. as discussed in a previous paragraph, flux distributions with no peptide hydrolysis cannot be considered as relevant. this suggests the extreme flux distributions a-f in and are not biologically relevant. on the contrary, two extreme flux distributions are consistent with the stoichiometry of the system for the and treatments. these two flux distributions consist in optimizing peptide synthesis and hydrolysis after nadph oxidation. they corresponded to a ratio between peptide hydrolysis and synthesis v64v of in the treatment and of in the treatment. these ratios are equal or lower to the maximum ratio of v. however, in, treatment peptide synthesis was expected to be higher than in ctrl treatments since in mammals protein synthesis is reported to increase with increasing protein intake. curiously, on the contrary, in both flux distributions g and h, the total mammary protein synthesis decreased for the treatment when compared to the treatment: v equals or mmol/h/half udder in and or mmol/h/half udder in. study of quantitative contributions of precursors for plausible extreme flux distributions as a further investigation to check the relevance of distributions g and h in the and treatments, we studied the quantitative contributions of precursors of output nutrients. this study was inspired by the usual techniques in the field of nutrition - or any domain concerned with organ studies. these techniques consist in computing yield rates to elucidate how an input nutrient may contribute to the composition of an output product, for instance to clarify what proportion of glucose, acetate or alanine taken up by the mammary gland can be recovered in the milk components or oxidized and recovered in co released in blood. to formalize this issue, we first selected carbon as the component according to which the contributions of precursors were to be computed. then, in order to determine how much carbon introduced into the system through a given input flux can be recovered in the rate of production of an output metabolite, we introduced a precise model for the allocation of nutrients in measured outputs ). as shown in tables and the nutrient allocation corresponding to both flux distributions g and h in the and treatments provides evidence that glucose is the unique precursor of lactose synthesis. this contradicts studies of dairy cows suggesting that glycerol and perhaps amino acids contribute to lactose synthesis, and that fifteen percent of lactose carbon could not derive from glucose. both models have empty flux through the reactions oaa → pyr, oaa → g3p and g3p → g6p. model shows strong nadph oxidation whereas model has zero nadph oxidation. both models have empty flux through the reactions oaa → pyr, oaa → g3p and g3p → g6p. model shows strong nadph oxidation whereas model has zero nadph oxidation. this precise analysis of the origin of carbon in lactose synthesis with aio suggests that extreme flux distributions g and h have to be rejected, although both flux distributions are consistent with flux variability criteria. we concluded that none of the extreme vertices of the set of plausible flux distributions could be considered biologically relevant with respect to the model and data at hand. notice that the optimization of any linear combination of metabolic fluxes is either reached by these extreme distributions or reached by the infinite number of flux distributions lying in a face of the simplex. this suggests that the functioning of the mammary gland cannot be uniquely modeled by the optimization of any linear combination of metabolic fluxes involved in the current model. discriminate treatments despite mammary gland flexibility: maxima of aio on the complete polyhedron of flux distributions previous studies suggest that the response of the mammary gland cannot be modeled uniquely by the optimization of a linear objective function of fluxes. however, there are several nonoptimal flux distributions that satisfy the literature-based information that we have used so farb. more generally, many flux distributions compatible with the additional constraints, including the condition on carbon precursors for lactose can be shown for both and treatments. topological arguments prove that an infinite set of flux distributions exist. nonetheless, without an idea about the exact shape of the space of feasible fluxes, we cannot select a plausible point within this space. in other words, the available knowledge appears insufficient to determine uniquely a flux distribution of nutrients among the different branches of the proposed model. to understand the functioning of the mammary gland and despite this difficulty, we introduced a method to estimate the variability of nutrient allocation among pathways on a carbon basis ) by computing the range of aio coefficients. as these coefficients are nonlinear functions of flux variables, computing these min and max over the complete space of plausible flux distribution requires solving nonlinear optimization problems. our dedicated algorithm detailed in the method section scaled properly to the real case that we are studyingc. interestingly, with this approach, we did not favor any internal functioning of the system since we parsed all objective functions, in order to have a complete description of the space of plausible flux distributions. the min-max tables for allocations of nutrients in the different pathways provided a clearer view of nutrient utilization within the mammary gland. unlike the functioning of the extremal distributions and shown in tables and the contribution of carbon from glucose to lactose carbon was quite variable within the full set of plausible distribution. indeed, as shown in tables and glucose was the precursor of 85% to 100% of the carbon lactose in the and treatments. therefore, there exist flux distributions such that glucose is not the only precursor of lactose synthesis, and these distributions are consistent with quantitative estimations of the ratio of carbon glucose recovered in lactose. 1 1 1 1 1 1 1 1 a local-search algorithm allowed us to compute the minima and maxima of each aio coefficient for the two treatments,. these tables allow discriminating the response of the mammary gland to the two treatments without requiring selection of a flux distribution for reactions in the metabolic network. treatment is characterized by a lower proportion of glucose which is oxidized in co than in. amino acid input corresponded to positive balances between amino acid net uptake and amino acid and utilization in milk protein. amino acid output corresponded to negative balance between, amino acid net uptake and utilization in milk protein. 1 1 1 1 1 1 a local-search algorithm allowed us to compute the minima and maxima of each aio coefficient for the two treatments,. these tables allow discriminating the response of the mammary gland to the two treatments without requiring selection of a flux distribution for reactions in the metabolic network. treatment is characterized by a lower proportion of glucose which is oxidized in co than in. amino acid input corresponded to positive balances between amino acid net uptake and amino acid and utilization in milk protein. amino acid output corresponded to negative balance between, amino acid net uptake and utilization in milk protein. more generally, the intervals of distributions of nutrients in different pathways were used to compare the effects of the and treatments. biologically, in comparison to the treatment, the treatment was characterized by a lower proportion of glucose which is oxidized in co and a larger ratio used for lactose synthesis. this hypothesis can be sustained since the intervals of distribution of carbon from glucose in co and lactose almost did not overlap in the and treatments. more precisely, in the treatment, from 3961422= % to 5651422= % of glucose carbon was oxidized to produce co whereas this ratio ranged from % to % for the treatment. in addition, in the treatment, % to % of the glucose carbon was required to produce lactose carbon, whereas % to 72% of glucose was required during the treatment. altogether, this analysis allowed us to discriminate the effects of the different treatments whatever the internal functioning of the system may be: the treatment was characterized by a lower proportion of glucose oxidized in co than in. it appears to be a suitable strategy to analyze the metabolism flexibility without selecting a precise flux distribution or making any assumption on the internal metabolic fluxes. as a final study, we used an interior point exploration method to estimate the aio variability over the boundary of the simplex, rather than the complete space that was explored previously. as both tables appeared to be equal, we concluded that optimized aio are reached on the boundary of the simplex. however, not all extreme vertices of the simplex are relevant biologically and they do not optimize aio. this suggests that the flux distributions that optimize an aio coefficient are placed at the interior of the simplex faces. in the future, it may be interesting to study the biological significance with multi-objective approaches and check whether these can be considered as a “characteristic” point of these faces. impact of long-chain fatty acids oxidation over the model predictions as a last study, we studied the robustness of our conclusion with respect to changes in the modeling of input of long-chain fatty acids, c and c in the tricarboxylic acid cycle. in the study of the and treatments, such an input of lcfa was not considered in the system since they are not synthesized by the mammary gland and we hypothesized that they are not oxidized within the mammary gland, based on isotope measurements in studies in fed lactating goats and in nonruminants. this hypothesis was sustained by the measured carbon balance in the and datasets, which did not require increasing the prediction of co by introducing any lcfa in the tca cycle. nonetheless, in other contexts and models, lcfa oxidation was either measured in starvation condition or introduced in order balance carbon consumption and production with respect to co prediction. this lcfa oxidation hypothesis was for instance retained in models studying the dataset. to study the impact of the variability of fa oxidation, a ratio of long-chain fa was introduced in the tca cycle. for datasets, and which were compatible with a given ratio of lcfa oxidation, extreme flux distributions and aio coefficients variability were studied. our main conclusions are robust to the introduction of lcfa oxidation. interestingly, as shown in table the structure of the simplex generated by the diet is more complicated than the and treatments, with vertices for all hypotheses. this is due to the fact that r is highly constrained by measurements, so that this flux cannot vanish when r is optimized or when r is maximized. all the vertices for the -simplex contradict knowledge-based literature. r 19 r 15 r 14 r 8 r 64 r 63 r 13 → to study the impact of the variability of fa oxidation, a ratio of long-chain fa was introduced in the tca cycle. for datasets, and which were compatible with a given ratio of lcfa oxidation, extreme flux distributions and aio coefficients variability were studied. our main conclusions are robust to the introduction of lcfa oxidation. interestingly, as shown in, the structure of the simplex generated by the diet is more complicated than the and treatments, with vertices for all hypotheses. this is due to the fact that r is highly constrained by measurements, so that this flux cannot vanish when r is optimized or when r is maximized. all the vertices for the -simplex contradict knowledge-based literature. from this literature study, it appears that lcfa oxidation may depend on both the environmental and experimental contexts, and no single model can be favored yet. to study the impact of lcfa oxidation, we successively introduced a ratio of lcfa in the tca cycle, and assumed that 50% of c output are synthesized within the mammary gland. then we performed a complete study of the, and datasets for these hypotheses. for the treatment, by comparing the predicted and measured co quantities, we concluded that the hypotheses of 20% and 25% of fa oxidized in the tca cycle were not in agreement with our experimental data: the increase of co prediction was too large both when compared to measured co and when confronted to the measured increase of 9% of co deriving from long lcfa in lactating goats in an extreme starvation condition. similarly, the hypothesis of 0% of lcfa oxidation was not consistent with the treatment. for all remaining compatible pairs of model and datasets, we first studied the atp maximization hypothesis. in all cases, our results, shown in table suggest that atp maximization is not biologically relevant: atp balance was even larger than the atp balance of the models where no lcfa was introduced in the tca cycle. then we enumerated extreme flux distributions and studied their biological relevance. as shown in table for each ratio of lcfa in the tca cycle, the structure of the simplex of plausible flux distributions for the and datasets was similar to that shown in table. using similar arguments, no extreme distribution could be considered as biologically relevant. when studying the space of solution associated with the dataset, a more complex structure appeared, based on extreme flux distributions instead of explicit linear functions that can be optimized to obtain these extreme distributions are detailed in table 9d. more precisely, for the treatment, we still obtain a division according to the activation of the four fluxes oaa →pyr, oaa →g3p, g3p →g6p, peptide hydrolysis. nonetheless, an intricate phenomenon appears. indeed, measurements imply that oaa →g3p is very constrained. it cannot vanish at the same time as oaa →pyr or when nadph oxidation is at its maximal value. therefore, for several optimal conditions related to nadph oxidation, oaa →pyr, g3p →g6p, peptide hydrolysis, we need to decide whether oaa →g3p is slightly non-zero, or whether oaa →g3p is assumed to vanish although the optimal is not exactly reached for the initial flux considered. despite this difference of structure between the dataset and the and datasets, our analysis suggested that extreme distributions were not biologically relevant, independently from the ratio of lcfa oxidation or the dataset under study. finally, comparing the variability of the aio coefficients for the and treatments when oxidizing 10% of lcfa in the tca cycle still suggested that the treatment is characterized by a lower proportion of glucose which is oxidized in co and a larger ratio used for lactose synthesis. altogether, this study suggests that the main characteristics of the and the treatments are robust and could be elucidated despite lacking information on the precise internal behavior of lcfa oxidation. discussions analyzing the distribution of nutrients in a metabolic network to study the flexibility of a metabolism at the organ level an important challenge in applicational fields of metabolism studies at the organ level is to understand how the components of inputs are transformed into some expected outputs, under some assumptions about the functioning of a system. to that end, great use is made of comparisons between yield rates describing the allocation of input nutrients within the set of outputs. nonetheless, to allow a precise comparison of nutrients, these studies require insights on the distribution of matter components across the output of each reaction involved in the system. such information may be provided by several experimental marking techniques that make use of carbon isotopes. nonetheless, their application to mammals is challenging. to elucidate how input nutrients are allocated among the output nutrients of a metabolic system despite experimental limitations, we have introduced novel methods which refine the flux balance analysis of a metabolic system related to an organ of a large animal. our method can be seen as an extension of flux variability analysis, where the emphasize is put on efficiency or yield rate variability rather than on flux variability. as an example of application, we have studied the mammary metabolism in ruminants. compared to the conventional models available, our stoichiometric model describes very precisely the variations of energy consumptions as atp, allowing us to investigate some optimization hypotheses related to energy variations. as a methodological innovation, we introduced a method to estimate the allocation of inputs in outputs, that is, the ratio of transformation of each input nutrient into outputs, provided that we are given a flux distribution balancing the production and consumption of intermediary metabolites. two aio computations are possible: one may fix a linear combination of flux variable to optimize, leading to an extreme flux distribution whose allocations are precisely computed with our method. alternatively, we can reason with regard to the possibly infinite complete convex space of plausible flux distributions, in the spirit of fva analysis, without solving the set of equations. in this case, our method estimates the minima and maxima of each allocation coefficient within the complete space. the latter point is where our main algorithmic innovation lies. indeed, aio coefficients are not linear with respect to the flux variables v. therefore, they have no reason to be extremal for flux distributions corresponding to vertices, edges or faces of the simplex, unless the gradient of every component of the aio is a nonzero function. therefore, we have introduced two algorithms allowing us to compute the extremal values of aio coefficients when the flux space is bounded. the fastest algorithm requires the exhibition of an analytic expression of the aio matrix. if this is not possible, a local-search algorithm can be used. these algorithms can be used to check whether the extremal values are reached for flux distributions lying on the boundary of the simplex of plausible flux distributions. this approach applied to a stoichiometry model taking atp variations into account permitted a better understanding of ruminant mammary gland metabolism in comparison to previous studies based on similar models and datasets. for instance, we were able to characterize the differences in the effects of two treatments, such as the quantitative range of the proportion of glucose which is oxidized in co or used for lactose synthesis. this information was derived from the carbon composition of each metabolite, expert knowledge which is easily accessible and is therefore compatible with the experimental limitations regarding mammals. optimization strategies within tissue or organs our study allowed us to revisit optimization-based hypotheses on the functioning of the mammary metabolism. we have provided evidence that flux distributions corresponding to an optimal production of energy cannot describe appropriately the metabolism of the mammary gland, as would be the case for bacterial metabolisms, which tend to optimize biomass-related functions. the main reason is that the atp balance is too large and variable among the responses to different treatments, confirming previous observations. as a consequence, the underlying hypotheses used to drive previous studies of mammary metabolism have to be carefully reformulated. as an alternative, we have hypothesized that atp balance remains nonoptimal and almost constant in response to several treatments, and we have introduced in the model a recent estimation for this quantity. our goal was to check whether the observed responses to the system could be explained by the optimization of a linear combination of fluxes, that is, an extremal vertex of the simplex of plausible flux distributions. this hypothesis was rejected, first because it led to non-realistic orders of magnitude for some fluxes such as peptide hydrolysis, and second because the precursors of some components such as glucose were not biologically relevant. notably, it was necessary to trace back the quantitative origin of the nutrients to complete this analysis and reject all extreme flux distributions, illustrating the advantage of the aio approach. therefore, no optimization of a linear combination of flux variables could be found to uniquely describe the metabolism of the mammary gland, a multicellular organ, as an extreme flux distribution of our model. to overcome this limitation of fba-inspired analyses at the mammalian tissue or organ level, we studied the variability of aio coefficients by introducing the computation of min-max ranges of aio. this led to a general overview of the effects of treatments without precluding any steady state internal behavior. benefits of studying the variability of the allocation of input in output and future model refinements the study of the variability of aio, that is, intervals of allocation of nutrients, made it possible to distinguish the metabolism of the mammary gland in two different nutritional conditions corresponding to an increase in protein supply. notably, the intervals of allocation of glucose in co and lactose were different in the two nutritional conditions, reflecting the mammary gland’s metabolic flexibility. nonetheless, several improvements can be made to the model. first, the simplex of plausible distributions could be reduced by introducing knowledge currently available on the kinetic bounds for enzymatic activities, introduced in previous numerical models. this may for instance prevent some non-constrained behaviors shown in tables and such as the possible nonzero contribution of acetate and β-hydroxybutyrate carbon to lactose. the main issue in this area will be to aggregate the enzymatic knowledge to fit with the format of the reactions in the model, such as the reaction r which transforms g3p to pyr under the regulation of five different enzymes, together with nad+ and atp availability. therefore, proposing relevant maximal values requires the use of advanced methods for model reductions. a second improvement of the model is dependent on the production of additional observations to clarify the set of possible behaviors of the system. particularly, the contribution of acetate and β-hydroxybutyrate to lactose shown in tables and remains questionable in ruminants since it suggests that the carbon of acetyl-coa could leave the tca cycle of the mitochondria through the oaa → g3p pathway. this could be explained by the fact that the carbon allocation in the model did not take account of the carbon positions. however, one isotope model and one study of gluconeogenesis in man suggested that labeled carbon within glucose came from -acetate. it remains challenging to obtain such precise information for large ruminants such as dairy cows. nevertheless, the model could be improved by modeling with additional constraints the few items of knowledge we have on the flux distribution of positioned carbons. the last improvement of the model consists of including constraints that do not correspond to rule-based metabolic equations. more precisely, in several numerical models, the effects of external fluxes were introduced to predict the response of the mammary gland in a consistent way. among external fluxes, we noticed the regulatory effect of long-chain fatty acids on the activity of enzymes involved in fatty acid synthesis, although the mechanisms are not clearly understood. although some extensions of the fba formalism to dynamic regulation exist, the data availability was insufficient to apply the extensions in the model. therefore, existing regulations of enzyme activities cannot be included directly in the stoichiometric model. to overcome this limitation, we plan to investigate the effects of external fluxes by introducing additional constraints in the model, to check whether the nutrient output in the milk can be predicted from the nutrient input.  we have introduced a method in the framework of flux-balance analysis of a metabolic network. as a main novelty, our approach allows studying the variability of efficiencies of a metabolic model provided with input-output measurements. more precisely, our approach allows a quantitative estimation of the minimum and maximum proportions of the carbon quantity of each input nutrient which is recovered in each output component of the system. the main innovation is to propose a method which does not require determining the quantitative distribution of nutrients between the branches of the system. to that end, we have performed a parsing of the space of flux distributions which are compatible with both the model stoichiometry and input-output measurements. this method was applied to study the response of the mammary gland to several treatments. it allowed us to distinguish two different metabolic responses of the system, corresponding to two nutritional situations and accurately reflecting metabolic flexibility. overall, our method appears to be configured to study the variability of the yield rates of a metabolic system at the multicellular or organ level without making any hypothesis on the internal behavior of the system. 
alternatively, an infinite family of compatible set of flux rates may have to be studied when the constraints raised by observations are insufficient to identify a single flux distribution. for example, in nutrition, biochemical models have been developed to study the metabolism of the mammary gland in relation to the synthesis of milk components. the precursors of output nutrients are traced back with analyses similar to the computation of yield rates. applied to a model of the metabolism of the mammary gland, our method made it possible to distinguish the effects of different nutritional treatments, although it cannot be proved that the mammary gland optimizes a specific linear combination of flux variables, including those based on energy. by solving nonlinear optimization problems, we introduced a method to study the variability of aio coefficients when parsing the space of flux distributions that are compatible with both model stoichiometry and experimental data. flux balance analysisflux distributions explorationyield variabilitynutritional model in these models, a compatible flux distribution is first selected. it allows a quantitative comparison of the respective contributions of precursors to the production of a set of nutrients by a metabolic network, regardless of the choice of the flux distribution within the different branches of the network. altogether, our study indicated that the mammary gland possesses considerable metabolic flexibility. we formally modeled the quantitative allocation of input nutrients among the branches of the metabolic network. our method enables to study the variability of a metabolic network with respect to efficiency. however, the computation of the quantitative contributions of precursors may lack precision, mainly because some precursors are involved in the composition of several nutrients and because some metabolites are cycled in loops. when studying metabolism at the organ level, a major challenge is to understand the matter exchanges between the input and output components of the system.
 when studying metabolism at the organ level, a major challenge is to understand the matter exchanges between the input and output components of the system. for example, in nutrition, biochemical models have been developed to study the metabolism of the mammary gland in relation to the synthesis of milk components. these models were designed to account for the quantitative constraints observed on inputs and outputs of the system. in these models, a compatible flux distribution is first selected. alternatively, an infinite family of compatible set of flux rates may have to be studied when the constraints raised by observations are insufficient to identify a single flux distribution. the precursors of output nutrients are traced back with analyses similar to the computation of yield rates. however, the computation of the quantitative contributions of precursors may lack precision, mainly because some precursors are involved in the composition of several nutrients and because some metabolites are cycled in loops.  we formally modeled the quantitative allocation of input nutrients among the branches of the metabolic network. it corresponds to yield information which, if standardized across all the outputs of the system, allows a precise quantitative understanding of their precursors. by solving nonlinear optimization problems, we introduced a method to study the variability of aio coefficients when parsing the space of flux distributions that are compatible with both model stoichiometry and experimental data. applied to a model of the metabolism of the mammary gland, our method made it possible to distinguish the effects of different nutritional treatments, although it cannot be proved that the mammary gland optimizes a specific linear combination of flux variables, including those based on energy. altogether, our study indicated that the mammary gland possesses considerable metabolic flexibility.  our method enables to study the variability of a metabolic network with respect to efficiency. it allows a quantitative comparison of the respective contributions of precursors to the production of a set of nutrients by a metabolic network, regardless of the choice of the flux distribution within the different branches of the network. flux balance analysisflux distributions explorationyield variabilitynutritional model 