## Abstract

Diabetes and obesity present a mounting global challenge. Clinicians are increasingly turning to mechanism-based mathematical models for a quantitative definition of physiological defects such as insulin resistance, glucose intolerance and elevated obesity set points, and for predictions of the likely outcomes of therapeutic interventions. However, a very large range of such models is available, making a judicious choice difficult. To better inform this choice, here we present the most important models published to date in a uniform format, discussing similarities and differences in terms of the decisions faced by modellers. We review models for glucostasis, based on the glucose–insulin feedback control loop, and consider extensions to long-term energy balance, dislipidaemia and obesity.

## 1. Introduction

Energy homeostasis requires the coordination of several metabolic fluxes. In broad outline, these fluxes comprise (i) the use of nutrients as metabolic fuels, (ii) the disposition into the storage of any nutrient supply surplus to immediate demands and (iii) the mobilization of nutrients from these stores when required (e.g. Frayn 2003). Among the hormonal factors regulating these fluxes, insulin occupies a central position due to its wide-ranging effects on blood glucose levels, cell growth and appetite, as well as internal energy stores (Brook & Marshall 2001). The latter effect is predominantly anabolic: insulin stimulates glycogenesis, the accumulation of triglycerides in adipose tissue and amino acid uptake into muscle, while it inhibits gluconeogenesis, lipolysis and (in muscle) proteolysis (Brook & Marshall 2001). One could regard the insulin control loop as the fundamental ‘first-order’ endocrine control, with other factors (thyroxine, cortisol, glucagon, GH, ILG-1, adipokines, catecholamines) providing a ‘second-order’ adjustment to the precise physiological circumstances (Schwartz *et al.* 1987; Raju & Cryer 2005).

The coordination of metabolic fluxes at the systemic level requires an interplay between various organs and tissues, which each make different contributions to the overall energy balance. The breakdown of this coordination leads to abnormalities in energy stores and the blood levels of glucose, lipids, ketone bodies and electrolytes (e.g. Salway 2004). This web of mutually exacerbating abnormalities may stress the insulin control loop until its ability to compensate for the deviations is finally overcome. Consequently, hyperglycaemia is a common endpoint, whether or not the initial defect occurred in insulin production, secretion or sensitivity. The resulting syndrome of diabetes mellitus is associated with a variety of long-term complications, many of which are thought to arise from hyperglycaemia via protein glycation and oxidative stress (Brownlee 2001; Fajans *et al.* 2001).

The clinical management of hyperglycaemia (and attendant abnormalities in blood lipids and lipid stores) requires a thorough understanding of the interactions between the liver, pancreas, fat stores and blood nutrient levels. An explicit representation of this complex dynamical system can aid the interpretation of diagnostic tests and help formulate therapeutic interventions. To this end, a large range of various mathematical models have been proposed in the literature. The purpose of this paper is to bring together these mathematical models and present them, for the ease of comparison, in a uniform format. Thus, we present the various models in a standardized notation (summarized in table 1) and a uniform graphical format (*viz.* dynamical charts). We emphasize connections between the models in terms of different choices the modeller has to make, and discuss these choices with reference to the underlying physiology.

We stress the central ideas behind the models and the key insights that can be gained from them, as well as the relationships between models and their physiological grounding. The present paper thus complements previous reviews that explicitly list systems of equations for the models available in the literature (for such reviews see van Riel (2004); Boutayeb & Chetouani (2006); Makroglou *et al.* (2006)). We hope that the survey will prove useful to modellers as well as to endocrinologists and clinicians who wish to apply mathematical models in their research and require a general orientation.

## 2. Overview of energy homeostasis models

The main physiological and endocrine interactions underlying energy homeostasis in humans are depicted schematically in figure 1. This chart shows how the various models proposed in the literature fit together (generally in dynamics charts, state variables are shown as boxes, rates or fluxes as double arrows going into or out of these state variable boxes and flows of information (i.e. regulatory dependencies) as dotted arrows, while sources or sinks are shown as clouds). Plasma glucose and plasma insulin are at the heart of all these models, and therefore the survey will start with ‘glucocentric’ models that focus on these two factors. More elaborate models include the non-esterified fatty acid (NEFA) pool, as well as internal stores (glycogen for carbohydrate, adipose tissue for lipids, protein for amino acids); moreover, dynamically more refined models differentiate between interstitial pools and the general blood compartment as well as time delays for certain physiological response components. These elements will be discussed in more detail below.

Mathematically, the models typically assume the form of a system of (generally nonlinear) ordinary differential equations (ODEs); thus, each state variable (box) in figure 1 corresponds to a first-order ODE. Section 3.2 below discusses extensions to delay-differential equations (DDEs) and integro-differential equations (IDEs).

## 3. Models of glycaemic regulation

Insulin and glucose are at the heart of the dynamics chart (figure 1) and form the fundamental glycaemic homeostatic feedback loop. Only two state variables are needed for a model limited to this fundamental loop, and more advanced models can be viewed as variations of this simple two-state model.

### 3.1 Basic glucostasis

At the heart of glucostasis lies the basic glycaemic feedback loop represented in figure 2. The simplest mathematical model of this loop deals only with the plasma levels of insulin and glucose (Bolie 1961). The following ODEs are basic balance equations:(3.1)(3.2)where *I*^{[p]} is the plasma level of insulin and *G*^{[p]} denotes the plasma level of glucose. Insulin is secreted by the pancreatic *β*-cells at a rate *Φ*_{I}, treated as a function of *G*^{[p]} alone. Insulin is cleared from the plasma with a physiological half-life equal to , where *λ*_{I} is the insulin clearance rate constant. More generally, it is reasonable to suppose that *Φ*_{I} also depends on *I*^{[p]}; Kulkarni *et al.* (1999) showed that in *β*-cell insulin receptor knockout mice, first-phase insulin secretion in response to glucose is ablated. It is known that *β*-cells express an insulin receptor, although the regulation of *β*-cell proliferation may be the primary function (Okada *et al.* 2007). A case could also be made for an explicit dependence on *t* of the term *Φ*_{I}, to account for the parasympathetic influence on *β*-cells.

The gastrointestinal tract is the main exogenous source of glucose. Endogenous sources are glycogen mobilization and hepatic gluconeogenesis from muscle protein or muscle glycogen via pyruvate or lactate (Frayn 2003). The basic model lumps all these processes together into a single term *Φ*_{G}, whose explicit dependence on time *t* accommodates the dependence on the exogenous input regime. Glucose transport into the tissues is represented by *Ψ*_{G}. This is also a lumped term, combining (i) *insulin-independent* transport that depends only on *G*^{[p]}, (ii) *insulin-dependent* transport that depends on *I*^{[p]} as well as *G*^{[p]} (Frayn 2003) and (iii) renal excretion of glucose; this flux is important in hyperglycaemia (Watkins 2003).

The mode of action of this negative feedback mechanism is intuitively clear (cf. figure 2): increased *G*^{[p]} promotes insulin secretion, and insulin in turn promotes the clearance of glucose from the plasma, stimulating both the formation of glycogen storage and the usage of glucose as a metabolic fuel in the insulin-dependent tissues (Brook & Marshall 2001). In the simple two-state model, these are sources and sinks outside the scope of the model, whereas more comprehensive models treat the internal stores as explicit state variables. The dynamics of the basic model is illustrated in figure 3.

#### 3.1.1 Analysis of the basic model: quantification of insulin sensitivity

Clinical management of diabetes is aimed at a restoration of *G*^{[p]} to normal values, or at least prevent excursions beyond the range of concentrations found in healthy subjects (Watkins 2003). This is reasonable in view of the fact that diabetic complications (such as vascular disease, neuropathy, retinopathy, nephropathy) are thought to result from hyperglycaemia via a small number of common pathways (Watkins 2003).

Suppose that the exogenous glucose input is held constant, so that the glucose gain term loses its explicit time dependence: . The plasma glucose equilibrium value is defined implicitly by(3.3)(from equations (3.1) and (3.2)). This equation may have zero, one or more solutions and the precise functional form of the flux terms determines whether any of these is an attractor of the system's dynamics. On the basis of the physio-endocrinological interpretation of the flux terms, it is reasonable to assume that is monotone increasing in , while is monotone decreasing in both its arguments and is monotone increasing in both its arguments. Under these assumptions, it can be shown that if equation (3.3) has a solution , it is unique and locally a stable attractor. The extent of the hyperglycaemic excursion (e.g. the ‘area under curve’ or AUC, which is the time integral of following a glucose challenge, where denotes the normoglycaemic value) depends on the slowest mode of the system (the fast mode is essentially *λ*_{I}): a faster slow mode implies a smaller excursion. The slow mode corresponds to the *insulin sensitivity* of the system, given by the following formula (to be evaluated at the equilibrium point):A diminished ability of the pancreatic cells to respond to a glucose challenge may be due to a reduced *β*-cell mass, or reduced responsiveness of individual *β*-cells or both (Watkins 2003). This diminished pancreatic responsiveness corresponds to a reduction of the term . A major factor in the aetiology of diabetes is *insulin resistance* (Watkins 2003), a reduced effectiveness of insulin on some of insulin's target tissues due to a reduced sensitivity to insulin. In the model, insulin resistance is reflected by a reduction of the second term, which will normally be dominated by . Further details on the use of these ideas in the analysis of glucose tolerance tests can be found in the papers by Toffolo *et al.* (1980), Steil *et al.* (1993), Saad *et al.* (1994) and Prigeon *et al.* (1996).

#### 3.1.2 Functional specification of the flux terms

The system, equations (3.1) and (3.2), leaves the flux terms unspecified. Simple obvious choices include linear models of the form(3.4)where *Φ* is some flux dependent on a factor *X* and *Φ*_{0} and *α*_{X} are (possibly negative) parameters, and Hill-type expressions:(3.5)where *Φ*_{X}, *K*_{X} and *m*_{X} are positive parameters; *K*_{X} is the value of *X* such that ; and *m*_{X} governs how steeply the function rises about this midpoint. Fluxes that depend on multiple factors (*X*,*Y*) are often taken to be bilinear terms ∝*XY*. Inhibitory interactions are often modelled in this way, with the risk that the flux spuriously changes sign somewhere during the simulated process. A more prudent choice for an inhibitory interaction is multiplied Hill terms, e.g.(3.6)where *Y* represents the concentration of a factor that inhibits the flux, characterized by location and steepness parameters *K*_{YX} and *m*_{YX}, respectively.

The justification for such specifications is limited. It is hoped that they are qualitatively correct and it is tacitly hoped or assumed that the details will not affect the outcome appreciably. This may be reasonable in view of the uncertainties that accrue on a more mechanistic approach, which split the lumped terms into numerous terms representing processes in various different tissues. Even if each term could be modelled accurately, reliable parameter estimation may be difficult. This provides some rationale for the cruder, simpler approach. Another difficulty is that the molecular transport machinery is embedded in cellular regulation. To indicate the problems encountered, we consider the example of insulin-dependent glucose transport.

##### Insulin-dependent glucose transport

A simple model for insulin-dependent glucose uptake is the bilinear model . If we take into account all relevant sub-cellular processes, such as intracellular insulin signalling, glycolytic pools of hexoses and trioses and their regulation, we end up with a dynamical system with greater than approximately 100 state variables. For example, Sedaghat *et al.* (1992) formulated a model with 20 state variables for the insulin/GLUT4 part of the system alone (GLUT means glucose transporter, GLUT4 being the insulin-responsive variant; see Silverman 1991). The purposes at hand dictate which of these approaches is most suitable. For instance, in the analysis of a glucose tolerance test, it suffices to know that GLUT4 translocation to the cell membrane in the presence of insulin happens quasi-exponentially with a relaxation time of approximately 5 min; this reduces the 20 state variables in the Sedaghat model to just one. On the other hand, when one is interested in how a drug affects this relaxation time, the 20-dimensional model is appropriate.

The influence of intracellular metabolites can be appreciated by considering the following simple, but mechanistically reasonable, expression for glucose transport into a cell (see Keener & Sneyd (1998) for a derivation):(3.7)Here *ϑ*, *K* and *K*_{d} are positive parameters; *n*_{GLUT} is the number of GLUT molecules translocated to the cellular membrane; and *G*^{[i]} is the intracellular concentration of glucose. Let denote the rate at which glucose is converted into glucose-6-phosphate (the first metabolic conversion that immediately follows uptake). Here the vector * w* collects the intracellular factors that affect this rate. If

*is fixed,*

**w***G*

^{[i]}is found by solving ; this flux can then be computed from

*n*

_{GLUT}and

*G*

^{[p]}. However, the situation is more complicated since any changes in

*G*

^{[i]}generally lead to compensatory changes downstream in glycolysis and ATP generation. This means that

*will change as well, and this must be taken into account to determine the shape of the response curve, which describes as a function of*

**w***G*

^{[p]}. A detailed analysis shows that the asymptotic maximum of this curve is determined by metabolic demand rather than

*n*

_{GLUT}, and that the midpoint of this curve is an increasing function of the ratio between metabolic demand and

*n*

_{GLUT}. If we approximate the actual response curve by a Hill-type function, we have that the ‘midpoint’ parameter

*K*

_{X}is (i) an increasing function of metabolic demand and (ii) inversely proportional to

*n*

_{GLUT}. In tissues that express GLUT1, which is not insulin responsive,

*n*

_{GLUT1}does not depend on

*I*

^{[p]}and the midpoint parameter can be expected to rise linearly with the demand flux. By contrast, in tissues that express GLUT4,

*n*

_{GLUT4}depends on

*I*

^{[p]}and the dependence of the midpoint parameter on demand becomes more complicated, being insensitive at low demand and very sensitive at high demand. These different responses necessitate a decomposition of

*Ψ*

_{G}into several terms, corresponding to the fluxes borne by tissues expressing different glucose transporters.

An important consequence of the role of intracellular regulation in this model is that it explains why the effectiveness of insulin is reduced under hyperglycaemia or hyperlipidaemia. According to the model, insulin levels do not affect the maximum flux, but only act to reduce the midpoint parameter. This allows the effective regulation of insulin-dependent glucose uptake as long as plasma glucose levels remain sufficiently low. The critical level where effective control starts to fail is lower when metabolic demand for glucose is lower, for instance when *β*-oxidation of fatty acids provides a substantial flux of AcCoA.

### 3.2 Glycaemic loop models with time delays in cellular responses

Cellular responses are not instantaneous and thus introduce delays into the system. Pancreatic *β*-cells and liver cells introduce such delays in their responses to glucose and insulin, respectively (figure 4). These delays can introduce dynamical instability, i.e. oscillations in the dynamics.

#### 3.2.1 ‘Hard’ delay models

Pancreatic *β*-cells introduce a time delay as the insulin secretion rate takes some time to adapt to the plasma glucose levels. From a modelling perspective, one option is to represent this delay explicitly(3.8)where *τ*_{I,d} is the response delay incurred by the pancreatic *β*-cells as they respond to changes in the plasma glucose levels. A hard delay *τ*_{I,d} in the pancreatic *β*-cell response was investigated by Bennette & Gourley (2004) and Li & Kuang (2007).

A generalization of the delay replacement (equation (3.8)) is as follows:(3.9)where is a new auxiliary variable, defined by(3.10)where *ω* is a convolution kernel function. The hard delay, equation (3.8), is retrieved by taking the kernel *ω*(*s*) to be the *Dirac delta function* .

#### 3.2.2 ‘Soft’ delay models

Many delay models proposed in the literature use the *gamma kernel*(3.11)The advantage of this kernel is that it can be implemented by extending the dynamics with a finite number of state variables. Thus, *n* auxiliary state variables are defined, with the following dynamics:(3.12)(3.13)and *x*_{n}(*t*) is used instead of . For *n*→∞, the gamma kernel approaches the hard Dirac kernel. The effects of delays can often be effectively mimicked with a relatively low *n*. At least for small *n*, the mathematical analysis of such a soft delay approximation can be easier to deal with than the corresponding hard delay formulation, for instance in local stability analysis.

There is also a delayed response in hepatic glucose production from glycogenolysis. Tolić *et al.* (2000) considered a soft hepatic time delay (*n*=3) *τ*_{G,d} for as input to *Φ*_{G}(3.14)whereas Engelborghs *et al.* (2001) investigated the hard delay case. More complex models incorporate both the pancreatic delay and the hepatic delay. These have been studied by Li & Kuang (2001) and Li *et al.* (2006). Explicit time delays constitute a fairly phenomenological approach to cellular response lags. Li *et al.* (2006) pointed out that a physiologically realistic model would have to describe the kinetics of intracellular pools (e.g. the hexose phosphate pool and the triose phosphate pool) as well as the activation/inactivation kinetics of key enzymes such as glucokinase, glucose 6-phosphatase and phosphofructokinase. Inasmuch as a physiologically realistic model would add a finite (and preferably modest) number of additional state variables, it would introduce a soft rather than a hard delay. One could thus argue that an *n*=3 soft delay is a more reasonable ersatz for a physiologically realistic model than an explicit delay.

#### Other delay kernels

De Gaetano & Arino (2000) considered a two-ODE minimal model extended to a DDE model with a generalized delay on *β*-cell insulin secretion. For the glucose kinetics they assume the Roy specification (equation (3.22)) with , whereas for the insulin secretion rate they adopt the specification (with positive parameter *η*)(3.15)Finally, they assume the following *uniform kernel*:(3.16)Based on these assumptions, De Gaetano & Arino (2000) showed that the dynamics admits no more than one equilibrium with both *I*^{[p]} and *G*^{[p]} positive. Mukhopadhyay *et al.* (2004) showed that this result obtains for this model in the more general case, where *ω* is a probability density function whose mean exists; the Dirac kernel, the gamma kernel and the uniform kernel all belong to this class. Li *et al.* (2006) studied the same model, but replacing specification equations (3.22) and (3.15) by more general functions.

### 3.3 Extended models for glucostasis

The simple two-state variable model is based on the implicit assumption that other dynamic degrees of freedom relax more quickly than those included in the simple model. However, this assumption is not warranted in general. Modellers have thus turned to models with additional state variables, as illustrated in figure 5.

#### 3.3.1 Endocrine dynamics

The model (equations (3.2) and (3.1)) lumps exogenous and endogenous glucose inputs into a single term *Φ*_{G}. The endogenous processes (glycogenolysis and gluconeogenesis) are controlled by insulin and various other endocrine factors, the most important of which is glucagon. Making this regulatory loop implicit, the positive flux term is replaced by an exogenous and an endogenous term,(3.17)where *A*^{[p]} is the glucagon concentration in the blood plasma (the term still lumps together glycogenolysis and gluconeogenesis). The exogenous input is a *forcing function*: it depends on *t*, reflecting the fact that the time course of glucose supply is controlled by external factors. The term can be calculated as the rate at which glucose is absorbed from the diet or supplied by infusion, divided by the volume of the *glucose distribution space* (e.g. Roy & Parker 2006). To close the system, an ODE for *A*^{[p]} must be supplied, e.g.(3.18)where is the physiological half-life of glucagon and is the production function dependent on *G*^{[p]} and *I*^{[p]}. The dependence on *G*^{[p]} is included because glucagon is secreted in response to low glucose levels (Brook & Marshall 2001) whereas the dependence on *I*^{[p]} accommodates the inhibitory effect of insulin on glucagon secretion (Greenstein & Wood 2006). A quasi-steady-state approximation for *A* is , which provides a justification for the dependence of *Φ*_{G} on *G*^{[p]} in the two-state variable model described previously. If the desired time resolution is no finer than , the original formulation can be retained.

#### 3.3.2 Interstitial space: tissue compartments

The physiological response to insulin may appear to lag behind the response expected on the basis of the insulin plasma level *I*^{[p]}. Part of this lag can be explained from the distribution kinetics of insulin. To account for this, one may add one or more tissue compartments to the model, also known as *remote compartments* (Roy & Parker 2006). For each compartment, an additional ODE must be supplied to describe the *interstitial fluid concentration* of insulin in the *ℓ*th tissue compartment(3.19)This equation models a passive exchange process with transport parameter (e.g. Sturis *et al.* 1991; Tolić *et al.* 2000). To balance these exchange fluxes, equation (3.1) must be modified, as follows:(3.20)where the sum is taken over all additional tissue compartments accounted for in the model. Glucose kinetics (equation (3.2)) must also be adapted to reflect the fact that the endogenous fluxes respond to the *interstitial* insulin concentration rather than the plasma level. Thus, equation (3.2) becomes(3.21)Here, *ℓ*=0 denotes the hepatic tissue compartment, where endogenous glucose generation from glycogenolysis and gluconeogenesis is realized (muscle glycogenesis yields three-carbon intermediates that are converted into glucose in the liver). The index *ℓ* is included in the usage terms to reflect the possible variation in insulin sensitivity among tissue compartments. A model with three interstitial compartments (following the central plasma compartment) for glucose and two for insulin is analysed in detail in the excellent paper by Silber *et al.* (2007). Their model illustrates another use for interstitial insulin compartments: *viz.* to simulate the slow-release effect from a bolus injection (bolus injections are represented mathematically as Dirac inputs; such an input at time *t*=0 can equivalently be implemented by a non-zero initial condition).

When the transport coefficient is assumed to be equal for all compartments ( ), the interstitial insulin concentrations will likewise converge to a common value. Discounting the initial transient, a single state variable then suffices. The resulting model, figure 5*b*, is a three-ODE model known as the *standard minimal model*, first proposed by Bergman *et al.* (1981). Roy & Parker (2006) proposed simple functional specifications for the glucose fluxes:(3.22)where *Φ*_{0} is a baseline glucose release rate from the liver; the coefficient *α* expresses the inhibitory effect of insulin on gluconeogenesis (Gibney *et al.* 2003; Salway 2004); *β* corresponds to insulin-independent glucose usage; and *γ* accounts for insulin-dependent glucose usage. In a diabetic subject with , the plasma glucose concentration then relaxes towards the value , which is the subject's *baseline glycaemia*.

Of course, additional interstitial state variables can be introduced for other soluble factors (e.g. the model described in §4.1 accounts for interstitial fatty acids). Simpler models omit such compartments. The justification is analogous to the above argument for glucagon dynamics: thus, for any substance *X*, if is sufficiently large, the tissue compartment lag will not be appreciable and (where *X* is any substance) can be replaced by .

A generalized interstitial/plasma exchange model favoured by some authors (e.g. Roy & Parker 2006) is as follows:(3.23)where if *x*>0 and if . This model has transport coefficients and , which may have different values, for transport into and out of the interstitial compartments. The baseline parameters are positive constants that ensure that transport does not take place when the concentration in the donor compartment is below the baseline.

## 4. Integrated models of carbohydrate, fat and protein metabolism

A key aspect of energy homeostasis is the coordinated use of various metabolic fuels. The glycaemic feedback loop is central to this coordination, since insulin also serves as a regulator of the use of lipids and protein as metabolic fuels (figure 1; for the sake of conciseness, various endocrine factors that modulate insulin's action are omitted). The link with the usage of fatty acids as a metabolic fuel is particularly important, and has been the focus of attention in integrated models, while the link with protein metabolism has received less attention.

### 4.1 Integrated lipid/glucose models

The *glucose–fatty acid cycle* is central to the interactions between carbohydrate and lipid metabolism (Randle *et al.* 1963). Through its effects on the metabolic control of the glycolytic pathway, the oxidation of fatty acids reduces the uptake and oxidation of glucose, with the consequence that glucose usage is reduced compared with that expected at given concentrations of insulin and glucose in the plasma; this manifests itself as *insulin resistance* (as discussed in §3.1.1 and 3.1.2). The glucose load following a meal will activate insulin secretion, which stops the release of fatty acids from the adipose stores and promotes the disposition of plasma NEFAs towards these stores via low-density lipoproteins, while glucose and insulin drop to baseline levels between meals, promoting lipolysis and the usage of NEFAs (Frayn *et al.* 1993; Frayn 2003).

#### 4.1.1 Fatty acid usage

The key characteristics of the glucose cycle have been incorporated in the model proposed by Roy & Parker (2006). This model is based on the standard minimal model (see §3.3) modified to have two interstitial compartments for insulin and an additional interstitial compartment for NEFAs, which communicates with a plasma NEFA compartment, a new state variable *F*^{[p]} (see the dynamics chart, figure 6). One of the interstitial insulin state variables drives glucose disposition (marked ‘ins.-dep’ in figure 6) whereas the other drives disposition of NEFAs in adipose storage (a sink/source in this model). There are no endogenous sources of insulin, so insulin gains are due to external inputs (there is an additional production term that causes insulin degradation to change sign at a ‘baseline value’). The model includes the standard dietary gain terms for plasma glucose and plasma NEFA. Fluxes that depend on a hormone concentration and a substrate are modelled as bilinear terms (e.g. ).

Noteworthy elements in the model are two modulatory interactions. First, Roy & Parker (2006) have included the inhibitory effect of plasma NEFA on glucose disposition, as a bilinear gain term (again they include an additional gain term that causes a sign change at baseline values). The model thus accommodates a form of insulin resistance. In addition, plasma glucose (*G*^{[p]}) and plasma NEFA (*F*^{[p]}) modulate lipolysis (NEFA mobilization from adipose stores); this term is of the form (with *G*_{0} a positive parameter) and represents modulatory interactions that are experimentally observed (Knight & Iliffe 1973). Roy & Parker (2006) obtained good model fits to experimental data, which include (i) decay of plasma NEFA under euglycaemic hyperinsulinaemic clamp, (ii) glucose disposition under intralipid infusion and (iii) decay of *I*^{[p]} and *G*^{[p]} in response to an *intravenous glucose tolerance test*.

#### 4.1.2 The glucose cycle with internal store dynamics

A model similar to the Roy & Parker (2006) model was proposed by Maas & Smith (2006); a dynamics chart is shown in figure 7. The model differs from the previous model in that it explicitly represents internal reserves (glycogen and triacylglyceride (TAG) stores) as state variables. Both insulin and glucagon are represented as state variables. There are no interstitial compartments. The model contains an additional state variable that evolves on a slow time scale. Slow dynamics will be discussed below in §5.2.1.

Noteworthy features of this model are the inclusion of the inhibitory effect of insulin on glucagon secretion (modelled as a bilinear term ), and the division of energy usage into a component that is obligatory glucose fuelled and a component that can use either glucose or NEFAs. In the latter component, the ratio of glucose usage over NEFA usage is assumed to be proportional to the corresponding ratio of plasma levels, . The model also includes glycogen autoregulation, which is a physiological phenomenon where the rate of glycogen accumulation diminishes as the glycogen store size approaches a maximum capacity (Laurent *et al.* 2000). Maas & Smith (2006) represented this by making the corresponding term proportional to the difference between the current glycogen store size and the maximum value.

### 4.2 The need to integrate glucocentric models with protein metabolism

Muscle protein represents a considerable energy store. Under conditions of prolonged starvation or disuse atrophy, this store is mobilized and the amino acids derived from muscle proteolysis either serve as metabolic fuel or are converted into glucose or fat (Karl *et al.* 1976; Salway 2004). Protein breakdown is connected to the glycaemic feedback loop since insulin is involved in the regulation of proteolysis (figure 1). Proteolysis is inhibited by insulin and promoted by cortisol and triiodothyronine (Karl *et al.* 1976; Fereday *et al.* 1998). Moreover, alanine is a substrate for gluconeogenesis. In view of these important links, it is perhaps surprising that protein metabolism is virtually ignored by the foregoing glucocentric models. However, this can be justified. To explain this, we first need to discuss the interface between carbohydrate and protein metabolism.

Amino acids can serve as metabolic fuel: muscle breaks down branched-chain amino acids, the liver covers half its energy requirements from amino acids and glutamine is an important metabolic fuel for mucosal intestinal cells (Frayn 2003). However, the use of protein for oxidation is a limited portion of total energy requirements in a well-fed subject (Frayn 2003). Second, while alanine certainly serves to carry *nitrogen* to the liver (which is cleared as urea), the *carbon* in alanine may derive primarily from pyruvate generated by glycolysis of muscle glycogen and glucose taken up by the muscle (Garber *et al.* 1976*a*,*b*; Karl *et al.* 1976). Alanine is thus just a tricarbon shuttle in a pathway from muscle glycogen reserves to plasma glucose via gluconeogenesis in the liver. The justification of models of the type considered in §3 is that these effectively subsume the alanine pathway, together with other gluconeogenic precursors, in the overall flux of *glucose* derived from a ‘lumped’ source, which accounts for hepatic as well as extra-hepatic glycogen reserves.

Nevertheless, the representation of the glucogenic precursor pathways would be appropriate in more detailed models that explicitly represent muscle and liver glycogen stores as separate state variables. Besides the alanine cycle, such models would also account for the *Cori cycle* in which muscle releases the gluconeogenic precursors lactate and pyruvate (Salway 2004). It has become clearer in recent years that these stores follow different dynamics and that extra-hepatic glycogen serves an important role as a postprandial buffer of the carbohydrate load ingested with the meal (Taylor *et al.* 1993, 1996; Meyer *et al.* 2002). In fact, Woerle *et al.* (2002) found that almost a third of glycogen formation following a meal may be via the indirect pathway, via a gluconeogenic precursor such as alanine, pyruvate or lactate, some of which are temporarily stored as extra-hepatic glycogen. Models that treat liver and muscle glycogen stores as important separate components are appropriate for the analysis of these recent studies. Moreover, such models seem certain to be relevant to disorders of energy homeostasis, since the dynamic buffer function of muscle glycogen is impaired in type 2 diabetes (Carey *et al.* 2002). Moreover, extra-hepatic glycogen stores play a role in the response to a positive energy balance (see §5.2.2).

### 4.3 Glucoplastic versus ketoplastic carbon reserves

Glycogen reserves and muscle protein together constitute reserves of moderately to highly oxidized carbon, in contrast to lipid reserves, which store carbon almost exclusively in a more reduced form. This distinction is important not only from the point of view of energetic density per unit weight (Frayn 2003), but also because mammals cannot synthesize glucose (or any other metabolite derived from the Krebs cycle intermediates succinyl-coA through to oxaloacetate) from fatty acids; while carbon atoms will readily pass through *β*-oxidation and the Krebs cycle, there can be no net gain since two carbon atoms are lost as CO_{2} for every two atoms derived from fatty acids (Salway 2004). The distinction between the two kinds of carbon is physiologically important since the *ketoplastic* carbon atoms of lipid reserves cannot be turned back into *glucoplastic* carbon atoms, whereas the latter can still be converted into ketoplastic carbon (the glycerol carbons of lipid reserves remain glucoplastic, of course; moreover, since one anabolic pathway branches off between the CO_{2}-evolving steps, there are amino acids such as glutamine that must be reckoned as one-fifth ketoplastic, and thus protein is in fact a mixed store of glucoplastic and ketoplastic carbon; similarly, citrate and isocitrate are one-third ketoplastic). Consequently, ketoplastic reserves are less versatile as anabolic precursors (i.e. building blocks for growth): most of the carbon in the lipid reserves will have to be used as fuel—lipid reserves are ‘deferred catabolism’. On the other hand, carbon stored as glycogen can to some extent be ‘deferred growth’, built up in times of shortage of high-quality nutrition (minerals, nitrogen). Deferred growth is not only relevant to juveniles, but also to adults: while adults no longer grow, the females do supply building blocks to growing offspring. In sum, the rate *at the whole-body level* of conversion of glucoplastic to ketoplastic carbon represents a key life-history trade-off; if this rate does not, over prolonged periods of time, match the rate at which carbon is oxidized to CO_{2}, obesity ensues.

## 5. Long-term dynamics

The models considered thus far were concerned with the physiological time scale at which the system absorbs acute glucose challenges. This is the minutes-to-hours time scale that is relevant for the response to short-term starvation between meals and for nutrient disposition following a meal. However, the processes underlying the development of obesity and diabetes (and their interplay in insulin resistance syndrome, or *metabolic syndrome*) take place on a longer (days-to-years) time scale. Here we discuss various modelling approaches to this longer time scale.

### 5.1 Diabetic aetiology

Diabetes can develop as a primary dysfunction of glucostasis or secondarily when the demands placed on the glycaemic control system overstress its capabilities. We discuss two modelling approaches to the long-term compensatory mechanisms.

#### 5.1.1 Regulation of the pancreatic *β*-cell mass

The basic glycaemic feedback loop has what engineers call *finite gain*: the loop does not effectuate the perfect regulation of plasma glucose levels. From equation (3.17) we havewhere is the fixed external glucose load. Thus, the function is different for different glucose loads, and so is the glucose equilibrium point, via equation (3.3). However, a slow-acting second feedback mechanism operates, which ensures (at least under physiological conditions) that the long-term average plasma glucose is regulated to a value that does *not* depend on the average load. The mode of action is essentially what a control engineer would call *integrating control* (cf. Milhorn 1966). This slow mechanism is based on the dynamics of the pancreatic *β*-cell mass(5.1)where is the *β*-cell mass; is the *β-cell proliferation function* that relates the rate of *β*-cell proliferation to the plasma glucose level (glucose is known to potentiate the effect of IGF-1 that stimulates *β*-cell proliferation; Hügl *et al.* 1998); and is the *β*-cell death rate. The specific pancreatic *β*-cell growth rate not only depends on plasma glucose, but has recently been shown to depend critically on insulin as well (Okada *et al.* 2007). Thus, a more realistic model would make dependent also on *I*^{[p]}.

To connect *β*-cell mass dynamics to the glycaemic loop model, the following replacement is made:(5.2)where *Φ*_{I} is the insulin secretion rate per pancreatic *β*-cell. The resulting three-state variable model is depicted in figure 8. Like *Φ*_{I}, the cell-specific function *Φ*_{I} is assumed to be monotone increasing. With the assumptions as stated in §3.1, becomes a monotone decreasing function of . This function can be substituted into equation (5.1) to obtain a scalar ODE for the slow time scale.

Let us assume, in the first instance, that *μ*(.) is monotone increasing (i.e. glucose promotes *β*-cell proliferation). It is then easy to see that, on the slow time scale, is a monotone decreasing function of with a single zero corresponding to a stable global attractor. Thus, will tend to the long-term equilibrium value regardless of the glucose load (of course, variations in the glucose load will determine the extent of variations of *G*^{[p]} around this value). This effect is typical of integrating control.

##### Catastrophic failure of the long-term loop.

Topp *et al.* (2000) speculated that glucose tends to *inhibit* *β*-cell proliferation at high concentrations. Thus, becomes a hump-shaped function, which first increases and then decreases (without becoming negative). As a result, now has two equilibrium points, the lower of which is unstable. Consequently, sustained hyperglycaemia can induce net *β*-cell loss, which promotes further hyperglycaemia. The result is a catastrophic feed-forward loop. Thus, as normal hyperglycaemia episodes become more frequent or longer, *β*-cell mass dynamics may be pushed past this catastrophic point: Topp *et al.* (2000) analysed this transition in detail and discussed the mechanism as a possible route to diabetes. This analysis raises the possibility of *β*-cell mass involution leading to *non-autoimmune*, fulminant type 1 diabetes mellitus; interestingly, precisely such a syndrome has been recently described by Imagawa *et al.* (2000).

#### 5.1.2 Evaluating therapeutic regimes

De Winter *et al.* (2006) took a different approach to the long-term dynamics of diabetic aetiology. To describe the deterioration of *β*-cell function and insulin sensitivity, they introduced two functions, denoted by *B*(*t*) and *S*(*t*), corresponding to the remaining population of fully functional *β*-cells and insulin sensitivity, respectively. The long-term dynamics of plasma glucose is described by the following IDE:(5.3)where angled brackets indicate long-term averages; *μ*_{G} and *μ*_{I} are rate constants; denotes the optimal plasma glucose concentration; and and are two forcing functions describing the effect of therapeutic intervention. To assess the success of the therapeutic regime an output function is required; for this De Winter *et al.* (2006) used(5.4)where *H* represents damage (glycosylated haemoglobin) and *μ*_{H} is the haemoglobin turnover constant. The usefulness of this approach hinges on the choice of two empirical functions (*B*(*t*) and *S*(*t*)) that represent the progressive loss of functionality; De Winter *et al.* (2006) took these to be logistic sigmoid curves. A conceptual difficulty with this approach is that the therapeutic intervention may have an effect on the disease progression, i.e. on the functions *B*(*t*) and *S*(*t*) (since therapy affects and , which may change the stress on *β*-cells and insulin sensitivity). In De Winter *et al.*'s (2006) formulation, these functions are autonomous. This dynamic feedback can be taken into account by augmenting the system with ODEs for *B* and *S*.

### 5.2 Adiposity dynamics in positive energy balance

When energy intake exceeds energy expenditure for a prolonged period (even by only a few per cent), the accumulation of stored energy will occur. Excess carbohydrate and amino acids are largely converted into lipids (although net lipogenesis is likely to play a major role only if the ingested diet is unusually rich in carbohydrates; Aarsland *et al.* 1997; Frayn 2003). Thus, a food intake that supports a positive energy balance over the long term results in the laying down of fat reserves. When the ability of adipocytes to absorb this supply is taxed to capacity, elevated plasma fatty acids levels may result. Since NEFA usage is regulated in a simple donor-control fashion by the availability of plasma NEFA (Izzekutz *et al.* 1967), glucose usage will be reduced and sustained elevated levels of plasma glucose and insulin will result. Moreover, sustained insulinaemia that extends beyond the normal postprandial peak could contribute to insulin receptor downregulation (Okabayashi *et al.* 1989), further reducing insulin sensitivity. The key step in this route to diabetes via dislipidaemia is the stress on the assimilatory capacity of the adipocytes.

#### 5.2.1 Adipocyte proliferation

Maas & Smith (2006) modelled the rate at which plasma NEFA is absorbed into adipose reserves as a bilinear term , where is the number of adipocytes in the body. Thus, the model simulates the elevated lipidaemia and hyperinsulinaemia when the adipocyte count is low in comparison with the energy excess in the dietary intake. Adipocyte resistance to the action of insulin is represented here by the quantity ; this parameter may be under hormonal control, since adipocytes secrete *resistin* that reduces the efficacy of glucose transport into adipocytes (Faust *et al.* 1978; Steppan *et al.* 2001) and may likewise affect insulin-stimulated NEFA uptake.

Adipocyte proliferation, i.e. an increase of , can compensate for the elevated *I*^{[p]} and *F*^{[p]}. Maas & Smith (2006) modelled this by taking to be a state variable with slow dynamics (‘adipocyte count’ in figure 7). Maas & Smith (2006) modelled the rate of proliferation phenomenologically, letting depend on the mean adipocyte lipid content via a very steep Hill function. If *Q*_{TAG} denotes the total TAG content of adipocytes, the mean adipocyte lipid content is equal to (this is represented by the quotient symbol in figure 7). This model is physiologically plausible: individual adipocytes greater than a certain size (and hence lipid content) are not observed, suggesting that adipocytes that have reached their maximum girth emit a signal that stimulates the proliferation and/or differentiation of preadipocytes. This is the *critical fat cell size hypothesis* (Hausman *et al.* 2001). In addition to a peripheral paracrine mechanism, central modulation via a neuroendocrine axis may also play a role (Steppan *et al.* 2001). The effect of the steep Hill function is to keep the ratio close to the upper size limit of adipocytes as long as positive energy balance persists.

In mathematical terms, it can be shown that the Maas & Smith (2006) model implies that TAG reserves over the long term exhibit a sigmoid (not a first-order) relaxation towards an asymptotic value that is proportional to , where the quantities in angle brackets 〈.〉 indicate long-term averages and is assumed to be kept constant (or very nearly so) by the steep ‘adipostat’ function. When glucagon levels are low, as they are likely to be under conditions of positive energy balance and insulinaemia, TAG reserves grow without bound as approximately as long as is kept at a fixed value. However, the increase in *Q*_{TAG} over the long term is already mathematically determined. The slow dynamics of *Q*_{TAG} is fixed by the average positive energy balance in the long term, which is imposed as a forcing function in this model. This mathematical determination is very interesting in physiological terms, since it means that is effectively governed by the imposed energy imbalance. To the biologist, this may seem counter-intuitive since the effect seems to turn causality on its head. This paradox is resolved by the interplay between the fast and slow time scales.

The positive correlation between adiposity and is an interesting model property; since must in the long term reflect the imposed dietary excess, the quantity must be greater than normal during adiposity gains. In other words, hyperinsulinaemia or hyperlipidaemia is a *necessary* concomitant of adiposity gains in this model. Moreover, they combine multiplicatively; is essentially an overlap integral. In normal subjects, *F*^{[p]} quickly falls as soon as *I*^{[p]} rises (Frayn *et al.* 1993). Thus, since the overlap is normally quite small, it may be increased appreciably by shifts in the timing of NEFA clearance even when mean levels do not change by much.

The Maas & Smith (2006) model explains why weight gains occur so readily after a diet. If adiposity before the diet is and after the diet , then the rate of fatty acid disposition into storage will be times elevated, *ceteris paribus*, when compared with normal subjects who have always been at adiposity , because the number of adipocytes in the post-diet subjects is elevated by that factor, due simply to the pre-diet adaptation of to via the steep Hill function, combined with the model assumption that adipocytes do not die. This elevated lipogenesis predisposes to renewed adiposity gains, especially if the subjects eat more to compensate for this enhanced flux into storage. These predictions are more likely to be relevant in middle-aged and older subjects, as a result of a second term that Maas & Smith (2006) included in the dynamics of *Q*, besides the steep adipostat Hill function. This second term is proportional to the body growth rate (which they impose as a forcing function). As a result, the effects noted here are masked in young adult life since the adolescent increase in tends to keep the ratio low.

#### 5.2.2 Adaptive response to positive energy balance

In the Maas & Smith (2006) model, the energy imbalance is imposed from without, and the response is an adaptive increase of , which allows the adipose tissues to absorb the energy surplus. However, it is now well established that signals emanating from adipocytes (adipokines) play a central role (Kennedy 1953; Benoit *et al.* 2004). A positive imbalance between intake and expenditure can be redressed by reducing the intake, increasing the expenditure or both. Food intake, basal metabolic rate and thermogenesis are centrally coordinated in the hypothalamus; figure 9 depicts major signalling pathways. The main sensory input centre is the arcuate nucleus (ARC), which contains neurons that are sensitive to the plasma levels of insulin, glucose, ghrelin and the adipokine leptin. These neurons modulate food intake as well as the thyroid hormone axis (Chen *et al.* 1993; Spanswick *et al.* 1997, 2000; Kim *et al.* 2000; Barsh & Schwartz 2002; Dobbins *et al.* 2002; Nowak *et al.* 2002; van den Top *et al.* 2004).

The hypothalamus also regulates muscle glycogen synthesis through the autonomous nervous system (ANS; Perrin *et al.* 2004). As noted above, extra-hepatic glycogen stores can serve as buffers that temporarily absorb the energy excess (in particular, its carbohydrate component) while preadipocyte proliferation takes place. Leptin stimulation (via the ANS) of glucose transport into extra-hepatic stores may also be involved (Kamohara *et al.* 1997; Haque *et al.* 1999; Shiuchi *et al.* 2001). In hyperleptinaemic rats, adipose TAG reserves all but disappear, in contrast to pair-fed controls, while energy intake and body weight were the same in these controls (Chen *et al.* 1993), suggesting that energy expenditure must have been nearly equal in the two groups (intracerebroventricular infusion of leptin results in increased expenditure, but the effect is much stronger after 3 days than after 14 days; Halaas *et al.* 1997). Energy reserves at the end of the experiment must therefore have been comparable despite the marked difference in adiposity, pointing to the involvement of glycogen stores—particularly extra-hepatic stores, given that liver weight did not change in similar experiments on mice (Levin *et al.* 1996).

To account for these centrally coordinated compensatory changes in energy intake and expenditure, the Maas & Smith model can be extended as shown in figure 10. The plasma leptin level is included as a state variable and treated as the main adiposity signal driving the central response (Benoit *et al.* 2004). The central nervous system (CNS) modulates food intake and energy expenditure. The long-term behaviour of this modified Maas & Smith model depends on the mathematical specification of the secretion term in the leptin kinetics, as well as on the specification of the input/output behaviour of the CNS. It is worthwhile to discuss this question of formulating an appropriate model in some depth, since it is connected with several very topical issues in obesity research.

##### Leptin as an adiposity signal

Leptin is secreted by adipocytes. The leptin plasma concentration correlates well with percentage body fat (Considine *et al.* 1996; Friedman & Halaas 1998), and is thus generally regarded as an adiposity signal (Benoit *et al.* 2004). The long-term average is proportional to the leptin secretion rate, which is a product of and the secretion rate per adipocyte. If the latter is a basal value, independent of adipocyte size, we would have , where is the plasma volume. By contrast, we obtain if the leptin secretion rate is proportional to the size of the adipocyte (which is on average). Again, if each adipocyte signals in proportion to the assimilatory flux it is conducting, we would obtain . Which, if any, of these options is the most realistic model?

Leptin gene expression is regulated by the hexosamine pathway that converts fructose 6-phosphate (Fruc6P) into UDP-*N*-acetylglucosamine (UDP-GlcNac); this product donates GlcNac moieties to transcription factors thus promoting the transcription of leptin mRNA (Wang *et al.* 1998). The flux into the hexosamine pathway is increased when Fruc6P accumulates due to NEFA availability (which slows down glycolysis) and/or hyperglycaemia (which promotes the influx of glucose). Larger adipocytes express more leptin (Maffei *et al.* 1995). Leptin expression is also regulated at the posttranslational level via the mTOR-mediated pathway, which is activated by free amino acids (Roh *et al.* 2003). The secretion of leptin is regulated by intracellular ATP (Levy *et al.* 2000). In keeping with these nutrient-sensing regulatory mechanisms, plasma leptin levels decrease during dynamic weight loss (Rosenbaum *et al.* 1997; Velkoska *et al.* 2003), whereas leptin rises during hyperinsulinaemic clamp (Boden *et al.* 1997) and following food intake (more markedly so in obese subjects). Plasma levels slowly fall during sleep (Yildiz *et al.* 2004), yielding a diurnal rhythm that has been shown to be entrained to meal timing rather than an endogenous clock (Schoeller *et al.* 1997). Collectively, these observations suggest, in general, that leptin secretion by adipocytes reflects assimilatory activity in these cells and, in particular, that adipocytes operating near their maximum storage capacity have high levels of leptin expression.

#### 5.2.3 The Friedman adipostasis model

Friedman & Halaas (1998) proposed a conceptual model in which the CNS permits a positive energy balance to persist while remains well below a critical value that they call the *leptin set point*. As leptin levels start to approach this threshold, the modulatory effects on food intake and energy expenditure come into effect, so that a neutral energy balance is attained at plasma leptin levels in the vicinity of the set point.

The leptin set point in this model is equivalent to an adiposity set point provided that we have . Thus, the model explains the existence of an apparent ‘body weight set point’ particular to the individual and subject to genetic variation in the population (Keesey & Hirvonen 1997). Furthermore, the Friedman model predicts that obesity can be caused by a reduced ability of adipocytes to produce leptin, or an impairment of the transport of leptin to the brain; in either case the effect will be that the brain ‘perceives’ adipose reserves to be smaller than they are. Alternatively, obesity can be caused by a reduced sensitivity to leptin in the CNS. In the first scenario, obesity with normal leptin levels ensues, whereas leptin levels will be elevated in the second scenario.

To represent the Friedman model mathematically, we start with the following expression for the energy balance:(5.5)where * x* represents an internal state of the hypothalamus; the assimilatory influx of energy; is the expenditure of energy (expressed in the same units as

*Q*

_{TAG}to keep the notation as simple as possible); the multiplier represents the hypothalamic modulation of food intake;

*represents environmental factors such as food availability that codetermine the assimilatory flux; and*

**u***represents physiological or endocrinological variables that affect energy expenditure, such as the body growth rate, pregnancy and physical exercise.*

**v**Suppose now that the hypothalamus receives inputs carrying information about (i) *Q*_{TAG} and/or , which is reasonable in view of the foregoing discussion, as well as (ii) * v*. Then, formally,

*can be recovered from these inputs; say where*

**u***is at least locally uniquely defined by the balance equation (5.5). This shows that the hypothalamus can regulate*

**Υ***Q*

_{TAG}towards an ‘adiposity set point’, which can be, in general, any function of

*and/or*

**u***. This is the*

**v***adipostasis property*. Two caveats apply: first, this regulatory ability is limited by obvious physiological constraints, namely the maximum food intake rate and the minimum expenditure rate. Second, it should not be thought that the hypothalamus is supposed to solve algebraic equations. Our formal representation is simply intended to highlight the fact that the information inputs mentioned suffice to obtain the adipostasis property.

The above argument is fairly abstract; let us consider two simple examples. First, take * x* to be one-dimensional with first-order relaxation towards the leptin input(5.6)where and are positive constants, and let the food intake multiplier

*η*be a steeply descending Hill-type function of

*x*(cf. Velkoska

*et al.*2003) with midpoint parameter . Then, as long as the quantity fluctuates in the interval [0,1], adiposity fluctuates around . The latter is the apparent set point. This pattern is typical: apparent set points arise as a compound of rate constants and sensitivity (coupling) constants.

Although very simple, this minimal implementation of the Friedman hypothesis already shows how one can make sense of Keesey & Hirvonen's (1997) remarks about the body weight set point and variability of this set point within the population. The model states that individuals with a lower-than-average -value have a higher apparent adiposity set point. The coupling constant represents hypothalamic sensitivity to leptin. Thus, if the hypothalamic leptin receptor is knocked out, would be zero and the model predicts that the individual becomes massively obese (i.e. goes to its maximum ), in keeping with the findings of Maffei *et al.* (1995).

In the simple model, the adiposity set point is , a constant independent of ambient conditions (* u*) such as food availability. Alternatively, one might propose that the set point proportional to the time-average caloric density of the available diet(5.7)where

*θ*and

*ν*are empirical constants. This is an intuitively attractive and physiologically plausible model, at least qualitatively. The caloric-density set point model gives rise to first-order adipose density dynamics in times of positive energy balance (see Sousa

*et al.*(2008) for a discussion of similar dynamics). The model can be justified on the basis of the adipostasis property. For example, the simple Friedman model can be extended by assuming that the midpoint parameter

*ξ*is governed by the available inputs listed above. We can recover the functional dependence by solving the following equation for

*ξ*:(5.8)Again, there is no suggestion that the hypothalamus explicitly carries out the algebra; the point is to show how it can nonetheless affect the equivalent mapping by showing how we can represent this mapping in terms of our notation.

### 5.3 Negative energy balance and starvation

The Friedman model deals with the adaptation to positive energy balance. There are also adaptations to negative energy balance. This response consists of coordinated regulation of energy expenditure, the mobilization of energy stores and hyperphagia once food becomes available. Both hypothalamus and peripheral tissues partake in this response (Keesey & Hirvonen 1997; Frayn 2003).

Energy expenditure is reduced as leptin levels fall; peripherally, leptin affects insulin secretion (Lam *et al.* 2004) and glucose usage (Haque *et al.* 1999; Shiuchi *et al.* 2001) while centrally, there is a shutdown of the leptin–TRH–TSH–thyroid axis (cf. figure 9; Rosenbaum *et al.* 1997; Kim *et al.* 2000) as well as the leptin–growth hormone axis (Tannenbaum *et al.* 1998). The shift in the insulin/glucagon ratio stimulates the Cori cycle (§4.2).

The key problem in the mobilization of energy stores is the coordinated use of the three types of energy reserve, *viz.* glycogen, fat and muscle protein. Hepatic glycogen stores are rapidly depleted (Nilsson & Hultman 1973), leaving adipose stores and protein. Various adaptations favour the use of the former: gluconeogenesis from lactate via the Cori cycle allows obligatory glycolytic tissues to use energy derived from fatty acids, which are also converted into ketone bodies used by the CNS as an alternative to glucose, whereas the residual obligatory usage of glucose is sustained by protein breakdown via the alanine cycle (§4.2; Salway 2004). The CRH–ACTH–cortisol axis regulates proteolysis in muscle, counteracting the effects of insulin and IGF-1 (Brook & Marshall 2001); this axis governs the terminal phase of starvation, when adipose stores have run out.

The long-term dynamics of *Q*_{TAG} is thus governed by (minimized) obligatory energy requirements, and the phenomenological approach that led to equation (5.7) may therefore be inappropriate, although Kooijman (2000) argued that this equation can be used to describe starvation as well as energy surplus.

## 6. Conclusion and outlook

Since the purpose of this paper was to survey the existing literature, we have followed the literature's glucocentric bias, starting with basic models of glucostasis and the various ways in which these models are extended with additional dynamic degrees of freedom. However, the problems presented by diabetes, obesity and controlled nutrition of, for example, cachexic subjects call for models that integrate glucostasis with adipostasis and myostasis at the long-term time scale (the ‘developmental’ or ‘life-history’ time scale). We have focused attention here on the interplay of glucostat and adipostat. Other aspects that are important on the developmental time scale, but which have not been emphasized in this review, include interactions with lipoprotein metabolism and cholesterol homeostasis and the interplay between amino acid metabolism, glucoplastic carbon metabolism, myostasis and the regulation of the overall growth rate.

An open question is how one can include all these aspects in one overarching model, and, indeed, whether one should. The 16-state variable model shown in figure 1 is already prohibitive in terms of parameter identifiability, and could easily be extended with many more humoral factors and compartments. Perhaps the approach taken in the analysis of the Friedman model and the adipostasis property (§5.2.3) will prove to be fruitful in general. Thus, to describe dynamics on the life-history time scale, one starts with balance equations, augmented with a minimum of additional physiological state variables needed to describe the regulation of the flux terms in response to, for example, dietary challenges. What is as yet missing is a systematic way of deriving the dynamics of these additional state variables at the slow time scale from the underlying detailed dynamics at the fast physiological time scales.

## Acknowledgments

R.P. gratefully acknowledges the support of the Thai government and the MOAC Doctoral Training Centre at Warwick University. H.A.B. thanks the EPSRC and BioSim.

## Footnotes

- Received May 20, 2008.
- Accepted June 10, 2008.

- © 2008 The Royal Society