A computational model of the glucagon/insulin-driven liver glucohomeostasis function, focusing on the buffering of glucose into glycogen, has been developed. The model exemplifies an ‘engineering’ approach to modelling in systems biology, and was produced by linking together seven component models of separate aspects of the physiology. The component models use a variety of modelling paradigms and degrees of simplification. Model parameters were determined by an iterative hybrid of fitting to high-scale physiological data, and determination from small-scale in vitro experiments or molecular biological techniques. The component models were not originally designed for inclusion within such a composite model, but were integrated, with modification, using our published modelling software and computational frameworks. This approach facilitates the development of large and complex composite models, although, inevitably, some compromises must be made when composing the individual models. Composite models of this form have not previously been demonstrated.
The construction of predictive models of complete biological systems that can be used to develop the understanding of integrated biological function is a key aim of systems biology. Biological modelling efforts have followed two general approaches. The first is simulation, where a usually restricted set of processes is modelled with fidelity and detail. This places heavy demands on biological data, in both quality and quantity . Alternatively, physiological behaviours are modelled with the aim of capturing the underlying control principles governing the behaviour of interest, with only key biological mechanisms characterized and parametrized in detail . The most comprehensive modelling effort to date—the Hunter/Noble model of the heart [3,4]—provides a hybrid between these two approaches.
Mathematical modelling of biological systems has been conspicuously successful in two particular circumstances: small models of specific, well-characterized system elements, e.g. in electrophysiology ; or larger models based on high-throughput data and extensive, thoroughly understood simplification methods, e.g. in metabolic flow analysis . Models that are both comprehensive and detailed, system-level models, covering a variety of interacting phenomena crossing many scales, are few , and those that exist have been extremely time-consuming and expensive in effort.
This paper presents the development of a model of glucose homeostasis, based principally on liver function, employing a novel ‘engineering’ approach for generating composite, multi-scale models. The approach allows for the construction of a flexible, large-scale system model by combining component sub-models that may be individually modified, and their number and/or complexity increased or decreased as appropriate. This allows some components to be characterized in greater detail than others, and has the flexibility to enable component substitution to refine characterizations when focus on the role of particular biological mechanisms is required.
Glucose homeostasis is an example of rein control . Insulin and glucagon, both generated in the pancreas, instruct the liver to, respectively, store and release glucose, thereby achieving glucose homeostasis. Insulin activates the storage of glucose in hepatocytes as a polymer, glycogen . Glucagon triggers the release of glucose from glycogen, which is broken down and built up by two enzymes, glycogen phosphorylase and glycogen synthase, with active and inactive forms differentiated by phosphorylation. The degree of activation of these enzymes is determined by the equilibrium between other enzymes, which also have inactive and active forms, in a series of phosphorylative epicycles . Insulin and glucagon affect these enzyme equilibria via intracellular second messengers including cyclic adenosine monophosphate (cAMP), calcium, inositol trisphosphate (IP3), protein kinase A (PKA) and others [10–13]. Insulin secretion by the pancreas is triggered at about 4 mM blood glucose, while glucagon secretion is inhibited at about 7 mM . In addition, liver glycogen storage is controlled directly by glucose and its phosphorylated forms , by the glucose-sensitive use of glucose in other organs of the body , and by non-pancreatic hormones (not modelled in this paper).
The principal pathology of glucose homeostasis is diabetes. A key test of models of the liver–pancreas interaction in achieving glucose homeostasis is type 2 diabetes, where the liver has reduced sensitivity to insulin. Models can be tested by reference to continuously measured blood glucose trajectories . One feature of such measurements is the presence of blood glucose oscillations with periods ranging from half an hour to 2 h, termed ultradian insulin oscillations (e.g. [18,19]). Much of the focus in trying to understand the genesis of these oscillations has been on pancreatic (insulin production) dynamics . However, the role of the liver has been of increasing interest in recent years  and modelling liver–pancreas feedback has a relatively long history .
Liu & Tang  and Liu et al.  have modelled glucose homeostasis focusing specifically on the dynamics of the insulin signalling pathway along with the action and kinetics of key enzymes involved in the process of glucose mobilization by glycogen phosphorylase and glycogen synthase in the liver. These authors construct a ‘flat’ model consisting of approximately 30 differential equations with 70 parameters. Their model is parametrized from the literature and from fitting to high-level data (also obtained from the literature). These authors' model performed well in some circumstances in qualitatively matching experimental insulin production data. The model could also produce stable oscillations known to occur in both blood glucose and insulin concentrations, although these did not qualitatively match some features of the empirically observed oscillations.
This paper shows how an engineering approach, which exploits modularity in a fundamental way, can be employed successfully to construct a detailed, multi-scale model of a major physiological system: the glucose homeostasis system. The model can be tuned and validated against existing data, and is capable of generating interesting and novel outputs, which potentially are of clinical interest. Our model is constructed to be amenable to further elaboration through refinement of existing components and the addition of new components, without requiring major reorganizations at each step. Problems in physiology usually involve many parts of the organism. Measurements taken either in vivo at the level of the organism or in vitro at the level of the cell inevitably are affected by connected systems and any potential treatments would also affect the whole organism. This integrated aspect of physiology makes it important to be able to model many interconnected systems.
Modularity has significant advantages for the construction of complex system models. (i) Modularity aids understanding by presenting a system in distinct functional units. (ii) Different researchers can work on each module separately, distributing effort and expertise. (iii) Modifications can be made to one module without affecting others. If desired, a module can be replaced entirely. (iv) Modules may be re-used as part of other projects and a library of models may be gradually accumulated. (v) Algorithms based on division of a system into elements based on differing spatial or temporal scales can be used to integrate modular components . (vi) Different components can be encoded in different computational environments, languages and tools depending on developer preferences. (vii) Modules may be independently validated or tested.
The incremental development of a single, gradually expanding model, with each refinement being made only when the existing elements have been well validated, is necessarily slow. The present approach allows models to be added together to create a complex composite model for tuning, validation and generation of results for analysis. The strategy deliberately makes use of models of individual elements that are already available in the literature. This has been made possible by a software framework [26–28] which allows models to be tied together in a model management system. The models are linked together after being wrapped using a composite model description language (CMDL), an XML-based language that can be used to specify composite models . The language supports specification of structure and implementation details for composite models, along with the interfaces provided by each sub-model. At the level of computational implementation, where it is inconvenient or costly to link distinct components into a single large ‘flat’ model, problems such as inconsistencies in structure or time scales between components can be managed by an adaptation of the numerical technique known as ‘waveform relaxation’, originally developed for parallel computation . This is a numerical technique for computing solutions to a system of ordinary differential equations (ODEs) that enables the integration of independent sub-systems (or modules) by using the outputs of sub-systems as inputs to others and vice versa. The algorithm makes it possible to integrate models of different mathematical forms independently of the solution technique that might be employed for individual components.
2. Building a composite model for glucose homeostasis
The composite model brings together a series of sub-models, some taken from or based on established models and others developed ab initio, of aspects of the physiology of glucose homeostasis. The current choice of component models is determined by the minimal physiology that needs to be represented to describe the basic features of the glucose homeostasis system and key elements of its control. Some sub-models represent the physiology at a greater level of detail than others.
Seven component models describe the major elements of glucose homeostasis, with multiple feedback relationships between sub-models, and an overarching feedback cycle: liver → blood → pancreas → liver. This cycle is moderated via the blood by dietary intake of glucose and active consumption of glucose in organism activity. The model aims to mimic quantitatively the behaviour of glucose regulation following external stimuli to the liver both in a healthy state and with diabetes. The interchange between glucose and glycogen is influenced by four factors: glucagon via calcium, glucagon via cAMP, insulin and directly by cellular glucose levels. These factors impinge upon a central model that represents the computation carried out by competing enzyme phosphorylation equilibria. These nested feedback structures, and the inputs and outputs of each component model, are represented in figure 1.
Five of the seven models represent aspects of hepatocyte physiology. The models treat the whole body of hepatocytes across the liver as a single cellular component, and the parameters reflect this. Complexities introduced by metabolic zonation  and inter-hepatocyte communication through direct cell–cell communication channels will be considered in future models. The two remaining components are a simple model of the pancreas and a model mediating between pancreas and liver representing blood transport.
Each component model is first described independently, together with, for the more complex models, some example results obtained when the model is considered separately. This allows independent verification of the component's behaviour in a context which does not depend on its embedding within the larger system. The mathematical descriptions of each sub-model and its parameters are given in the electronic supplementary material (appendix A) and all models are available at http://www.compbio.org./models.html. A summary of the models is given in table 1.
2.1. Glucagon receptor model (A)
This module represents the activation of a G-protein-coupled receptor by a hormone stimulus. The G-protein-coupled receptor is an important element of many signalling pathways, and this module can, therefore, be used as a component of several signalling systems. Receptor activation of the G-protein causes GDP for GTP exchange, which results in conformational changes in the G-protein and separation of the Gα sub-unit from Gβγ sub-units. Both GTP-bound Gα and Gβγ are able to regulate effectors [30,31].
Activation of the α sub-unit is considered under the cAMP model. The glucagon receptor model describes activation of the βγ sub-units of the G-protein and is based on Nauroschat & an der Heiden  and Riccobene et al. . We have added the known effect that calcium increases receptor inactivation . The processes modelled are: ligand–receptor binding and dissociation, receptor sequestration and desequestration and its dependence on receptor phosphorylation state, receptor phosphorylation and its dependence on active G-protein and ligand-binding, G-protein activation and inactivation and its dependence on calcium and phospholipase C (PLC), and the production of PLC by active G-protein.
Model results in figure 2 show the characteristic initial rise then fall seen in G-protein responses to glucagon stimulation (inset), which compares well with the results in Bridgette et al. . A sustained stimulus produces a transient response. A second glucagon stimulus reveals desensitization of the G-protein receptor .
2.2. Calcium model (B)
This describes the calcium signalling pathway activated by IP3, taken from Höfer . It is a simplification of Höfer's model in which all switch-like response curves are unified as a Hill function with a common coefficient, as discussed in Hetherington et al. . The model includes IP3-dependent calcium entry through the cell membrane, calcium- and IP3- dependent release of calcium by the endoplasmic reticulum (ER), and the ER and membrane calcium pumps. The parameters were obtained from Höfer  as explained in Hetherington et al. .
Höfer's model was developed to describe the calcium oscillations induced by vasopressin and epinephrine. However, glucagon generates only a single calcium transient . Figure 3a shows a single calcium transient recorded experimentally, while figure 3b shows results generated by linking the glucagon receptor model to the calcium model and illustrates the transient IP3 response generated by 5 nM glucagon together with the consequent single calcium transient.
The G-protein and calcium models are coupled by PLC/IP3 and calcium. The G-protein model produces an output of PLC that generates IP3, the input to the calcium model. The scaling between PLC and IP3 is assumed to be linear and is determined by the scaling parameter GSPLC, a component of the composite model (see electronic supplementary material, appendix B). This was set as 1 µm active PLC generates 100 µm IP3, so that, given the range of values of PLC produced by the receptor model, the range of values of the IP3 variable is as expected in calcium oscillation modelling [34,36]. The calcium output generated by the calcium model is fed back as an input into the G-protein model.
2.3. Cyclic adenosine monophosphate model (C)
This model, developed ab initio, represents activation of the G-protein α sub-unit of the glucagon receptor model (A), which has different dynamics from the βγ sub-unit associated with the calcium pathway described above. Key modelled processes are the activation of the receptor, the production of cAMP by the active receptor complex, the elevated rate of cAMP production by activated receptor, the activation of PKA and the potential nuclear localization of PKA, which will become important when the current composite model is extended to cover the transcriptional changes that result from hormonal stimulation of hepatocytes. Like the calcium component, the cAMP model is based largely on Hill function response dynamics. Receptor–ligand binding dynamics are assumed fast, so that the ligand-bound receptor is not a separate quantity, but rather obtained from the total receptor population by multiplication by an appropriate function, L(t) (electronic supplementary material, equation A 3.3). Parameter values were based on inspection of our experimental data: an example of experimental measurements of cAMP production in hepatocytes in response to glucagon is shown in figure 4. Details of the methods used in our experiments are given in the electronic supplementary material (appendix D).
2.4. Insulin model (D)
The current composite model uses a very simple, ab initio model for the insulin receptor and associated pathway. This model declares that the controlling species that tends to increase the amount of active (α-form) glycogen synthase kinase (GSK) has a threshold response to insulin with a relaxation time delay. The representative variable for this pathway's effects on glycogen is the proportion of inactivated GSK-3. The inactive β-form is the phosphorylated form, so that an increase in the GSK variable results in increased glycogen synthesis activity.
2.5. Blood model (E)
The blood model was developed ab initio to describe the movements of glucose between the blood, the liver and the pancreas. Blood glucose is transported into the hepatocyte, where it is converted into glucose-6-phosphate. The membrane glucose transporter (GLUT) is passive, but selective, transporting only glucose, and not glucose-6-phosphate. As a result, the transport appears to be active, with glucose apparently pumped into the cell with a slower transport time constant than might be expected for a strictly passive GLUT. In this and the glycogenolysis models (F—see §2.6), a single intracellular glucose pseudo-concentration variable represents the effective total glucose and glucose-6-phosphate. This illusory-pump model fits the data for a perfusive radiolabelling experiment in pigs in Munk et al. , where it was found necessary to consider separately a forward clearance rate into cells (the effective pumping) and a reverse rate out of the cells. The apparent active pumping of glucose actually represents the glucose to glucose-6-phosphate equilibrium—the GLUT pump itself is known to be an effective rapid passive transporter. The radiolabelling methods of Munk et al.  cannot distinguish between intracellular glucose and glucose held in the cell as glycogen or in the phosphorylated form so that an apparent increase in cellular glucose pseudo-concentration above the concentration in the surrounding medium will reflect glucose held in other forms. The advective effects of blood circulation are neglected—the blood is treated as a single effective volume, in contact with liver and pancreas and mediates communication between the two via the blood glucose concentration variable gB (see electronic supplementary material, appendix A).
An important input to this model, which is also an input to the composite model, represents feeding, exercise and the basic metabolic drain of the (unmodelled) rest of the body. This exogenous rate of input/output of blood glucose is denoted by M. We define M = 0 to be zero net glucose supply/demand, so that M = 0 can be considered to be a situation where a person is given a glucose drip that exactly matches the constant metabolic demand of the resting state. An increase in M represents an increase in blood glucose, as, for example, after feeding. A fall in M represents a fall in blood glucose, as, for example, during fasting. This device allows us to avoid representing the rest of the body as a separate component.
The model involves a quantity, denoted vf, that represents the effective volume ratio of liver cytoplasm to blood volume, so that a unit concentration fall of cellular glucose in the liver produces vf units of blood glucose concentration increase. The transmembrane glucose transport parameters are based on the values in Munk et al. , where the steady-state ratio of cellular pseudo-concentration to blood concentration is somewhat over unity, and the time scale of response (k2 in the study of Munk et al. ) is of the order of 1 min. The selection of appropriate values for vf is explained below as this choice is complicated by the process of model composition.
2.6. Glycogenolysis model (F)
Many factors control glycogen breakdown and synthesis. This model incorporates direct control by cellular glucose and glucose-6-phosphate, based on the model in Cardenas & Goldbeter , control by calcium ions, based on the model in Gall et al. , control by cAMP via PKA, and by insulin via GSK, based on the representations in Bollen et al.  and King . The model includes direct control of glucose by glucose, in addition to pancreatic hormone control. In order to reconcile these four competing inputs, one could, in principle, model the full enzyme kinetic scheme consisting of glycogen phosphorylase and synthase, their phosphatases and kinases, and the additional phosphatases and kinases involved. This phosphorylation cycle cascade is, however, extremely complicated and not sufficiently well characterized experimentally for detailed dynamical modelling to be possible. Instead, a simplified, pragmatic approach has been chosen. The overall decision system, whereby the protein cycles that compute the behaviour in terms of the cell's second messengers, rather than being modelled in detail, is represented using fuzzy logic . A Boolean logical expression, based on each of the four inputs, determines whether glycogen build-up and breakdown are on or off. This Boolean expression is extended to the domain where the four inputs are not themselves Boolean variables, but, rather, intermediate degrees of cellular activity in the unit interval 0–1. The model is completed by defining a simple saturating function to represent the saturating kinetics of the enzymes, and specifying response times for these activities to the inputs.
In developing the model, the logical expressions were explored to achieve a model consistent with known behaviour. Parameter values for this model were based on Gall et al.  and Cardenas & Goldbeter , obtained from the enzyme kinetics database BRENDA .
2.7. Pancreas model (G)
To complete the glucohomeostatic feedback cycle, we need a pancreas model responsible for the production of the hormones glucagon and insulin. In order to retain a focus on hepatic physiology, we developed ab initio a particularly simple pancreas model, in which the release of glucagon or insulin follows time-delayed threshold responses. This is modelled using a logarithmic response to the ratio of blood glucose to a fixed reference level, defined so that there is a zone within which the rate of glucagon or insulin release is small, but which rises steeply as the perimeter of this zone is approached. An additional factor not modelled directly is the timescale of the circulation between pancreas and liver, which is represented using only the pancreatic and liver response times.
Figure 5 shows the consequences of an imposed activity cycle obtained by coupling this simple pancreas model to a similarly simple liver model, incorporating the actions of glycogen phosphorylase, assumed here to be controlled only by glucagon, and glycogen synthase, assumed here to be controlled only by insulin. Positive activity values represent glucose comsumption by the body and negative values glucose intake from feeding. The example assumes two cycles in a 24 h period.
Oscillations in blood glucose concentration about the reference level at the system's natural frequency are induced by the periodic activity-feeding (glucose demand–supply) function M(t). However, additional, high-frequency oscillations are superimposed on this cycle. In the full composite model, these will be identified with the ultradian glucose oscillations . It is interesting that these oscillations appear even in this very coarse-grained model as an emergent property of pancreas–liver feedback dynamics. The frequency of the oscillations observed here and in early runs of the composite model is faster than observed experimentally for ultradian oscillations. This arises from values for enzyme kinetics taken from the BRENDA database, which are probably inappropriate (see §4).
3. Model tuning and validation
3.1. Tuning the model parameters
The seven models described above must now be composed together, matching each input to the appropriate output, as indicated in figure 1. The scalings of time and of the transferred variables must be made consistent, and, for the most part, follow from the units indicated above for each component model (table 1). Scaling factors introduce parameters to the composite model that are not present in any of the original component models. For example, one model may use the number of molecules and another molarity. A volume parameter must be introduced to produce appropriate scaling between such component models. Where the additional parameter occurs only through multiplying combinations of other parameters it can, if desired, be removed by a reparametrization of the composite model. The choice is either to have additional parameters or reduced parameters that are less easily identified with biological properties. Details are given in the electronic supplementary material (appendix B).
Model tuning was done manually: first using the most appropriate published data and then modifying parameters, while maintaining consistency where parameters were not independent, to fit with available experimental data. We began with parameters selected for individual models. Details of the procedure used for tuning the composite model are given in the electronic supplementary material (appendix C).
In initial runs of the model using the preliminary parameter estimates, the liver was cleared of glucose extremely rapidly, with ultradian glucose oscillations displaying an excessively high frequency, but with the relative amplitude of the first peak to other peaks of the oscillations, and shape, being broadly correct (figure 6).
As discussed in the electronic supplementary material (appendix B) competition between components and mixing of enzyme and substrate will be slowed by diffusion in the intracellular space. This was accommodated by slowing down the relevant processes by reducing rate constants or, equivalently, increasing characteristic times (reciprocals of rate constants). The model was then shown to fit published data well (figure 7).
The tuning was also checked by comparison of time scales with direct experiment in which the rate of glucose release from cultured hepatocytes was determined. Figure 8 shows the release of glucose from primary, cultured hepatocytes that had been loaded with glucose by the addition of insulin to glucose-containing culture medium overnight. At time zero, the bathing medium was replaced by medium lacking glucose and the glucose released into the medium monitored over the next 8 h in the absence and presence of glucagon. Glucose release was relatively slow and occurred with a time course consistent with the model using the adjusted parameter values. The production of cAMP by glucagon-treated hepatocytes referred to above and illustrated in figure 4 similarly displayed a relatively slow time course, consistent with the model using the adjusted parameter values.
3.2. Model validation
Individual models were fitted from small-scale data and the degree of fit of the composite model compared with the behaviour of the overall glucose homeostasis system. Parameters were chosen for each component model and the models validated individually against appropriate experimental observations. Scaling parameters needed for model composition were identified from the definition of the models' units and appropriate data where possible. Where parameters remained to be determined, these were tuned as part of the validation of the composite model, based on higher scale observations of the whole system as described above. The component models' parameters were not further tuned to obtain an optimal fit to experimental results for the whole system. We have required only that the model gives a reasonable description of the key qualitative features of the whole system and shows order-of-magnitude quantitative fidelity.
Validation was based on the following criteria:
— existence of a homeostatic response to a range of glucose challenges;
— existence and correct behaviour of ultradian oscillations under appropriate conditions; and
— match to observed experimental behaviour in response to a glucagon challenge .
A fundamental validation step is to examine glucose variation in the presence of a challenge to homeostasis. This was achieved by either inserting (feeding) or removing (exercise) excess glucose from the blood (represented by the activity function M(t) in the blood glucose model (E)—see electronic supplementary material, appendix A). This approximates to an oral glucose tolerance test (OGTT), although, for modelling convenience, a sudden-onset, subsequently sustained glucose input is used: M(t) = M, a constant (mM s−1), for t ≥ τ0, the time of onset. This matches the conditions of continuous enteral nutrition employed by Simon et al.  and figure 7 shows that the tuned model (figure 7b) qualitatively reproduces their observations (figure 7a) reasonably well.
Since cAMP control of glycogen phosphorylase via the α sub-unit of the G-protein is probably dominant in hepatocytes [42,43], in order to restrict the number of potential parameters and variables to be explored, we have inactivated the G-protein βγ sub-unit (A) and calcium (B) modules for the results presented below, removing that entire arm from the model structure (figure 1)—see electronic supplementary material (appendix A, equations (A 6.9) and (A 6.10)). The modular structure of the composite model makes such adjustment extremely simple.
Figure 9 explores the behaviour of the tuned model to different challenges. The first is the consequence of a small glucose input (M = 5) or demand (M = −5), which produce a raised or decreased steady-state glucose concentration and a steadily rising or falling glycogen level (figure 9a,b). There is a small shift in the blood glucose level, because the homeostatic machinery stimulates glycogen synthesis or breakdown, just enough to cope with the glucose supply or demand. For a negative glucose input, reflecting an increase in demand or reduced glucose supply (figure 9b), glucagon rises to mobilize glucose from glycogen stores and blood glucose first falls to a new, lower stable state. This is maintained until the glycogen store is exhausted and glucose homeostasis then fails, with blood glucose falling away to life-threatening levels, mimicking extreme glucose deficit.
For high levels of glucose input (M =10), blood glucose can be maintained at physiological levels provided that the enzyme GSK retains a normal sensitivity to insulin. Figure 9c shows that even if GSK resistance to insulin is increased (from a normal value of tI = 0.5 to tI = 0.8 in model D; electronic supplementary material, appendix A 4) normal glucose homeostasis can still be maintained by the direct glucose control mechanism because a rise in cellular glucose concentration stimulates glycogen synthesis despite the reduced sensitivity to insulin that prevents a fully hormone-driven response. While this could be another indication of the robustness of the system, it nevertheless suggests that small reductions in insulin sensitivity, such as might occur early in the diabetic response, can be accommodated through the direct glucose control mechanism and therefore could be masked, potentially delaying diagnosis of type 2 diabetes.
Figure 9d shows that if GSK resistance to insulin is lower than normal, so that insulin sensitivity increases, then, for M = 10 and an insulin sensitivity of 0.1, all parameters oscillate. By the time the liver responds to changes in pancreatic hormones, the blood glucose level has been driven too far in the opposite direction and the pancreas response remains out of phase with the liver response, producing varying glucose levels. In extreme cases, large oscillations result in acute temporary hypoglycaemia. Once the glucose concentration drops to such low values, other physiological factors not modelled would play an important role, so the subsequent course of the model after the first extreme oscillation is not valid.
It is of interest to consider other challenges to the glucose homeostasis system, in particular a glucagon challenge, as proposed in Lockton & Poucher . This provides a useful test of model reliability, since this virtual experiment was not performed until after the model content and parameters were settled. Figure 10 compares results from the model with those obtained by Lockton & Poucher , who challenged fasting subjects with a bolus of glucagon and monitored the rise in blood glucose, taken to be the consequence of glycogenolysis. The model produces a transient increase in blood glucose (figure 10a) that matches their results (figure 10b). They note also that the initial glucagon challenge is accompanied by a transient rise in insulin. The model (figure 10a) also shows a transient rise in insulin after an initial glucagon challenge. Additionally, the model reproduces the transient hypoglycaemia following glucagon challenge noted in Lockton & Poucher's experiments, although the time scale of the model response is about twice as fast. The results of Simon & Brandenberger  make it clear that there is substantial variation in time scales from individual to individual. However, the discrepancy in time scales between model and experiment is most probably due to the buffering effect of other tissues, not modelled. Detailed additional tuning of model parameters does not seem warranted given the satisfactory degree of qualitative similarity in the responses of model and experiment.
This paper presents an effective approach to building large-scale system models that allow sub-models to be composed, using existing models where possible, while ensuring that model couplings are syntactically correct as well as physiologically sensible, to create a complex composite model for tuning, validation and generation of results for analysis. The composite model of glucose homeostasis presented here builds on the modelling efforts of previous workers. The model crosses biological scales from individual molecules up to whole organism phenomena. As a result of this work, a number of issues arise.
Modularization requires the biology to be broken up into components. The restrictions of experimental work, imposed by technical limitations, can effectively split a system into components so that biologists rarely analyse a complete system. However, precise division of the biology into non-overlapping components is not always possible either experimentally or in models. This means that some modules may generate duplication in the sense that the same process participates in more than one component model, and these processes may need to be instantiated separately in the different model components. This is most likely when using previously published models to construct a composite model. However, for both modelling and physiological reasons, this may prove justifiable. For example, it may be reasonable to use a relatively simple model to obtain inputs for a module that it affects only in some minor features, and a more detailed model of the same component to obtain inputs to a different module that is affected significantly. This strategy may provide advantages in terms of computational time. Apparent duplication also may be necessary when common molecules are used in different biological pathways (modelled separately), but which do not result in cross talk owing to physical separation, or also potentially to scaffold proteins which contribute to the structuring of the intracellular environment. For example, one model's input function may be a function of another model's state variables. A typical example for ODE modelling is where species A exists in two compartments, each with its own model, and the rate of flux depends on the concentration of A in both compartments. One can express the function for flux as a separate term in each model, and pass the concentration in each compartment to the model of the other in an iterative process . However, it is essential to maintain consistent parameter choices.
In this study, the assignment of values to model parameters revealed the importance of using data acquired under experimental conditions appropriate to the model. Parameter values need to be determined under realistic, rather than idealized, experimental conditions. It is important to use effective rate constants that incorporate delays owing, for example, to diffusion of molecules through the cytoplasm and competition for substrates. Values deposited in databases such as BRENDA are usually derived from measurements made on purified enzymes under ‘ideal’ (supersaturated) conditions, in a stop flow apparatus where mixing is virtually instantaneous. These values are more than likely to be unphysiological.
The primary advantage of an integrative approach to building complex system models lies in the reuse of existing modelling work, which speeds the process of model building and exploits existing intellectual capital. With such a ‘model-based engineering’ approach, which scales modelling techniques to the substantial challenges faced, the full potential of multi-scale modelling in systems biology can be realized. This approach brings its own problems, as detailed above. However, the findings of our model make it clear that composite models constructed in this way can generate interesting biological hypotheses, potential experiments or theoretical advances at a much better effort-to-benefit ratio than would be possible with ab initio detailed modelling. The detailed experimental data that are necessary to underpin simulation are usually generated by focusing on particular elements of a process, according to the experimental techniques available. It is rarely possible to analyse large numbers of variables and parameters simultaneously, and obtain a uniform level of quantitative detail at all relevant levels. Thus, a pragmatic, modular integrative approach that allows components to be individually brought to a quality suitable for use in larger constructs is advantageous.
Once insights have been obtained, it may be possible (and instructive) to design simplified and abstracted models to obtain these features. However, such models will then be ‘stand-alone’ and will inevitably lose many of the subtle interplays between sub-systems that a ‘systems biology’ approach to physiology seeks to capture. This latter aspect will almost certainly be of critical importance when these models are extended and refined enough to be of use in a diagnostic clinical setting (a key long-term goal). Simplified and abstracted models will not have the comprehensiveness and flexibility for this.
The model itself demonstrates direct glucose control of the system, i.e. apancreatic stabilization of blood glucose, illustrating the general concept of robustness through redundancy in evolutionary theory. Functional glucohomeostasis has been observed after pancreatectomy [42,43]. If a suitable experimental model can be found, investigation of whether this is due to direct glucose control would be of value.
In the companion paper , the behaviour of the composite model is explored more fully.
This work was funded as a Beacon Project of the UK Department of Trade and Industry (DTI). We thank the DTI for their support. T. S. was funded by a CoMPLEX EPSRC Doctoral Training Centre Studentship. We have enjoyed useful discussions with Jonathan Ashmore and Rajiv Jalan.
- Received March 15, 2011.
- Accepted May 26, 2011.
- This journal is © 2011 The Royal Society