A Strategy for Development of Realistic Mathematical Models of Whole-Body Metabolism

When realistic mathematical models of whole body metabolism eventually become available, they are likely to add entirely new dimensions to the understanding of the integrated physiological function of the organism, in particular the mechanisms governing the regulation of transitions between different physiological states, like fed-fasted, exercise-rest and normal-diseased. So far the strategy for whole body modelling has primarily been a bottom-up approach where the central problem is an apparently insurmountable barrier of complexity involved in defining and optimising the huge number of parameters. Here we follow a top-down strategy and present a complete mathematical framework for realistic whole body model development. The approach proposed is modular and hierarchical and whole body metabolism is taken as the top level. Next are the organs, where the sum of the contributions from the individual organs must equal the top level metabolism. This hierarchy can be extended to lower levels of organisation, i.e. clusters of cells, individual cells, organelle and individual pathways. Exploiting this hierarchy, metabolism at each level forms an absolute constraint on the contributions from lower level. Importantly, these constraints can in many ways be defined experimentally through mass balance and flux data. Furthermore, the constrained approach allows the lower level models to be developed independently and subsequently adapted to the whole body model. The paper describes the process of whole body modelling in practical terms, centred on a mathematical framework, devised to allow whole-body models of any complexity to be developed. Furthermore, an example of sub-model incorporation in the whole-body framework is illustrated by adapting an existing erythrocyte model to the whole body constraints. Finally, we illustrate the operation of the system by including two sets of whole-body data from humans, reflecting two different physiological states.


Introduction
The development of mathematical models that can describe the overall metabolism of a human will not only be a milestone for biosimulation but is also likely to add completely new dimensions to the understanding of the dynamics of the integrated metabolic and physiological function in health and disease.
The human and other mammalian organisms are naturally compartmentalised, composed of a number of organs with specialised physiological and metabolic functions.Similarly, the organs are composed of several types of cells, often with sub-populations which are specialised to perform specific tasks, like the perivenous hepatocytes, which exclusively perform the glutamine synthesis of the liver [1].Thus, from a modelling point of view, it is important to realize that the body and the organs work in a modular, but concerted, manner.A well-known example of concerted action at the organ level is the interplay be-tween the liver and muscle, termed the Cori cycle, where the muscle delivers lactate and alanine to the liver, which converts lactate to glucose and alanine to urea plus glucose.Subsequently, the glucose may return to the muscle as energy substrate [2].Another more recently observed organ-organ interaction is the functional coupling between brain and muscle, where lactate derived from the muscles during intense exercise is used extensively by the brain as a supplement to the normal substrate, glucose [3].Analogously, during fasting the liver produces significant amounts of ketone-bodies which are excellent substrates for the brain and heart [4,5].
The key outcome of whole-body metabolism modelling will be the piecing together of fragmented biochemical and physiological information to an integrated understanding of the overall metabolic function, and in particular the possibility to interrogate the model about dynamic aspects of metabolism, which for practical reasons are very difficult or impossible to approach experimentally.Thus, whole-body models will allow testing of hypotheses about the complicated regulatory mechanisms that govern metabolic flows, and allow the study of transitions occurring when the body switches from one mode of operation (a physiological state) to another, e.g. during the transition from rest to exercise, from the fed to the starved state, and indeed from a normal to a diseased state.In particular, this can be expected to provide crucial information about dysmetabolic conditions such as type 2 diabetes [6,7] and the intimate involvement of mitochondrial metabolism of this disease [8].
Mathematical models of metabolism have so far mainly been targeted at relatively small sub-networks in the individual cells, like for example erythrocyte metabolism [9,10] or the cell cycle machinery [11].Even at this intracellular level there is extensive modularity in the form of subcompartments like the mitochondria, the nucleus, the Golgi apparatus, the endoplasmic reticulum, etc.This modularity has been observed and utilized extensively by researchers working to establish models of the individual cells, or analyse their behaviour.Examples of such modular approaches are modular kinetic analysis [12][13][14] and modular control analysis [15][16][17].It is very likely that the application of these methods could provide a good basis for guiding modelling and experimentation, also in the context of whole-body modelling.There are some examples of modelling at the wholebody level like pharmacokinetic studies describing the uptake, conversion, and excretion of drugs [18] and pharmacokinetic modelling where drug metabolism at the whole-body level was elucidated [19,20].Several models of type 2 diabetes, of varying complexity, have also been developed [21,22], and multi-organ models have been used in the study of amino acid metabolism [23,24].Also, minimal whole-body models have proven useful for interpreting the result of glucose tolerance tests, but are of a very abstract character, often containing no more than a few variables [25].A recent approach exploits the fact that glucose and insulin flux estimates can be obtained from experiments using minimal modelling in conjunction with isotope labelled glucose [26].These flux estimates can be used to perform parameter estimation on submodules of the whole system using a forcing approach [27].Another purely top-down approach to modelling of whole-body metabolism has been attempted recently [28].Models of the heart have been constructed bridging several levels of complexity [29].While these models have been successful in describing the pumping function of the heart, their integration in a whole-body context is lacking.
In principle, a whole-body model could be constructed simply by assembling models of all the individual processes involved.This is often the modelling strategy pursued in engineering, like for example robotics [30].However, the reality of biochemical systems is that most quan-titative aspects of the processes and regulatory systems are unknown at the molecular level, as is also the exact conditions in vivo under which they are operating.This makes modelling of whole-body function according to a strict bottom-up approach impractical.Even if enough information was available about relevant individual pathways, the assembling of the necessary, very large amounts of data is not likely to produce a whole-body model that mimics the actual observed overall human physiology.It would seem that this was the lesson learned from the "silicon cell" approach, where models constructed according to this approach could not be made to work appropriately [31][32][33][34]-as examined in [33] one is often required to change in vitro enzyme activities many fold in order to mimick in vivo behaviour of yeast glycolytic fluxes.One likely reason for these difficulties is the inevitable accumulation of errors, which will arise during the measurement of the individual parts of the system.Another potential source of problems is the transferability of in vitro measurements to in vivo conditions [35][36][37].One might plan on improving the performance by extensive parameter optimisation.However, as illustrated in Figure 1, the huge number of parts and parameters in such a model will lead to a combinatorial explosion that will render the task of parameter fitting very difficult, if not impossible.On the other hand, the silicon-cell approach has the great advantage of building directly on a vast amount of biochemical knowledge.This sets out a straight-forward, rational, method for formulating biochemically realistic model structures.It also allows the modeller to apply standard biochemical methods and data that would have to be discarded with a strict top-down modelling strategy.Consequently, an optimal modelling strategy would combine the best sides of the top-down and silicon cell approaches, and lend from the developments of modular approaches already utilized in the modelling community, in for example engineering where hierarchical modelling has been applied successfully for years [30,38,39]; the purpose of the present paper is to propose such a strategy.
One of the top-down aspects of our strategy is that each dynamic organ model is constrained by net metabolic changes across the organ as defined by arteriovenous differences observed experimentally.Another is that the sum of the model outputs of all organs must equal the total turn over of the whole body.Similarly, the individual organ model can form absolute constraints for the modelling of the cells making up the organ.This principle of constraint can be applied recursively to individual organelles within cells as illustrated in Figure 1.This procedure complements the silicon-cell approach by circumventing the problems of error accumulation encountered here, and ensures an over all realistic model behaviour in the predefined physiological state.The bottom-up aspect of the strategy is that the complexity of the individual organ models can be chosen at will, so that any available silicon-cell model can be included in the metabolic whole body model.In essence, the top-down constraints act like "the cement" that combines individual silicon-cell building blocks.
In practical terms, the modelling strategy developed in this communication, exploits the natural compartmentation into organs to facilitate parameter estimation: The integrated function of each individual organ can in principle be studied in vivo by applying the Fick principle [40] by multiplying arterio-venous concentration differences with blood flow across the organ.This experimental approach has been applied extensively in the study of organ physiology in the past and more recently also to study of brain energy metabolism as well as muscle-brain interactions providing data on the global metabolism of the human brain in vivo [41][42][43] at a time resolution of less than five seconds [44].Combined with stable isotope infusion to the subject during the experiment yet another level of information on the dynamics of metabolism may be obtained [45].
In this paper a formalism that allows this kind of modular, constrained metabolic whole-body modelling to be performed, is presented.Its use is illustrated on metabolic data describing two well-defined physiological states, fed and fasted.Finally, the mathematical tools for adapting specific organ models as modules of the whole-body model are presented.

The Modelling Framework
The whole-body model is a set of ordinary differential equations describing the dynamics of a large number of metabolite concentrations in many different, specific locations of the body.The whole-body model consists of dynamic modules, the organ models, which are connected metabolically by connecting compartments, typically the blood and the interstitial compartments.The modules are the anatomical entities, heart, liver, brain, kidney etc, but also more disperse structures like the erythrocytes are treated as an organ.The connecting compartments are assumed to be metabolically inert, although this may not be strictly correct.An example of how a whole-body model may be structured is shown in Figure 2(d).
Below we introduce the nomenclature, the mathemati-cal structure of the model, and the specific constraints on both level 0 (whole organism; cf. Figure 1) and level -1 (individual organs).

Definitions
For each of the organ models, the chemical amounts inside organ i is denoted by the vector m i.The time evolution of chemical amounts in organ model i is described by the set of ordinary differential equations, d d f .An example of the organs that could be described by f i and m i is given in Figure 2(a).
The organ models interact through the connecting compartments as illustrated by ellipses in Figure 2(b).The chemical amounts in these connecting compartments are denoted by the matrix n.Rows correspond to chemical substances, columns to compartments.Chemical substances can be everything from ions, over metabolites to signalling molecules; in the remainder of the article, we will refer to this broad class of substances as "metabolites".The dimension of n is M C  , where M is the number of metabolites that are present in the connecting compartments, and C the number of such compartments.The constant matrix t D  n has the same structure, and contains the experimentally determined rates of change of metabolite m in the connecting compartment c in the physiological state .
To keep track of the connections between organ models and connecting compartments, we introduce D i , the exchange matrix.This matrix is constructed from the identity matrix; diagonal elements, which correspond to compartments that have no physical contact to organ i, are changed to zeros.The matrix is finally augmented with a single row of zeros, to eliminate the storage subset of b (see below).These connections between organ models and connecting compartments are illustrated in Figure 2(c The matrix has the same structure as D i , but both reversible reactions and missing reactions have a zero diagonal element.If there is only inflow from a compartment, the diagonal element is changed to -1, if there is only outflow it is left unchanged.Example 1 in the discussion illustrates the definitions introduced above.

Model
A whole-body model takes the form of a set of ordinary differential equations, ODEs, describing both the time evolution of the chemical amounts of the metabolites being exchanged between the organs and/or the surroundings through the connecting compartments, as well as nternal conversions inside individual organs: i C Thus, a whole-body model is composed of a set of submodels, each describing the metabolic function of an organ.The whole-body model is obtained by describing the interactions between these organ models via the con-necting compartments.The dynamics of the body's overall metabolism is simulated by numerical integration of Equation 1.

Constraints
In order to insure that the whole-body model yields physiologically realistic simulations of a chosen physiological state , we establish the following constraint on the model at level 0: Copyright © 2012 SciRes.OJAppS In other words, the concentrations of the metabolites in the connecting compartments (the ellipses in Figure 2) should change at fixed rates that match the experimental observations in the physiological state .
In addition, it is required that the functions of the individual organ models match reality.This means that their input-output relations must fit the experimental determinations.These functional organ constraints are given by i.e. there should exist an m i such that the rates of metabolite transport and storage buildup match the measured values in the physiological state .An example of functional organ constraints is given in Section 3.4, Example 2.
In addition to these nontrivial constraints on the behaviour of the individual organ models, we impose the following technical constraints on i  b and

 
, i i b m n , which are based on blood flow direction and spatial organisation of the organs in the blood circulation.Uptake by an organ can only occur from the in-flowing blood and vice versa for the out-flowing blood, Also, no organ model should exchange metabolites with connecting compartments, with which it is not physically connected:

The Process of Model Development
The two main tasks of whole-body model construction are the experimental determination of functional organ constraints i  b , and the adaptation of organ models f i to these functional organ constraints.Once an organ model has been adapted to the functional organ constraints, it can be plugged into the corresponding whole-body model, and the whole-body model is operational when all the organ models specified in the whole-body model structure have been plugged in.
As indicated above, the functional organ constraints are calculated from simultaneous measurements of the blood flow and the concentration differences between arterial and venous blood across an organ in each of the specific physiological states .For each of these states, the resulting functional organ constraints must be internally consistent in terms of mass balances, i.e.Equation 2 should be fulfilled.
The adaptation of an organ model to measured metabolite fluxes is a non-trivial problem of constrained optimisation, where the constraints are given by Equation 3. We demonstrate the feasibility of this crucial step in the following section, where we give two explicit suggestions for adaptation of organ models to functional organ constraints.Our approach is related to that of [46] and [34].

Matching the Fluxes of an Organ Model to a Specific Physiological State
We propose a flux-matching strategy based on the formalism of extreme currents, which provides an efficient way of describing all possible steady state fluxes in a chemical reaction system.See chapter II, Section C in [47] for details.The kinetics of a chemical reaction system is generally given by d dt    m f ν v , where is the stoichiometric matrix of the reaction network, and is the velocity vector, which indicates the velocities of each of the reactions in the network.As before, f is the organ model and m is the vector denoting the chemical amounts inside the organ.(Note that the organ-denoting indices i have been left out, as we here deal with a single organ at a time.)The possible steady states of the reaction network can be found as the null space of v, i.e. the space of all vectors v for which .We will match the organ model to the flux matrix such a steady state; this steady-state constraint is more restrictive than the constraint imposed by Equation 4, since strict stationarity is required for all non-storage metabolites in the organ.
Geometrically, the part of velocity space where   ν v 0 is a polygonal cone (colloquially: an inverted pyramid of infinite height) extending from origo, and it is generally of lower dimension than the velocity space.This steady-state part of velocity space is known as the current cone, and the velocity vectors defining the edges of the cone are known as the extreme currents ("current" is short for steady-state velocity).Any steady state velocity vector can be represented by a point in the current cone and expressed as a "convex combination", or linear combination with non-negative coefficients j k , of the extreme currents E k : This representation need not be unique, however.
The functional organ constraints in the matrix  b only define the input-output relations of the organ model, so many different steady-state flux distributions v will match.For each reaction r in v matching an element of  b , the combined steady state and functional organ constraints are given by: Copyright © 2012 SciRes.OJAppS where , k r are the elements of E k .If the constraint cannot be met, then this implies that the reaction network (the "stoichiometry") of the organ model makes it impossible to have a steady state with the input-output relations imposed by the functional organ constraints, and the network must be changed accordingly.

E
The constraint is linear.If a quadratic optimisation criterion is given, then the efficient method of convex quadratic programming can be used to find a global optimum [48].The optimisation criterion would typically express how close the steady-state flux v of the modified organ model is to that of a particular steady-state flux w  in the original organ model, e.g.
. If K is included in the parameter optimisation, then this objective function defines the optimal steady-state as the one with the smallest relative flux changes compared to the original model, i.e. the flux distribution pattern is preserved while the flux magnitudes are allowed to change.If, on the other hand, K is fixed at 0 K  1 K  , then the optimal steady state is the one with the smallest absolute flux changes compared to the original model.Example 2 in the discussion illustrates the use of the flux-matching strategy presented above.

Fixing an Organ Model at Its Optimal Steady-State Flux
Once an optimal steady-state flux has been found, it must be ensured that the kinetics of the organ model allows it to operate at this steady-state flux.This can be done in several ways; we will now indicate a particularly simple one [46].The rate of a chemical reaction r can be written as where r is the velocity parameter of the reaction, K r is a vector containing the intrinsic parameters of the reaction (i.e.all remaining parameters), and c is the concentration vector given by m and n.The splitting of the parameters into velocity parameters and intrinsic parameters makes it straight forward to fix a model at a particular steady-state flux v and a particular set of concentrations c, which we want to become steadystate concentrations of the modified organ model.For each reaction, the velocity parameter is now given by Operationally, this corresponds to adjusting the enzyme activities ( max V values) of the organ model, so as to meet the fluxes of the optimal steady state.The use of the method is exemplified in 3.4, Example 2.

Adapting an Organ Model to Multiple Physiological States
For multiple physiological states, the problem is essentially the same as above, only now, instead of a single set of transport matrices, we use several sets of  b .In the different states, internal chemical amounts m i will in general be different, but parameters and K r are required to be the same between states.r  With several physiological states the problem can, in general, no longer be solved by adjusting only the velocity parameters r  , so the method must be expanded to include all parameters p and concentrations  c in the system.
As before, each state, , gives a set of constraints on the convex coefficients in the corresponding state: for all states,  , and elements of b, b r .In addition, the in- dividual rate expressions must match the flux-distribution: The optimisation would typically be done using the objective function where ,orig k are the reference parameters, are the reference concentrations, and p ,orig k c  par and met is a set of user-defined weights that determines the relative importance of the different classes of data.

C C
Equations 11 are now nonlinear, and this might render the optimisation problem computationally very demanding.It is, however, possible to reduce the optimisation problem to a polynomial programming problem if certain parameters are not included in the optimisation-in the typical case, Hill coefficients would need to be excluded from the optimisation, but max V parameters, K M parameters and  c vectors could be included.If the optimisation is performed in this way, the global optimum can always be found [49,50].

Discussions
We have developed a mathematical formalism that allows for realistic modelling of whole-body metabolism of any complexity.The formalism takes the form of a modelling framework, where models are established by a top-down approach designed to encompass detailed biochemical sub-models, which are typically constructed by a bottom-up silicon-cell like process [31][32][33][34].A key feature of our approach is the exploitation of the natural modular structure (individual organs) of the body, which communicates in terms of chemical mass flux via connecting compartments, primarily the blood.The mass flux between organs is considered in terms of uptake, release and storage of substrates and products of metabolism, including inorganic ions.Other humorally transmitted substances like hormones and other signaling molecules can also be included on the same level as other substances.This structure of the framework also provides a realistic means of collecting the necessary experimental data at the organ level as arterio-venous concentration differences.To complete the experimental data sets it is also necessary to have organ blood flow information in order to quantify the net mass balance across the organ as well as tissue biopsies to quantify the balance of organ mass storage.For humans such data already exists for skeletal muscle [51][52][53] and brain [54,55] However, for the human brain, biopsy data is normally not available in vivo, albeit they may be obtained in special cases [56].In addition, similar data sets may be generated by the use of non-radioactive, stable isotope technique [45,53,57].In animal studies it will be possible to obtain complete data sets for all major organ systems, using also radioactive tracer techniques, albeit multi-organ data from the same animal will be technically very demanding.
While it is mathematically straightforward to expand this modular approach to more detailed structures, like cell types within an organ, or indeed the organelles of the cell, the major challenge in this direction will be to collect the necessary data.In contrast to the methodology proposed here for whole organs, the experimental techniques for probing input-output relations at the cellular or subcellular level will be highly dependent on the system studied.Even if no data can be obtained, the individual submodules would of course still be subject to wholeorgan constraints.The modular nature of the approach would make it natural to adapt many of the modular techniques already in use in systems biology for constraining the individual submodules [12][13][14][15][16][17].

Limitations
One major assumption of the modelling framework as it is presented here is that the connecting compartments are biochemically inert.Biochemical processes do indeed occur in blood plasma, and some enzyme activities are present, but we expect that the overall metabolic activity is negligible compared to that of the organs.If this assumption proves wrong in special cases, we suggest that models of the unexplained metabolic activity can be described as a new organ model connected to the compartment in question.
While we have chosen to focus the description on well defined physiological states, a natural extension of the proposed strategy would be to collect metabolite concentration and flux data during transitions from one physiological state to another, like rest to exercise, fed to fasted or normal to diseased states.Experimentally, this would involve time series sampled at a suitable time resolution.This data can serve as additional input to the organ model adaptation step, for example by using a forcing approach to ensure a more realistic behaviour of the full model.
At this stage our framework does not include modelling of the blood circulation.However, during transitions between states, e.g.rest to exercise, there are major flow changes.Such changes could, however, be modelled based upon experimental monitoring of blood flow by ultrasound Doppler sonography [44].

The Process of Generating Whole-Body Models
The workflow associated with the construction of the full model, applying the formalism detailed above, is illustrated in Figure 3.As shown the assembly and formulation is associated with a number of work processes.We suggest some tools needed for these work processes in the methods section.In the supplementary material the use of the formalism is exemplified by steady state data sets for human metabolism in the fed and the fasted state, and the adaptation of existing silicon-cell organ models to these data sets is exemplified with an existing erythrocyte model [9,10].
The formats necessary for the collaborative process are described in detail in Supplementary material.In short, whole-body model structures are represented in custom XML files with a hierarchical structure of physiological states, organs, connecting compartments and transported or stored metabolites.Organ models are represented in SBML [58] with additional identifiers that denote the name of the organ, and list the reactions that correspond to transport in and out of connecting compartments.
To make the proposed whole-body modelling approach operational, and add-on package to the SBtoolbox2 [59], that can perform the processes indicated by circles in Figure 3 would solve the problem.The construction of whole-body models is, however, entirely implementation specific, and the standardised file formats will allow alternative tools to be constructed.

Example 1: A Simple Whole-Body Model Structure
To illustrate the structures of the constraint and flow direction matrices introduced above, we introduce the simple model of the Cori cycle given in Figure 4.The model keeps track of extracellular glucose and lactate as well as glycogen storage, in that order.The connecting compartments are ordered as {arterial plasma, venous plasma}.The organs are ordered as: {muscle, liver, IO}.
The exchange matrices for this system become: The flux direction matrices become: Copyright © 2012 SciRes.OJAppS The amount matrix for the system in a given state could be: where A and V denotes arterial and venous plasma, re-spectively.
The functional organ constraint matrix for an organ has the following structure: where t denotes the amount transported per time to or from a compartment, and s storage per time in the organ.The matrix function liver has the same structure, only with functions instead of specific numbers.J is the blood flow through organ i [63], denotes the net rate of uptake of the corresponding organ, and , and denote the concentrations in the arterial, venous, and portal vein connecting compartments, respectively.Concentrations, which are irrelevant for the functional organ constraints, are not included.AA is quantified as alanine, FA as C18H36O2, glycogen as glucose residues, KB as β-hydroxybutyrate, TG as C57H110O6, and protein as alanine residues.A total blood volume of 5 L is assumed, and the blood distribution is 20% arterial, 70% venous and 10% other [64].Abbreviations: AA, amino acids; FA, fatty acids; KB, ketone bodies; RBC, red blood cells (erythrocytes); TG, triglycerides.state is the 36 h starved state where glycogen has been used up, and the organism is living off its own energy reserves, including breakdown of triglycerides for energy production and the use of muscle protein as substrate for the necessary hepatic gluconeogenesis, fueling brain, erythrocytes and other tissues needing glucose.The estimates are constructed by considering the normal physiological functions of the organs in these two states [60][61][62], while also imposing a number of level 0 constraints.The major constraints are energy consumptions of 3000 kcal/24 h in the post-absorptive state and 1500 kcal/24 h in the starved state (this applies for a normal, young human) and the corresponding O 2 and CO 2 exchange rates in the lungs.We also impose mass balance and a whole-body steady-state criterium, Equation 2 in the main text with t .For simplicity, all protein is represented as poly-alanine, all amino acids as alanine, all carbohydrate as glucose and glycogen, all fatty acids as stearic acid and all triglyceride as glycerol-tri-stearyl.The energy consumption of the biosynthesis reactions included in the tables are not accounted for, but the overall oxidation of all the substrates involved, including the storage molecules, is accounted for by standard mass equations of oxidation, e.g.CH

Adaptation of an Organ Model to Functional
Organ Constraints We have applied the flux matching method and the model adjustment procedure to adjust the erythrocyte model described by Mulquiney and Kuchel [9] to the functional organ constraints described in Table 1.The erythrocyte model includes all major aspects of the erythrocyte's metabolism: the pentose phosphate pathway (PPP), glycolysis and the 2,3-bisphosphoglycerate shunt (2,3-BGP shunt).The reaction network consists of 59 net reactions and 63 metabolites and other substances.After fluxmatching, the model is left with four active net extreme currents.One of these is glycolysis, another is glycolysis where the phosphoglycerate kinase reaction is bypassed by the 2,3-BPG shunt.In addition, two net extreme currents related to the 2,3-BPG synthase/phosphatase complex are active, albeit with very low net fluxes.Solution of the quadratic programming problem (cf.Matching the fluxes of an organ model to a specific physiological state above) yields the optimal distribution of flux between these four extreme currents, and dictates that all other parts of the reaction network must have zero net flux.Biochemically, this corresponds to closing down the PPP because no CO 2 production and no pyruvate consumption is allowed by the functional organ constraint RBC .This results in a 6% relative decrease of the hexokinase flux, a relative flux increase of 2% -4% from phosphoglucoi- somerase to aldolase, and in less than 1% relative flux decrease in the lower part of glycolysis.The 2,3-BPG shunt has almost unaltered relative flux, and the relative increase of the ATPase flux is 1%.To fix the model at this optimal steady-state flux, all values are changed proportionally.max V list of IDs as well as the initial concentration in each connecting compartment.Finally, each species contains a list of concentrations in each compartment for each defined physiological state.Initial conditions are required only for the model instance, and the physiological state concentrations are required only for the model constraints.<listOfSpecies> <species localname="A"> <id source="http://HuMMproject.org/ns/wbm/v0/metabolites">Aex</id><initialconcentration compartment="Art" value="2.3"/><initialconcentration compartment="Ven" value="4.5"/></species> <species localname="B"> <id source="http://HuMMproject.org/ns/wbm/v0/metabolites">Bex</id><notes>Some notes for species B...</notes> <initialconcentration compartment="Art" value="1.2"/><initialconcentration compartment="Ven" value="3.4"/></species> </listOfSpecies>

Organ Models
Organ models are SBML files that can be integrated in a whole-body model.An example organ model file can be found on the project website.All entities in SBML files such as compartments, species, and whole models have an optional annotation field.The model itself has an annotation field denoting the organ being modelled, and the physiological states it fulfills.

Figure 1 .
Figure 1.Hierarchical organisation of the body.The number of distinct components increases exponentially when a description of the body is made at lower levels of organisation: celltypes organs 10 n n   ; organelles ).The definition of connectivity is purely technical in nature; it is only taken into account during the construction of the organ models.The actual mass flows between organ model i and the connecting compartments is contained in the matrix function b i .The elements of b i are the transport and storage subset of f i , i.e. a reordering of the elements corresponding to reactions crossing the boundary between 1 M C   organ and connecting compartment.The extra compartment column keeps track of the build up or depletion of storage compounds in the organ model.The constraints on the fluxes that the functions b i have to fulfil are given in i  b , of identical dimension. refers to a specific physiological state.Due to the anatomical organisation of in-and out-flow to an organ there can be no uptake from the venous compartment and, likewise, no excretion to the arterial compartment.Blood flow directions thus create another technical constraint to be taken into account during the construction of an organ model.This direction of flow is shown as arrow heads in Figure 2(d).The information is contained in , the flux direction matrix.

Figure 2 .
Figure 2. Structure of a whole-body model with classes of structural elements being added from a to d. Boxes indicate organ models (a)-(d), ellipses indicate connecting compartments (b)-(d), lines between boxes and ellipses indicate connections (c)-(d), and arrows indicates mass flow directions (d).The dynamics of each organ is described by a kinetic model, and these organ models are coupled via the connecting compartments.GI tract, gastrointestinal tract; Port.V., portal vein; RBC, red blood cells; WBC, white blood cells. i

Figure 3 .
Figure 3. Workflow for the construction of whole-body models.Ellipses represent work processes, squares information or data structures.Greyed shapes indicate application or researcher specific processes/data structures; white represent data structures that requires a fixed format for collaboration.

1 .
Functional Organ ConstraintsAn example of a more detailed whole-body model structure is shown in Figure5.Tables1 and 2show estimates of the corresponding organ input-output relations in two physiological states, i.e. the tables provide functional organ constraints ( i  b matrices) for construction of organ models that match the two physiological states.One of the states is the early post-absorptive state, where glycogen is being synthesized in liver and skeletal muscle, triglycerides are being stored in adipose tissue, and protein is being build up in skeletal muscle.The other

Figure 5 .
Figure 5.The whole-body model structure of Tables1 and 2. Organs 5 -10 = {muscle, fat, brain, heart, kidney, other}.Port V. is the portal vein, and the RBC organ is the erythrocytes.Table1.Functional organ constraints for the early post-absorptive phase ( = pap).

Table 2 :
For all organs except the erythrocytes (RBC), i  b matrices; these are presented in a condensed form in Table1and