A composite computational model of liver glucose homeostasis. II. Exploring system behaviour

T. Sumner, J. Hetherington, R. M. Seymour, L. Li, M. Varela Rey, S. Yamaji, P. Saffrey, O. Margoninski, I. D. L. Bogle, A. Finkelstein, A. Warner


Using a composite model of the glucose homeostasis system, consisting of seven interconnected submodels, we enumerate the possible behaviours of the model in response to variation of liver insulin sensitivity and dietary glucose variability. The model can reproduce published experimental manipulations of the glucose homeostasis system and clearly illustrates several important properties of glucose homeostasis—boundedness in model parameters of the region of efficient homeostasis, existence of an insulin sensitivity that allows effective homeostatic control and the importance of transient and oscillatory behaviour in characterizing homeostatic failure. Bifurcation analysis shows that the appearance of a stable limit cycle can be identified.

1. Introduction

Complex biological systems involve many interacting phenomena which often need to be explored together to understand system behaviour. In a companion paper [1], we presented a composite model of glucose homeostasis using seven distinct validated models combined together in a model management system [2]. The power of this approach is that it permits the use of ordinary differential equation models from a wide range of sources and allows models to be re-used and replaced easily to test alternative hypotheses. There have been few comprehensive system level models of physiological systems covering many scales [3]. However, the ability to tie biological phenomena to system-level behaviour—particularly in clinical applications where most information is obtained from various parts of the system—will be important in the future.

The seven component models for glucose homeostasis which constitute the composite model presented in Hetherington et al. [1] are as follows: (i) glucagon receptor model for the activation of a G-protein-coupled receptor by a hormone stimulus, (ii) calcium model for the calcium signalling pathway activated by inositol triphosphate (IP3), (iii) cyclic AMP model representing the activation of the G-protein α-subunit of the glycogen receptor model, (iv) insulin model for the insulin receptor and associated pathway, (v) blood model to describe the movements of glucose between the blood, the liver and the pancreas, (vi) glycogenolysis model representing glycogen breakdown and synthesis, and (vii) pancreas model for the production of the hormones glucagon and insulin. The detailed models are given in the appendix A of Hetherington et al. [1]. All models are available at http://www.compbio.org./models.html.

The model was validated against data for the homeostatic response to a range of glucose challenges, to demonstrate the existence and correct behaviour of ultradian oscillations under appropriate conditions and to match observed experimental behaviour in response to a glucagon challenge [4].

Continuous measurement of glucose has demonstrated the presence of ultradian blood glucose oscillations [5,6]. Although originally considered to have their genesis in the pancreas, there has been increasing interest in the role of the liver. In the companion paper [1], we used the composite model to confirm that this phenomenon occurs as a result of feedback between pancreas and liver, first identified by Sturis et al. [7]. However, it has not been possible until now to undertake a more comprehensive analysis of the oscillatory behaviour of the composite system and this is developed in this paper.

Specifically, here we explore the behaviour of the composite model in response to a range of glucose challenges. We focus on insulin sensitivity because this is the variable that is most susceptible to lifestyle changes and most affected by diet. It therefore has the biggest impact on human health. In particular, we seek to ascertain the quantitative behaviour of the expected ultradian oscillations to demonstrate the power of using systems biology models.

2. Exploring system behaviour

In order to explore the system's oscillatory behaviour, the qualitative and quantitative behaviour of the blood glucose output was explored in the plane of parameter space described by the external glucose stimulus (M), and the threshold value for the Hill function describing the hepatocyte insulin response (tI). M is the glucose input to the blood cell glucose transport model (model E of Hetherington et al. [1]), and tI models the resistance of action by hepatocytes to insulin in the insulin model (model D of Hetherington et al. [1]). The correct response of the liver to insulin is important in the regulation of blood glucose. In particular, insulin resistance in which the response of the liver to insulin is below the normal level is a major cause of type 2 diabetes that accounts for over 90 per cent of diabetes cases globally [8].

In all cases, M(t) is modelled as a step-function, describing a sudden-onset, subsequently sustained glucose input. Also, a step function is the normal test for dynamic systems to demonstrate significant effects in a consistent manner. Instead of insulin resistance, it is more useful to refer to insulin sensitivity that is given by 1 − tI and this value is used in the graphs presented in this paper. Figure 1 shows three example outputs of the model for different values of tI: stable fixed-point solution, oscillatory solution and stable fixed-point solution reached via a transient overshoot. Our previous paper [1] gave a series of specific time course results using the model: fig. 6 showed the formation of the limit cycle, fig. 9 time courses for various values of glucose demand and fig. 10 response to a glucagon challenge.

Figure 1.

Sample model responses showing stable response (M = 7.5, tI = 0.6), transient with damped oscillations (M = 7.5, tI = 0.5) and stable limit cycle (M = 7.5, tI = 0.4).

The parameter tI is a derived one, so it is instructive to determine how it relates to normal and diabetic states. In most studies, insulin resistance was inferred from insulin or peptide measurements and was not directly measured [9]. Weyer et al. [9] give the following ranges for insulin: fasting healthy Caucasians (23 ± 9 μunit l−1), fasting diabetic (44 ± 21 μunit l−1); and 2 h after a meal healthy (100 ± 98 μunit l−1) and diabetic (242 ± 204 μunit l−1; a unit is 45 μg of insulin). We have assumed that our units of insulin are at the maximum normal levels produced by the pancreas (see electronic supplementary material of Hetherington et al. [1]) and have assumed that the 2 h data for a diabetic patient would be approximately equal to this. For our model, a value of tI of 1 would be approximately 240 μunit l−1 and tI of 0.5 approximately 120 μunit l−1 which is somewhat above the normal level.

A bifurcation analysis was undertaken to explore the changes in the behaviour of the composite model as M and tI were varied. The analysis was performed using the bifurcation analysis package, auto [10].

The bifurcation diagram in figure 2 shows the stable solutions of the model with respect to tI for three values of M = 0.5, 5 and 7.5 μM of blood glucose (as noted earlier the x-axis shows (1tI), such that increasing values indicate increased sensitivity to insulin). For low values of insulin sensitivity, there is a stable fixed point of blood glucose, the value of which depends on the value of M. As the sensitivity to insulin is increased, the blood glucose output undergoes a Hopf bifurcation at which an oscillatory solution (a stable limit cycle) emerges from the stable fixed point. As the insulin sensitivity is increased further, the mean blood glucose level (time-averaged mean) decreases but with an associated increase in the amplitude of the ultradian oscillations. Insulin sensitivity, therefore, has a significant effect on the stable blood glucose level resulting from a sustained glucose input, with insulin acting to lower that level. This, however, is at the cost of transient glucose (ultradian) oscillations.

Figure 2.

Bifurcation plot as tI is varied for three values of M showing when Hopf bifurcations occur. Solid line indicates a stable fixed point; filled circles indicate the minimum and maximum values of the limit cycle. Insulin sensitivity is plotted as 1 − tI.

The value of insulin sensitivity at which the Hopf bifurcation occurs also varies with the level of the external glucose challenge, M. For higher glucose inputs, the onset of the limit cycle occurs at lower values of insulin sensitivity. The larger the glucose input the greater the production of insulin by the pancreas and the less sensitive the liver needs to be to observe ultradian oscillations.

Figure 3 shows that at low insulin sensitivity (high insulin resistance), the level of sustained glucose input has a significant effect on the stable blood glucose level. However, at high sensitivities (low resistance), the mean level (time-averaged mean) is not significantly affected but transient glucose (ultradian) oscillations occur particularly at the higher sensitivities (figure 2).

Figure 3.

Final resting blood glucose concentration as a function of relative insulin resistance. For oscillatory solutions, the mean value is shown. The curves represent various glucose inputs, from bottom to top.

When the glucose input (M) is driven greater than 22.5 μM, in the model, the blood glucose level increases indefinitely. This does not occur in vivo where glucose in excess of that which can be converted into glycogen and stored is broken down to provide substrates for fat metabolism. The inclusion of the way in which excess glucose is handled is perhaps the most obvious extension to the composite model.

The model demonstrates the large transient first excursion in the ultradian oscillation, a feature observed in some of the available experimental studies [7]. Most studies, however, concentrate on the long-term behaviour under continuous glucose stimulus. Because the magnitude and duration of the transients vary with the model parameters, experimentally observed transients can provide information about the underlying physiology. This is illustrated in figure 4a,b. Fig. 9 of Hetherington et al. [1] shows example time courses for M = 10 tI = 0.8 (low insulin sensitivity resulting in no oscillations) and M = 10 and tI = 0.1 (high sensitivity resulting in oscillations).

Figure 4.

(a) Maximum glucose value reached as a function of relative insulin resistance. (b) Amplitude of initial transient relative to steady-state value defined as (xmax − xfinal)/(xfinal). M = 7.5, 10, 15 and 20. Glucose input values and code for lines as in figure 3 .

The maximum glucose level reached during the transient response to external glucose is controlled by the glucose input, and is largely independent of insulin sensitivity (figure 4a). However, the relative amplitude of the initial rise in blood glucose compared with the steady-state oscillations (figure 4b) does depend on the sensitivity of the liver to insulin. A high initial transient can therefore be used to identify insulin insensitivity that may have consequences for diagnosis and treatment of disease.

From these results, some general findings may be summarized as follows (figure 5 which gives a schematic of the qualitative behaviour of the model with respect to variation in the glucose input parameter, M, and the insulin sensitivity parameter, tI):

Figure 5.

The qualitative behaviours of the model as a function of insulin resistance and glucose input. When insulin resistance is high, homeostasis can be maintained by direct glucose control. As insulin resistance falls, homeostatic control is increasingly dominated by pancreatic release of insulin and glucagon, oscillations in blood glucose become increasingly prominent, and occur at progressively lower blood glucose levels.

Small glucose inputs that do not take blood glucose above the threshold for insulin release from the pancreas can be accommodated by direct glucose control, shown by the light grey region in figure 5.

For intermediate glucose inputs, even for a liver which is completely resistant to insulin (and hence not subject to pancreatic control), as found in severe type 2 diabetes or in type 1 diabetes, where the pancreas fails to release insulin, direct glucose control always results in glycogen synthesis to stabilize blood glucose levels (figure 4a with insulin sensitivity at or near zero). With extreme hyperglycaemia, it is known that blood glucose levels may be as high as 20 mM [11]. Abnormally, high blood glucose levels are diagnosed for blood glucose concentrations greater than 7 mM and restoration of a steady state after a glucose bolus is extremely slow in patients. A high, stabilized value of blood glucose characterizes a situation where the system cannot return to the glucose set point. This result may represent the situation when the liver is completely resistant to insulin, or, alternatively type 1 diabetes where no insulin is generated by the pancreas. Our model suggests that pancreatic control of blood glucose is possible as long as the sensitivity of the liver insulin receptor is at or above tI = 0.8 (1 − tI = 0.2 as shown in figure 3) when maximum blood glucose is within the range where the pancreas is sensitive.

For large, sustained glucose inputs, the model shows that blood glucose would increase without limit over time if there were no way of eliminating or using the excess. In reality, excess glucose is diverted into fat metabolism pathways.

For moderate glucose input and insulin resistance, pancreatic control results in good homeostasis, with low maximum blood glucose (the light grey region in figure 5). Insulin sensitivity above the normal range is characterized by underdamping of blood glucose oscillations.

Increased sensitivity to insulin may occur as a result of an autoimmune reaction, or pituitary or adrenal insufficiency. The onset of oscillations in a glucose tolerance test could provide an early indicator of insulin hypersensitivity. Currently, there is no clear consensus over the protocol for a glucose tolerance test. Most practitioners sample only once or twice following the administration of glucose, and there is no standard advice over the appropriate timing. In order to reveal the oscillations, it would be necessary to sample several times over 2 or 3 h to reveal the temporal picture. Although individual variation will be considerable, we suggest that sampling for 2–3 h at half-hourly intervals would reveal the temporal pattern of the response and indicate whether the response is oscillatory. Oscillations would alert the clinician early to the possibility of a change in insulin sensitivity.

A large transient excursion during the first oscillation, significantly larger than subsequent oscillations, also indicates insulin hypersensitivity. High insulin sensitivity can be a problem for individuals on insulin therapy for diabetes mellitus as a consequence of an autoimmune response. The consequences of insulin hypersensitivity are particularly marked in relation to fat metabolism [1214].

3. Discussion

The results we obtain from analysing the composite model suggest that analysis of the relationship between the glucose stimulus and the magnitude of the resultant blood glucose excursion may provide a more effective indicator of homeostatic capability than considering the glucose increase resulting from a single glucose input. In our model, effective homeostasis occurs over a region of glucose inputs where increasing glucose input does not result in increased blood glucose level. Further development of the model may reveal potential early clinical indicators of altered insulin resistance. This could be valuable because the model shows that unless insulin resistance increases by more than tI = 0.8 (1 − tI = 0.2), direct glucose control will compensate. Increased resistance indicates diabetes, whereas increased sensitivity may reflect an autoimmune response to insulin or pituitary or adrenal insufficiency. Therefore, both conditions benefit from early recognition and treatment. For example, the analysis of transients in glucose tolerance tests may have diagnostic power for the study of pathologies of glucose metabolism (figure 4b) which reveals that the magnitude of the transient response is determined by insulin sensitivity.

Finally, the model supports the proposal that ultradian oscillations result from system-level behaviour of feedback between liver and pancreas. Consider a situation in which glucose is fed into the blood at a constant positive rate M. This results in the indefinite growth of blood glucose concentration, other things being equal. However, when blood glucose concentration has reached a threshold level, release of insulin by the pancreas is triggered. This in turn results in an increase in phosphorylated glycogen synthase kinase (GSK), and hence in the increased rate of conversion of blood glucose into glycogen in hepatocytes (model F of Hetherington et al. [1], appendix A6). The increase in blood glucose concentration is then arrested, resulting in a stable equilibrium concentration, but may be reversed if GSK production is sufficiently sensitive to insulin (i.e. provided tI is sufficiently low—model D of Hetherington et al. [1], appendix A4 and eqn A4.1). If this reversal is fast enough, then blood glucose concentration can drop rapidly below the threshold at which insulin is released by the pancreas. In the absence of insulin, GSK then also drops rapidly, and the rate of conversion of glucose into glycogen is therefore reduced. Because glucose is still being ingested at a constant rate, blood glucose concentration may then again increase, and the cycle begins again.

This cycle of events leads to oscillations in blood glucose concentrations. These oscillations may be overdamped, or underdamped [15] or, result in a stable limit cycle as shown in figure 2 depending on M and other model parameters. For a given M, these outcomes depend essentially on the sensitivity of GSK to insulin. As the control parameter tI decreases, the insulin sensitivity increases and there is a sequence of bifurcations through resulting outcome blood glucose concentrations, i.e. from stable equilibrium but with high glucose level to overdamped oscillations leading to lower stable equilibrium, to the appearance of a stable limit cycle with lower mean concentration, leading finally to physiologically unsustainable underdamped oscillations. We can demonstrate this computationally using the model.

If we effectively remove the pancreas from the system—by setting the scaling factors between the pancreas and insulin models and pancreas and cAMP models to zero—then the system does not produce ultradian oscillations in blood glucose for values of the external glucose input (M) for which oscillations were observed in the complete system. Removing the pancreas from the system is effectively equivalent to setting the insulin sensitivity of the liver to zero. In the first case (no pancreas), no insulin is produced and, hence there is no response in the liver. In the second case (insulin sensitivity = 0), insulin may be produced by the pancreas but the liver is unable to respond to the hormone. The bifurcation analysis has shown that under these conditions oscillations do not occur (figure 2; bifurcation diagram).

If we incrementally increase the scaling between the pancreas and insulin models back to its original value, we then observe the reappearance of oscillatory solutions. The point at which this occurs depends on the glucose input (M) and the sensitivity of the liver to insulin (1 − tI); the higher the insulin sensitivity the lower the scaling value required to produce ultradian oscillations. These results do not depend on the value of the glucagon-scaling parameter because for positive glucose inputs the blood glucose level does not fall below the value required to trigger the production of glucagon in the pancreas model.

We can explore the role of the feedback of blood glucose on the pancreas by varying the parameters that govern the response of the pancreas to the changes in blood glucose. The production of hormones by the pancreas is modelled using threshold (or Hill) functions. The parameters (tIg and tLg) determine the level of blood glucose above or below the reference value at which the pancreas releases insulin or glucagon, respectively. Increasing these parameters represents the situation in which the pancreas is less sensitive to changes in blood glucose. As tIg is increased for a given M and tI, oscillations cease to exist as the pancreas does not respond to the initial rise in blood glucose with the release of insulin. The more sensitive the liver is to insulin the greater the range of tIg for which we observe oscillations (figure 6); the amount of insulin produced by the pancreas, even with a high glucose response threshold, is sufficient to trigger the regulation of GSK in the liver.

Figure 6.

The y-axis shows the value of tIg at which model oscillations cease. For lower values, the system oscillates.

4. Conclusions

This paper presents an exploration of the behaviour of the composite model of glucose homeostasis described in the companion paper [1]. We have explored the behaviour as a function of a range of glucose input values and also, because oscillatory behaviour is known to occur, we present a bifurcation analysis to provide quantitative information about the occurrence of these oscillations.

Analysis of the composite model and its results in this and the companion paper [1] reveals that direct glucose control, determined by the capacity of the liver to store glucose, plays an important role in the control of glucose homeostasis. Deficiencies in hormonal control generated by, for example, a reduction in insulin sensitivity, can be accommodated through the direct control mechanism over a surprisingly wide range of insulin sensitivities. The analysis suggests that the consequences of reductions in the sensitivity of the insulin receptor do not become significant until receptor sensitivity has fallen below 1 − tI = 0.2. The masking of clinically significant alterations in insulin sensitivity by the direct glucose-control mechanism has important consequences, because delays in recognizing the early signs of diabetes may have consequences for treatment.

The direct glucose-control mechanism depends on the capacity of the liver to store glucose as glycogen, which provides a significant buffer to alterations in blood glucose however they may be generated. Glycogen storage diseases form a significant body of contributory causes to defects in glucose homeostasis, particularly when the defect resides within the liver, rather than the muscles. Analysis of the part played by direct glucose control in glucose homeostasis should improve understanding of both causes and consequences of glycogen storage disease. The model presented here will provide a relatively simple way of testing the physiological consequences of metabolic defects that underlie glycogen storage disease because there should be directly demonstrable and measurable consequences for glucose homeostasis. Any reduction in the capacity of the liver to store glucose in the form of glycogen will, because of the significant role in glucose homeostasis, exacerbate any consequences of glycogen storage diseases.

The composite model presented here is one step in an ongoing process of developing more complete system models of liver function and enhancing modelling methodology. The current model already yields interesting results and provides a foundation for work, including other important elements of liver physiology. An analysis of glycogen storage disease would be valuable. Obvious additions to the core model include: an elaboration of glucose metabolism and storage to represent the physiological transition towards fat metabolism as cellular glucose levels increase and the representation of multiple hepatocytes with heterogeneous properties and linked to each other by gap junctions, to investigate the role of metabolic zonation in liver function. Our composite model could also be enhanced significantly by the inclusion of a more detailed model of the pancreas.


This work was funded as a Beacon Project of the United Kingdom 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 November 16, 2011.
  • Accepted January 19, 2012.


View Abstract