See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/12571934

"Mathematical Modelling of the Intravenous Glucose

Tolerance Test"

Article in Journal of Mathematical Biology · March 2000

DOI: 10.1007/s002850050007 · Source: PubMed

CITATIONS READS
362 2,802

2 authors, including:

Andrea De Gaetano
Italian National Research Council and Óbuda University
252 PUBLICATIONS 10,920 CITATIONS

SEE PROFILE

All content following this page was uploaded by Andrea De Gaetano on 05 June 2014.

The user has requested enhancement of the downloaded file.

J. Math. Biol. (2000) 40: 136}168

Mathematical modelling of the intravenous glucose

tolerance test
Andrea De Gaetano 夹, Ovide Arino
CNR, Centro Fisiopatologia Shock, Università Cattolica del Sacro Cuore,
Largo A. Gemelli, 8, 00168 Roma, Italy. e-mail: [email protected]
Laboratoire de MatheH matiques AppliqueH es, UPRES A 5033 IPRA, Pau,
France

Received: 22 June 1998 / Revised version: 24 February 1999

Abstract. Several attempts at building a satisfactory model of the

glucose-insulin system are recorded in the literature. The minimal
model, which is the model currently mostly used in physiological
research on the metabolism of glucose, was proposed in the early
eighties for the interpretation of the glucose and insulin plasma con-
centrations following the intravenous glucose tolerance test. It is com-
posed of two parts: the "rst consists of two di!erential equations and
describes the glucose plasma concentration time-course treating insu-
lin plasma concentration as a known forcing function; the second
consists of a single equation and describes the time course of plasma
insulin concentration treating glucose plasma concentration as
a known forcing function. The two parts are to be separately estimated
on the available data. In order to study glucose-insulin homeostasis as
a single dynamical system, a uni"ed model would be desirable. To this
end, the simple coupling of the original two parts of the minimal model
is not appropriate, since it can be shown that, for commonly observed
combinations of parameter values, the coupled model would not admit
an equilibrium and the concentration of active insulin in the &&distant''
compartment would be predicted to increase without bounds. For
comparison, a simple delay-di!erential model is introduced, is demon-
strated to be globally asymptotically stable around a unique equilib-
rium point corresponding to the pre-bolus conditions, and is shown to
have positive and bounded solutions for all times. The results of "tting
the delay-di!erential model to experimental data from ten healthy

*****
夹
Corresponding author.
Mathematical modelling of IVGTT 137

volunteers are also shown. It is concluded that a global uni"ed model

is both theoretically desirable and practically usable, and that any
such model ought to undergo formal analysis to establish its appro-
priateness and to exclude con#icts with accepted physiological
notions.

Key words: Glucose } Insulin } Dynamical systems } Di!erential

equations } Delay

1. Introduction

Since the sixties, the homeostasis of glucose, involving the secretion of

its controlling hormone insulin by the pancreas, has been the object of
several mathematical models [1, 3, 6, 12, 14, 16, 17, 23, 26, 27, 29]. With
the increased emphasis on derangements of the sensitivity of tissues to
insulin in diverse pathological conditions like diabetes, obesity, and
cardiovascular disease [8, 9, 11], the quantitation of insulin sensitivity
from a relatively non-invasive test procedure has acquired greater
importance in physiological research.
Of the two experimental procedures currently in use for the estima-
tion of insulin sensitivity in a subject, the intra venous glucose toler-
ance test (IVGTT) has appeal, over the euglycemic hyperinsulinemic
clamp [10], because of its ease of execution, especially in its unmodi"ed
form (without additional Tolbutamide injection [30]), and because of
the richness of the information its analysis yields. The test consists of
injecting I.V. a bolus of glucose and frequently sampling the glucose
and insulin plasma concentrations afterwards, for a period of about
three hours.
The physiological model which has been mostly used in the inter-
pretation of the IVGTT, since its publication in the early eighties, is
generally known as the &&Minimal Model'' [2, 28]: it is described below
(Eqs. 1}3). The model, as originally proposed by the authors, is to be
regarded as composed of two separate parts. The "rst part [2] uses
Equations (1) and (2) to describe the time course of plasma glucose
concentration, accounting for the dynamics of glucose uptake depen-
dent on and independent of circulating insulin; for this "rst part,
plasma insulin concentration is to be regarded as a known forcing
function. The second part [28] consists of Equation (3) and describes
the time course of plasma insulin concentration, accounting for the
dynamics of pancreatic insulin release in response to the glucose
stimulus; for this second part plasma glucose concentration is to be
138 A. De Gaetano, O. Arino

regarded as a known forcing function. The proposing authors speci"-

cally stated [19] that the model parameter "tting has to be conducted
in two steps: "rst, using the recorded insulin concentration as input
data in order to derive the parameters in the "rst two equations, then
using the recorded glucose as input data to derive the parameters in the
third equation.
However, the glucose-insulin system is an integrated physiologic
dynamical system and we would like to be able to describe it as
a whole, this uni"ed description being also conducive to a single-step
parameter "tting process.
We "rst o!er a formal study of the dynamical system obtained by
coupling the two parts of the minimal model. By doing so, this coupled
model will be shown to present some di$culties in describing math-
ematically the glucose-insulin system from a uni"ed point of view. We
want to underscore that by so doing we depart from the explicitly
stated objectives of the model's proponents: no criticism is therefore
implied regarding its ful"lling its original goals and regarding its
practical usefulness; indeed, our group uses the minimal model in the
routine evaluation of insulin sensitivity in clinical patients [5], and we
have introduced improved statistical techniques to make the deter-
mination of its coe$cients (notably those of its "rst part, relative to
glucose dynamics) more e$cient [7].
We then introduce a di!erent global dynamical model of the
glucose-insulin system and discuss its stability properties. We also
show some numerical results obtained from "tting this dynamical
model to plasma glucose and insulin concentrations measured on
healthy volunteers undergoing the IVGTT.

2. Materials and methods

2.1. Patient sample

Ten healthy volunteers (5 males and 5 females, anthropometric charac-

teristics reported in Table 1) participated in the study. All subjects had
negative family and personal histories for diabetes mellitus and other
endocrine diseases, were on no medications, had no current illness and
had maintained a constant body weight for the six months preceding
the study. For the three days preceding the study each subject followed
a standard composition diet (55% carbohydrate, 30% fat, 15% pro-
tein) ad libitum with at least 250g carbohydrates per day.
Each study was performed at 8:00 AM, after an overnight fast, with
the subject supine in a quiet room with constant temperature of
Mathematical modelling of IVGTT 139

22}24 3C. Bilateral polyethylene IV cannulas were inserted into ante-

cubital veins. The standard Intra Venous Glucose Tolerance Test
(IVGTT) was employed, without using additional Tolbutamide so as
to be able to use the recorded data for pancreatic secretion evaluation:
at time 0 (0) a 33% glucose solution (0.33 g Glucose/kg body weight)
was rapidly injected (less than 3 minutes) through one arm line. Blood
samples (3 ml each, in lithium heparin) were obtained at !30, !15,
0, 2, 4, 6, 8, 10, 12, 15, 20, 25, 30, 35, 40, 50, 60, 80, 100, 120,
140, 160 and 180 through the contralateral arm vein. Each sample
was immediately centrifuged and plasma was separated. Plasma glu-
cose was measured by the glucose oxidase method (Beckman Glucose
Analyzer II, Beckman Instruments, Fullerton, CA, USA). Plasma insu-
lin was assayed by standard radio immunoassay technique. The
plasma levels of glucose and insulin obtained at !30, !15 and 0
were averaged to yield the baseline values referred to 0.

2.2. Minimal model

As we have seen, the physiologic experiment consists in injecting into

the bloodstream of the experimental subject a bolus of glucose, thus
inducing an (impulsive) increase in the plasma glucose concentration
G(t) and a corresponding increase of the plasma concentration of
insulin I(t), secreted by the pancreas. These concentrations are mea-
sured during a three-hour time interval beginning at injection, after
which time interval it is found that the perturbed concentrations G(t)
and I(t) have essentially returned to normal. In order to describe the
time course of these concentrations, the minimal model of the glucose-
insulin kinetics has been proposed.
We will use the standard formulation of the minimal model,
renaming some parameters for ease of notation:
dG (t)
"![p #X(t)] G(t)#p G , G (0)"p (1)
dt   @ 
dX (t)
"!p X (t)#p [I(t)!I ], X (0)"0 (2)
dt   @
dI (t)
"p [G(t)!p ]>t!p [I (t)!I ], I (0)"p #I (3)
dt    @  @
where
G (t) [mg/dl] is the blood glucose concentration at time t [min];
I(t) [lUI/ml] is the blood insulin concentration;
140 A. De Gaetano, O. Arino

X (t) [min\ ] is an auxiliary function representing insulin-excitable

tissue glucose uptake activity, proportional to insulin concentra-
tion in a &&distant'' compartment;
G [mg/dl] is the subject's baseline glycemia;
@
I [lUI/ml] is the subject's baseline insulinemia;
@
p [mg/dl] is the theoretical glycemia at time 0 after the instan-

taneous glucose bolus;
p [min\] is the glucose &&mass action'' rate constant, i.e. the insu-

lin-independent rate constant of tissue glucose uptake, &&glucose
e!ectiveness'';
p [min\] is the rate constant expressing the spontaneous decrease

of tissue glucose uptake ability;
p [min\ (lUI/ml)\] is the insulin-dependent increase in tissue

glucose uptake ability, per unit of insulin concentration excess
over baseline insulin;
p [(lUI/ml) (mg/dl)\ min\] is the rate of pancreatic release of

insulin after the bolus, per minute and per mg/dl of glucose
concentration above the &&target'' glycemia;
p [mg/dl] is the pancreatic &&target glycemia'';

p [min\] is the "rst order decay rate constant for Insulin in

plasma;
p [lUI/ml] is the theoretical plasma insulin concentration at time

0, above basal insulinemia, immediately after the glucose bolus.

In Eq. (3), only the positive part of the term [G(t)!p ] is taken, i.e.

when G (t ) is greater than p the value is taken to be [G(t)!p ],
 
otherwise the value is zero. Also in Eq. (3), the multiplication by t is
introduced by the authors to express, as a "rst approximation, the
hypothesis that the e!ect of circulating hyperglycemia on the rate of
pancreatic secretion of insulin is proportional both to the hyper-
glycemia attained and to the time elapsed from the glucose stimulus
[28]. Multiplying by t in this way introduces the necessity of establish-
ing an origin for time, e!ectively binding this model to the IVGTT
experimental procedure.
Parameters p , p , p , p , p and p are usually referred to in the
     
literature as G , S , c, h, n and I respectively, while the insulin
 % 
sensitivity index S is computed as p /p .
'  

2.3. Dynamical model

In order to overcome the di$culties of the coupled minimal

model, another model for the glucose-insulin system is proposed. The
Mathematical modelling of IVGTT 141

physiological hypotheses underlying Eq. (1) in the minimal model

above have been retained, i.e. that disappearance of glucose from
plasma may be described as a "rst-order process, of rate partly depen-
dent on insulin concentration and partly independent of it. However, it
was desired not to introduce non-observable state variables if possible,
to represent delay explicitly when a delay is apparent, and to avoid
a non-autonomous form of the equation, di$cult to justify from
a physiological point of view and likely to contribute to model instabil-
ity. The questionable physiological assumption that the pancreas is
able to linearly increase its rate of insulin secretion with time, and the
related necessity of establishing an initial time point with respect to
which all biochemical events take place, are both avoided in the
proposed formulation. In this way an attempt is made to formulate
a model embodying the underlying physiological mechanism, without
associating it by necessity to the IVGTT experiment. The dynamical
model of the glucose-insulin system to be studied is therefore:
dG(t)
"!b G (t)!b I (t) G(t)#b ,
dt   
G(t),G ∀t3[!b , 0), G(0)"G #b (4)
@  @ 


dI (t) b R
"!b I (t)#  G (s) ds, I (0)"I #b b , (5)
dt  b t!b @  

where
t [min] is time;
G [mg/dl] is the glucose plasma concentration;
G [mg/dl] is the basal (preinjection) plasma glucose concentration;
@
I [pM] is the insulin plasma concentration;
I [pM] is the basal (preinjection) insulin plasma concentration;
@
b [mg/dl] is the theoretical increase in plasma concentration over

basal glucose concentration at time zero after instantaneous ad-
ministration and redistribution of the I.V. glucose bolus;
b [min\] is the spontaneous glucose "rst order disappearance rate

constant;
b [min\] is the apparent "rst-order disappearance rate constant for

insulin;
b [pM/(mg/dl)] is the "rst-phase insulin concentration increase per

(mg/dl) increase in the concentration of glucose at time zero due to
the injected bolus;
b [min\ pM\] is the constant amount of insulin-dependent

glucose disappearance rate constant per pM of plasma insulin
concentration;
142 A. De Gaetano, O. Arino

b [min] is the length of the past period whose plasma glucose


concentrations in#uence the current pancreatic insulin secretion;
b [min\ pM/(mg/dl)] is the constant amount of second-phase

insulin release rate per (mg/dl) of average plasma glucose con-
centration throughout the previous b minutes;

b [(mg/dl) min\] is the constant increase in plasma glucose concen-

tration due to constant baseline liver glucose release.

The above model describes glucose concentration changes in blood

as depending on spontaneous, insulin-independent net glucose tissue
uptake, on insulin-dependent net glucose tissue uptake and on con-
stant baseline liver glucose production. The term net glucose uptake
indicates that changes in tissue glucose uptake and in liver glucose
delivery are considered together.
Insulin plasma concentration changes are considered to depend on
a spontaneous constant-rate decay, due to insulin catabolism, and on
pancreatic insulin secretion. The delay term refers to the pancreatic
secretion of insulin: e!ective pancreatic secretion (after the liver "rst-
pass e!ect) at time t is considered to be proportional to the average
value of glucose concentration in the b minutes preceding time t.

Due to the delay, as initial conditions for the problem we have to
specify not only the level of glucose at time zero, but also its value at
times from !b to 0.

If < [ml/kgBW] is the volume of distribution of glucose, = [kg]
%
the weight of the experimental subject and D [mg] is the dose of
%
injected glucose, then b "D /(< =), and the model may be expressed
 % %
in terms of < instead terms of b .
% 
The term (1/b ) in front of the integral in Eq. (5) has been introduc-

ed so as to make the integral equal one for constant unit glucose
concentration, thus making b , pancreatic responsiveness, independent

of b , the time period of pancreatic sensitivity to plasma glucose

concentrations.
The free parameters are only six (b through b ). In fact, assuming
 
the subject is at equilibrium at (G , I ) for a su$ciently long time
@ @
('b ) prior to the administration of the bolus, then

0"!b G !b I G #b and 0"!b I #b G together imply
 @  @ @   @  @
I
b "b G #b I G , b "b @ .
  @  @ @  G
@
For model "tting, observations have been weighted according to
the usual IVGTT scheme [19]: glucose observations before 8 minutes
have been given a weight of zero (assuming glucose bolus distribution
Mathematical modelling of IVGTT 143

to be complete by 8 minutes), and insulin observations before the "rst

insulin peak (usually at 2 or 4 minutes) have also been given a weight of
zero. No points have been overweighted (i.e. all points have weight
either zero or one). Parameter values were obtained by weighted least
squares using a quasi-Newton minimization algorithm [22].

3. Results

The proofs of all the statements contained in this section are reported
in the Appendix.

I. Formal study of the minimal model

Recall that in the minimal model p is the target glycemia which the

pancreatic secretion of insulin attempts to attain (i.e. above which the
pancreas is assumed to secrete the glycemia-lowering hormone insu-
lin), whereas G is the measured baseline glycemia, which results from
@
the equilibrium between the pancreatic action to lower glycemia down
to p and the endogenous (liver) glucose production which tends to

raise glycemia. In general, G may be greater than p , and this is in fact
@ 
the case in [19] where the program to estimate the parameters of the
minimal model is described: in the reported series G "92 mg/dl,
@
p "89.5 mg/dl.

Proposition I.1 refers to the case in which glycemia G (t) returns to
basal G after the metabolization of the glucose bolus.
@
Proposition I.1. Suppose G 'p , lim sup G (t)'p . ¹hen
@  R 
lim sup X (t)"R.
R
Proposition I.2 refers to the case in which glycemia tends to drop
below the pancreatic target level p :

Proposition I.2. Suppose lim sup G(t)(p ; then G 6p .
R  @ 
Remark I.3. Proposition I.1 shows the model to exhibit a pathologic
behavior if lim sup G(t)'p . Proposition I.2 shows that if
R 
lim sup G(t)(p , then necessarily G 6p , contrary to what ob-
R  @ 
served. The two propositions above leave open the possibility that the
upper limit of glycemia, as t increases, exactly equals p . The analysis

of this boundary case is not easy: however, either one of the two
above cases would be produced for arbitrarily small changes in the p

parameter. In fact, more can be proven, i.e.:
144 A. De Gaetano, O. Arino

Proposition I.4. For any value p (G the system does not admit an
 @
equilibrium.

The above fact may be stated in another way:

Proposition I.5. If the subject is considered to be at steady state before

the glucose bolus, then G must be lesser than or equal to p .
@ 
The above objections to the qualitative behavior of the system (1, 2,
3) would indicate as a possible solution the imposition that the param-
eter p exactly equal the experimentally observed quantity G . This
 @
empirical operation, which would actually change the model to a dif-
ferent one, "nds little justi"cation. In fact, p is the unknown but true

value of a model parameter, while G is a measured quantity equal to
@
the sum of the true unknown value of the pre-injection equilibrium
state plus some (observation) error. In any case, it is interesting to note
that, in the case where p "G , the possible solutions, which start out
 @
at a greater value than G (due to the bolus injection of glucose) are
@
forced to pass below the G level before they can converge to G , and
@ @
once they pass below this value, they can never cross it again to become
greater than G . In other words, it is not possible for a solution to
@
converge to G from above or to oscillate in any way (damped or
@
otherwise) around G .
@
Proposition I.6. ¸et p '0, p "G , G(0)'G . Assume G(t), X(t)
  @ @
bounded. ¹hen there exists ¹'0 such that G(t)'G for all t(¹,
@
G(¹)"G , and G (t)(G for all t'¹.
@ @
A brief discussion of the properties of the minimal model in case
G (p follows the proof of Proposition I.6 in the Appendix.
@ 

II. Stability of the dynamical model

It can be shown that the dynamical model admits one and only one
equilibrium point with positive concentrations, (G , I ).
@ @
Proposition II.3. ¹he solutions +G (t), I (t), are positive and bounded.

Proposition II.4. ¹he time derivatives of the solutions are bounded.

It can be "nally shown that any solution to the original

system converges to (G , I ) and thus this state is asymptotically
@ @
stable.
Mathematical modelling of IVGTT 145

III. Numerical xtting of the dynamical model

The parameter estimates obtained for each subject by "tting the

dynamical model to the experimental data are reported in Table 1,
which also shows the sample mean, sample standard deviation, sample
standard error and coe$cient of variation for each parameter. The
corresponding sample correlation matrix is reported in Table 2.
It has to be noted that the parameter b was evaluated in every

experimental subject as an integer value, since the resolution of the
estimate of delay (in minutes) depends on the smallest increment in
time used in the evaluation of the integral, in this case one minute.
The "tting was good to excellent for all subjects: the R-squared
computed from the simple (unweighted) deviance lying between 0.862
and 0.990.
Figures 1, 2 and 3 report the observed and predicted time courses
for three prototypical subjects: one (Fig. 1) showing an apparently
exponential decay of insulin after its primary peak, one (Fig. 2) presen-
ting no secondary insulin peak but a clear non-exponential decay of
insulin towards the basal value, and "nally one (Fig. 3) showing
a well-de"ned secondary insulin peak.
We notice in passing that while the minimal model was originally
described using lUI/ml as the units for Insulin concentration, our
laboratory provides insulin concentrations in pM (pico moles/Liter):
1 lUI/ml being equal to 7.00 pM.

4. Discussion

The IVGTT is a reasonably simple experiment yielding a potentially

very rich data set, and, in the applications, parameters referring both to
tissue glucose uptake and to pancreatic responsiveness are of interest
[4, 20, 21].
From a dynamical point of view, the pancreas and tissues form an
integrated system with feedback regulations and it would seem desir-
able to have a model explicitly representing the whole system, which
could be "tted in a single pass to both glucose and insulin data, rather
than splitting the model into two subsystems and "tting separately
each one. In fact, for a model "tting simultaneously the two arms of the
control mechanism, the error variance would be a more appropriate
expression of the e!ective applicability of the assumptions underlying
both subsystems to the experimental situation. In other words, by
"tting simultaneously the two parts, we observe the coherent "tting of
the entire dynamical model to the entire set of observations. By
 146

Table 1. Dynamical model results: anthropometric characteristics and dynamical model parameter values found for each experimental subject, together with their
sample mean, standard deviation, standard error and coe$cient of variation. BW is body weight, LBM is lean body mass, FBM is fat mass, bas. gluc. is basal blood
glucose, bas. insul. is basal plasma insulin, R is the (unweighted data) coe$cient of determination.

Subject Sex Age Height BW LBM FM bas. bas. b b b b b b b b R

       
(years) (cm) (kg) (kg) (kg) gluc. insul. (mg/dl) (min\) (min\) pM/ (min\) (min) min\ pM/ (mg/dl)
(mg/dl) (pM) (mg/dl) pM\ (mg/dl) min\

1 m 35 172 72 56 16 69 71.3 170 0.0226 0.0437 2.57 3.80E-08 20 0.045 1.56 0.865
2 f 28 155 45 36.2 8.8 79 51.7 241 0.0509 0.2062 3.55 1.29E-07 14 0.135 4.02 0.955
3 f 25 162 61 48.4 12.6 74 29.4 208 0.0309 0.1817 2.96 6.99E-07 12 0.072 2.29 0.931
4 m 32 169 68 53.5 14.5 80 56.6 355 0.0084 0.1039 4.25 7.55E-05 8 0.073 1.01 0.985
5 m 23 179 65 55 10 74 45 216 0.0273 0.0275 2.77 1.10E-07 5 0.017 2.02 0.869
6 f 27 162 65 44.5 20.5 88 68.6 209 0.0002 0.0422 1.64 1.09E-04 23 0.033 0.68 0.953
7 m 25 170 66 53 13 87 37.9 311 0.0001 0.2196 0.64 3.73E-04 23 0.096 1.24 0.957
8 f 34 158 64 42.4 21.6 78 55.8 217 0.0565 0.0438 4.39 5.70E-06 19 0.031 4.43 0.99
9 m 42 172 78 61.2 16.8 70 43.8 156 0.0135 0.2972 5.92 3.51E-08 11 0.186 0.94 0.93
10 f 55 169 67 47.4 19.6 67 37.7 184 0.0159 0.0965 2.51 8.72E-08 14 0.054 1.07 0.862

Mean 32.6 166.8 65.1 49.8 15.3 76.6 49.8 226.7 0.0226 0.1262 3.12 5.64E-05 14.9 0.074 1.93 0.93
s.d. 9.8 7.3 8.5 7.4 4.4 7.2 13.6 62.1 0.0194 0.0938 1.49 1.18E-04 6.2 0.052 1.31 0.048
s.e. 3.1 2.3 2.7 2.3 1.4 2.3 4.3 19.6 0.0061 0.0297 0.47 3.73E-05 2 0.017 0.41 0.015
c.v. (%) 9.5 1.4 4.1 4.7 9 3 8.6 8.7 27.1 23.5 15.1 66.1 13.1 22.4 21.5 1.6
A. De Gaetano, O. Arino
Mathematical modelling of IVGTT 147

Table 2. Dynamical model: matrix of sample partial correlation coe$cients between

the free parameters.

b b b b b b
     
b 1.000

b !0.233 1.000

b 0.044 !0.110 1.000

b !0.213 0.398 0.290 1.000

b 0.594 !0.571 0.237 !0.628 1.000

b !0.062 !0.168 !0.081 !0.532 0.519 1.000


Fig. 1. Blood glucose (data: blank squares; model: continuous line) and plasma insulin
(data: solid diamonds; model: dashed line) time courses for subject 8. Both glucose and
insulin decays seem grossly exponential.

splitting the two subsystems we may estimate coe$cients for one

segment, by optimally "tting the relative data, which are not the ones
which would generate an optimal approximation to the whole data set
if the interaction between the two subsystems were allowed. The end
result is that by splitting the system we obtain an impression of success
because our error looks smaller, but we are in fact omitting an internal
coherency check. Consequently, we may assert that the (global) system
works in some way (i.e. with some parameters) which may be di!erent
from the best approximation to the real one. It may indeed happen
(like in "nding a value for b smaller than the value for G ) that
 @
the global system cannot work at all with the estimated parameter
values.
148 A. De Gaetano, O. Arino

Fig. 2. Blood glucose (data: blank squares; model: continuous line) and plasma insulin
(data: solid diamonds; model: dashed line) time courses for subject 6. While glucose
seems grossly exponentially decaying, the insulin time course shows appreciable
changes in decay rate, but no secondary peak.

Fig. 3. Blood glucose (data: blank squares; model: continuous line) and plasma insulin
(data: solid diamonds; model: dashed line) time courses for subject 7. To an evident
secondary peak in insulin plasma concentration, there corresponds a parallel variation
in blood glucose decay rate.

In the applied use of the minimal model, the feeling that the most
important piece of information is the S index may lead physicians to
'
"t only the "rst part of the minimal model (Eqs. 1, 2), disregarding
the second part (Eq. 3), using insulin as the &&"xed'' input data [18]:
using the minimal model in this way avoids formal problems due to
Mathematical modelling of IVGTT 149

feedback, but it fails to provide information on the pancreatic response

to circulating glucose, information which is implicitly contained in the
available measured insulin concentration curve. It is to be noted that
insu$cient interaction between insulin and glucose levels (e.g. in the
case of NIDDM or obese patients) makes estimation of the minimal
model parameters impossible in a certain number of cases. Taking into
account the absolute levels of glucose and insulin reached during the
experiment may help: an insulin resistant subject being de"ned as one
whose insulin levels are maintained higher than normal for compara-
ble glucose levels. It seems, in fact, that NIDDM subjects have a higher
insulin response, and a more prolonged one, as if the stimulus repre-
sented by glucose were higher in them than in normals. Following this
line of thought it would seem that it is the second part of the minimal
model (Eq. 3) that probably best discriminates this behavior, as if
insulin sensitivity were only part of the picture, pancreatic responsive-
ness being the other part: in other words, estimation of the parameters
of (Eq. 3) or of an equivalent formulation of glucose e!ect on insulin
seems needed as well.
Once the necessity of a single global dynamical model of the
insulin-glucose system is appreciated, the "rst likely candidate model
to be studied is the one obtained by coupling the three Minimal Model
equations (Eqs. 1}3). In the present work, a formal study of the
qualitative behavior of this model has been conducted, and three basic
results have been proved. The "rst result (Propositions I.1 and I.2) says
that under the stated set of model de"nitions, allowing baseline
glycemia to be greater than pancreatic target glycemia, the model
equations give rise to unbounded solutions except, at best, for a bound-
ary case (any small perturbation from which would reinstate the
unbounded solutions). The second result (Propositions I.4 and I.5)
says that the de"ning equations as proposed impose that baseline
glycemia be smaller than or at most equal to pancreatic target
glycemia, otherwise no equilibrium solution can be obtained. The
third result (Proposition I.6) poses a qualitative limit to the behavior of
the solutions in the case that b is exactly equal to G , that is, in
 @
this case there cannot be any solutions where glycemia progressively
decreases towards the target equilibrium value. If it is considered
that, at steady state, basal levels of insulin are observed and that
therefore pancreatic activity continues (i.e. the true pancreatic set-point
seems in reality to be lower than baseline glycemia), then the physio-
logic likelihood of the model seems to be somewhat diminished by
these results.
If, on one hand, the original minimal model was actually designed
to cope just with the three hours after the administration of the glucose
150 A. De Gaetano, O. Arino

bolus, making in"nite-time asymptotic behavior somewhat less impor-

tant, on the other hand the de"nition itself of the Insulin Sensitivity
index S is strictly valid only at in"nite time. Moreover, the steady state
'
conditions which the present aggregated model has di$culties with are
not only those established many hours after the administration of the
glucose bolus, but also those prevailing a few minutes before such
administration.
In addition, if the current pancreas model is retained, its
nonautonomous form, requiring a de"nite starting time, makes its
application di$cult to any experimental set-up di!erent from a theor-
etically instantaneous IV bolus glucose injection. For all of these
reasons, it seems that the system that we could obtain by piecing
together the available components of the minimal model should be
changed.
It must be noted that the minimal model incorporates some basic
ideas that should be kept under consideration in any further work.
First of all, the action of the pancreas on peripheral tissue utilization of
glucose is not immediate. In the minimal model this delay has been
represented as the progressive accumulation of the active hormone in
an intermediate compartment (X), and while other formalizations may
be preferred, the idea itself is important and should be retained.
Secondarily, a delay is also introduced, via the nonautonomous term
t in the third equation, on the actual insulin incretion. That a delay is
present is apparent from the classically biphasic shape of the insulin
concentration curve after a glucose bolus: while "rst-phase insulin
response occurs immediately (indicating the availability of readily
released hormone), second-phase insulin response appears over several
tens of minutes, indicating either the slow release of the hormone from
previously stored reserves (di!erent from those responsible for "rst-
phase insulin release) or the necessity of de-novo synthesis. The form
chosen to represent this delay is the reason of many of the analytic
di$culties of the aggregated model, and even from a physiologic point
of view this simpli"cation, requiring the pancreas to secrete second-
phase insulin with linearly increasing response to the glucose stimulus
as time progresses, can be improved. However, while it is easy to
improve on the description of the pancreatic response to glucose, it is
di$cult to do so without introducing additional parameters. In fact,
the most obvious alternative would be to model pancreatic response as
a nonlinear curve with respect to time, attaining asymptotically a max-
imum of responsiveness for long times. This curve needs, unfortunately,
two parameters to be described (the maximum responsiveness to the
unit glucose stimulus and the time to half-maximum responsiveness),
instead of the single slope parameter in the minimal model. Since the
Mathematical modelling of IVGTT 151

experiment on which the parameters are estimated consists only of

about twenty points for glucose and twenty for insulin, estimating nine
parameters instead of eight is a major experimental cost which we
would like to avoid, if possible. In fact, if at all possible, we would like
to reduce the number of parameters needed to model the experiment to
a smaller number than eight.
One possible form of a dynamical model for the glucose-
insulin system is that described by Equations (4) and (5). It depends on
six free parameters overall and may exhibit a secondary insulin peak.
The model admits only one equilibrium point, which represents the
resting, basal glucose and insulin values for the subject. The stability
of the model around this equilibrium point is assured, and num-
erical evaluation of its parameters presents no problems, with
small coe$cients of variation and good overall "tting to the
experimental observations. Solutions to this model are positive and
bounded.
While no claims are made that the dynamical model presented may
be more than a working approximation to a structure which diabetol-
ogists and physiologists might consider acceptable, its very existence
underscores two basic concepts which we believe ought to be taken
into account by any future improved models. First of all, the whole
self-regulatory system of insulin and glucose (at least) should be
simultaneously considered and can actually be modelled to yield
a reasonably close approximation to actual data. Second, the dynamic
properties of any suggested model ought to be formally investigated in
order to subject to closer scrutiny and meaningful numerical estima-
tion those models which are known structurally to possess desirable
overall characteristics, like stability, or positivity of solutions, or boun-
dedness of solutions and the like.

Appendix: Proofs

I. Formal study of the minimal model

We would like to establish some basic results "rst:

Proposition I.0.
i) lim inf I(t)7I ;
@
R
ii) lim inf X(t)70;
R
iii) lim sup G (t) 6G .
@
R
152 A. De Gaetano, O. Arino

If, furthermore, I(0)'I , then

@
iv) I (t)'I ∀ t70;
@
v) X (t)70 ∀ t70;
vi) if G (¹ )(G for some ¹'0, then G(t)(G ∀ t'¹.
@ @
Proof. These results are easily obtained by standard arguments on
di!erential inequalities, and the proof is skipped.

Proposition I.1. Suppose G 'p , lim sup G(t)'p . ¹hen

@  R 
lim sup X(t)"R.
R
Proof. If XPR there is nothing to show. Let us suppose then that
lim inf X(t) "nite.
R
Solve Eq. (3) to obtain


R
I(t)!I "p e\NR#p e\N R\S [G(u)!p ]>udu;
@   

solving Eq. (2) and substituting we get


R
X(t)" e\N R\S p [I (u)!I ] du,
 @

or


R
X (t)"p p e\NR eN\NS du
 



 
R S
#p p e\N R\S e\N S\Q [G(s)!p ]>sds du.
  
 
Now, we have two possible cases, p "p , or p 9p .
   
Introduce the auxiliary function z, de"ned as


t, if p "p

R  
z(t)" eN\NS du" eN\NR!1
 p !p
, if p 9 p
 
 
and, changing the order of integration in the double integral at right
hand-side above we obtain


R
X (t)"p p e\NR z (t)#p p e\NR\Q z (t!s) [G(s)!p ]> sds.
    

Mathematical modelling of IVGTT 153

Now, since X (t)#p '0 for all su$ciently large t, the derivative of

G is bounded superiorly, and under the assumptions made,  +t ,L1,
L
t PR, h'0, ¹'0 U ∀ n 3 -, G (s)7p #h for t 6s6t #¹.
L  L L
For t 6s6t #¹, both the "rst term and the factors of the
L L
integrand above are positive quantities. Now, z (t)71, and we get


t #¹
L
X (t #¹)7p p t h e\N RL>2\Q ds,
L   L t
L



X (t #¹)7p p t h eNQ ds
L  L
\2
and
X (t #¹ ) PR, as nPR.
L
Proposition I.2. Suppose lim sup G(t)(p ; then G 6p .
R  @ 
Proof. Suppose G 'p . By Eq. (3),
@ 
dI (t)
lim sup G(t)(p N "!p [I (t)!I ] N I (t)PI
 dt  @ @
R
NX(t) P0NG (t)PG : contradiction.
@
Proposition I.4. For any value p (G the system does not admit an
 @
equilibrium.
Proof. Letting dG/dt"dX/dt"dI/dt"0 we can easily see that the
only possible equilibrium state is (G , 0, I ). We claim that for p (G
@ @  @
this is not an equilibrium state, hence that there is no equilibrium. If
there were an equilibrium at (G , 0, I ), Eq. (3) would give
@ @
dI/dt"p (G !p )'0, a contradiction.
 @ 
Proposition I.5. If the subject is considered to be at steady state before
the glucose bolus, then G must be lesser than or equal to p .
@ 
Proof. Suppose we consider the subject's glucose-insulin system in the
steady state condition preceding the glucose bolus. Then all derivatives
in Equations (1), (2) and (3) are zero, t"0\, i.e. can be considered
arbitrarily small preceding the bolus, and X"0 from the initial condi-
tion of Eq. (2). It follows from Eq. (2) that I(0\)"I and from Eq. (1)
@
that G(0\)"G , as expected. Therefore Eq. (3) would imply that
@
G 6p .
@ 
Proposition I.6. ¸et p '0, p "G , G(0)'G . Assume I(0)'I . ¹hen
  @ @ @
 ¹'0 U G(t)'G ∀t(¹, G (¹ )"G , and G (t)(G ∀ t'¹.
@ @ @
154 A. De Gaetano, O. Arino

Proof. We will "rst establish the fact that G (t)PG as tPR. This is
@
clear if G (¹ )(G for some ¹, since if G (t)(G after some time,
@ @
eventually I (t)PI and X (t)P0, so that G (t)PG . However, assum-
@ @
ing G (t)'G ∀ t'0, since from proposition (I.0) lim sup G(t)6
@ R
G , then lim G(t)"G .
@ R @
Suppose now that G (t)'G ∀ t'0: we want to show that we
@
reach a contradiction, in which case the theorem is proved.
Let us study again the above integral. Since R X(s) G(s) ds is
bounded (following from Eq. (1) by formal integration  and using
 R X(s) G(s) ds(G (0)!G (t)6(G(0)!G ) e\NR ), then
 @

 
R R
G X (s) ds6 X (s) G(s) ds(#R,
@
 
which means that


R
X (s) ds(#R.

Integrating Eq. (2) we obtain

 
R R
X(t)"X(t)!X(0)"!p X(s) ds#p (I (s)!I ) dsN
  @
 

 
R R
X (t)#p X(s) ds"p (I(s)!I ) ds,
  @
 
but I (t)7I , and therefore the integral on the right is positive. Now,
@
lim inf X(t)"0, otherwise X would integrate to in"nity, so
R
 +t ,P#R U X(t )P0 and
L L

 
RL RL
X(t )#p X(s) ds"p (I(s)!I ) ds,
L   @
 
and as t PR we have
L

 
 
0#p X (s) ds"p (I (s)!I ) ds,
  @
 
so that the integral on the right is "nite since that on the left is "nite.
Moreover, letting tPR in the equation

 
R R
X(t)#p X(s) ds"p
(I(s)!I )ds,
  @
 
we have that X (t )P0 (not only its lim inf does).
Mathematical modelling of IVGTT 155

We notice in passing that the same result is obtained if p


is replaced with a function that tends to zero as t goes to in"nity.
In the present case, since p is independent of t, G (t)PG at least
 @
exponentially:

dG
6!p (G(t)!G ) N G(t)!G 6e\NR (G(0)!G ),
dt  @ @ @

so that with G (t )'G ∀ t, G(t)PG at least exponentially. This means

@ @
that p (G (t)-G )>t"p (G(t)-G ) t also goes to zero exponentially, so
 @  @
that from Eq. (3) we may write

dI
6!p [I(t)!I ]#d(t), with 06d(t)6C e\IR, 0(k(p .
dt  @ 

This di!erential inequality of course implies that [I (t)!I ]P0 at

@
least exponentially as tPR.
With the same reasoning, letting X be majored by a function which
decreases exponentially we deduce that X (t)P0 exponentially as
tPR.
Now, we will try to estimate more closely the rates of convergence
of the three functions G, X and I, to evidence a contradiction. For
the following of the present proof it is clearer to switch to the nota-
tion u (t)"G (t)!G , and v (t)"I (t)!I . Notice that in the context
@ @
of the present proof u (t)'0 ∀ t; we rewrite therefore the model equa-
tions as

du
"!p u!Xu!XG
dt  @

dX
"!p X#p v
dt  

dv
"!p v#tu
dt 

We integrate the three equations between t and R (knowing that they

decrease to zero at least exponentially and are therefore integrable).
Integrating, remembering that u (R)"X (R)"v (R)"0, and
changing all signs we obtain:

  
  
u(t)"p u(s) ds# X(s) u(s) ds#G X (s) ds
 @
R R R
156 A. De Gaetano, O. Arino

 
 
X(t)"p X (s) ds!p v(s) ds
 
R R

 
 
v (t)"p v(s) ds# su (s) ds

R R
We deduce some relations: from the second model equation we have

 
 
p X(s) ds"X(t)#p v (s) ds,
 
R R
and since all terms on the right-hand side are positive,

 
 p 
X(s) ds 7  v (s) ds.
p
R  R
Similarly, from the third equation we get

 
 
p v (s) ds7 su (s) ds,

R R
but since s7t, su (s)7t u (s), so

 
 t 
v (s) ds7 u(s) ds,
p
R  R
and incorporating the inequality above,

 
 p t 
X(s) ds 7  u(s) ds.
p p
R   R
From the "rst integrated model equation we get

 
 
u(t)7p u (s) ds#G X(s) ds,
 @
R R
and replacing for the integral of X


p t 
u (t)7 p #  G u (s) ds.
 p p @
  R
Now let


 p G
w (t)" u(s) ds, a"  @ .
p p
R  
Mathematical modelling of IVGTT 157

We may rewrite the last inequality as

dw w(t)
6!(p #at) w(t) or 6!(p #at),
dt  w(t) 

and integrating we obtain



w (t) R w(t) t
log 6! [p #as] ds or log 6! p t#a ,
w(0)  w(0)  2

that is

,
t
! p t#a
w (t)6w (0) e  2

which means that w (t)P0 faster than exponentially, or

 
 a
!p t! t

R a
!p t! t

u(s) ds6Ce 2 N u(s) ds6Ce 2 .

R R
But u is a decreasing function (for su$ciently large t), hence


R
tu (2t)6 u(s) ds,
R
and therefore
1 !p t!a t
u (2t)6 Ce  2 .
t
Calling r"2t, h"a/2, and taking t'1 we can write u (r)6C e\FP,

so that we can see u decreases faster than exponentially. Therefore,
since



u(t)7G X(s) ds,
@
R
we have that


p 
u (t)7G  v(s) ds,
@p
 R
so that



v(s) dsP0
R
faster than exponential, i.e.,



v(s)ds6C e\FR.

R
158 A. De Gaetano, O. Arino

On the other hand,


dv  v(0)
7!p v (t)Nv(t)7e\NR v(0)N v(s) ds7 e\NR.
dt  b
R 
Therefore



C e\NR6 v(s) ds6C e\FR :
 
R
contradiction.
Therefore it is not possible that G (t) remains above G for all times,
@
and at a "nite time the solution passes below G .
@
For completeness we also consider the case when G (p . We
@ 
notice that (G"G , X"0, I"I ) is an equilibrium point. If we start
@ @
from any initial data (G(0), 0, I (0)), we see that

dI
7!p (I (t)!I )Nlim inf (I(t))7I Nlim inf X (t)70;
dt  @ @
R> R>
dG
6![p !e] G (t )#p G ; t7¹Nlimsup G(t)6G .
dt   @ @
R>
We conclude that there exists a ¹ such that

dI
t5¹NG(t)6G #e(p N "!p (I (t )!I ), t7¹
@  dt  @

which implies that I(t)PI , X(t)P0, G(t)PG . In this case, therefore,

@ @
the behavior of the model is as expected. Further, we remark that we
may have both G (t)'G ∀ t, and G (t)(G for some time, depending
@ @
on the relative magnitude of the coe$cients, p 'min (p , p ) with
  
p '0 being a su$cient condition for G(¹)"G for some ¹. Notice
 @
that a secondary insulin peak can only occur when G (t)'p .

It is interesting to evaluate whether there exist time limits to the
occurrence of the secondary insulin peak. Indeed, the peak may only
occur when the second derivative of insulin is negative and its "rst
derivative is zero. Considering that G (t)'p 'G , we may write
 @
dI dI dG
"!p #p [G(t)!G ]#p t
dt  dt  @  dt

"p [1!p t] [G(t)!G ]!p t X(t) G(t)

  @ 
"!t [p p (G(t)!G )#p X(t)G(t)]#p [G(t)!G ]
  @   @
Mathematical modelling of IVGTT 159

So, requiring the second derivative to be negative, we have that the

time to the secondary peak, call it t , must satisfy the relation

[G(t )!G ] 1
t '  @ 7 .
 p [G(t )!G ]#X(t ) G(t ) p #X (t )
  @    
From the second model equation, using the di!erential inequality with
initial data 0, we have that
p
X (t ) 6 max  [I(s)!I ],
 06t6t p
 
@

and since G is exponentially decreasing, we can only have a "nite

number of peaks. As far as the second one is concerned,
p
X (t )6  max (I(0)!I , I(t )!I ),
 p @  @

and
1
t 7 .
 p
p #  max (I(0)!I , I(t )!I )
 p @  @

But from the third model equation we get, for 06t6t ,

dI
!p [I(t)!I ]#p [G (0)!G ] t ,
dt  @  @ 

since G (t) is decreasing, and therefore

 
p [G(0)!G ] t
I (t)!I 6 max I(0)!I ,  @  ,
@ @ p

so that
1
t 7 ,



p p [G(0)!G ] t
p #  max I(0)!I ,  @ 
 p @ p
 
which can be rewritten as

    
p \ p p [G(0)!G ] t \
t 7min p #  (I(0)!I ) , p #   @  .
  p @  p p
  
Solving explicitly for t in the second lower bound, we obtain
p p
  [G(0)!G ] t#p t !170.
p p @  
 
160 A. De Gaetano, O. Arino

Of the two roots, one is negative and is not applicable, hence

!p p p #(p p p#4p p p p [G(0)!G ]

t 7           @ ,
 2p p [G(0)!G ]
  @
and we can "nally write that the time t to the second peak must satisfy




p \
t 7min p #  (I(0)!I ) ,
  p @



!p p p #(p p p#4p p p p [G (0)!G ]
          @ .
2p p [G(0)!G ]
  @

II. Stability of the dynamical model

For notational convenience we rewrite the model as:

(M1)
dG
II.1) "!b G(t)!b I(t)G(t)#b ,
dt   

G(t),G ∀t 3 [!b ,0), G(0)"G #b ,

@  @ 


dI b 
II.2) "!b I (t)#  G(t#s) ds, I(0)"I #b b .
dt  b !b @  

The initial conditions translate physiologically into the requirement
that the subject be at his/her (presumably equilibrium) state (G , I ) for
@ @
at least b minutes prior to (and excluding) the bolus injection of the

glucose load. From this equilibrium state, the concentrations of glu-
cose and insulin are supposed to jump instantaneously to new values
determined by the amount of glucose administered and by the corre-
sponding "rst-phase pancreatic secretion of insulin.
Since we are dealing with a delayed di!erential equation, the state
of the system is properly de"ned by the present concentration of insulin
I (t), together with all present and past concentrations of glucose from
time (t!b ) to the present time t. The past concentrations of glucose

may be regarded as given by a piecewise continuous function c from
the speci"ed time interval [!b , 0] into the set of admissible glucose

values (the positive reals). We will employ the usual notation in the
domain of delay-di!erential equations [13], writing G (s)"G (t#s)
R
with !b 6s60.

Mathematical modelling of IVGTT 161

At equilibrium,
*"(G*, I*)2, 0"!b G*!b G* I*#b ,
  
0"!b I*#b G*.
 
We had assumed that (since the subject is assumed at equilibrium
before perturbation)
I
b "b G #b G I , b "b @ ,
  @  @ @  G
@
therefore


I
b @ G*#b G*!(b G #b G I )"0
G   @  @ @
@
G*
I*"I
@G
@


G*"G
and since we require G*70, then @.
I*"I
@
There is then only one equilibrium point with positive concentra-
tions, *"(G , I ).
@ @
In order to show stability of the system around the equilibrium
point, we can follow one of two approaches: linearize the system
around the equilibrium point and study the real parts of the eigen-
values of the resulting linear system [13,15], or show convergence to
the equilibrium for every solution by majoring and minoring at in"n-
ity. While the "rst approach is the more usual one, it leads, in our case,
to a condition for stability in terms of the parameters on which it is not
easy to improve. We will instead demonstrate stability for any set of
parameter values following the second approach.
We need to establish beforehand boundedness of the solutions.
Proposition II.3. ¹he solutions +G(t), I(t), of M1 are positive and
bounded.
Proof. By inspection of Eq. (II.2), since G , I and all b's are positive,
@ @
I (t) is positive for all times. Hence, by Eq. (II.1), G(t) is positive; further,
it cannot exceed max (G #b , b /b ). Hence G is bounded. Substitu-
@   
ting the bound for G in the integral in Eq. (II.2), we can see that I (t)
cannot exceed max (I #b b , b b /(b b )), and is bounded.
@      
Proposition II.4. ¹he time derivatives of the solutions are bounded.
Proof. Obvious from Proposition II.3 and the form of the derivatives
in Equations (4) and (5).
A key idea is the following:
162 A. De Gaetano, O. Arino

Proposition II.5. Any sequence of points in phase space along a solution

of M1 admits a convergent subsequence.

Proof. May show precompactness coordinate-wise.

For I (t) it is enough to consider that the values of the sequence are
reals in a bounded set (by Proposition II.3), close the set with +0, and
the bound, and use the Heine-Borel property of 1.
For G : the values of the sequence are functions in C and Proposi-
R
tions II.3 and II.4 imply that the sequence is bounded as a subset of
C (since both the norm of the function and the norm of its derivative
are bounded). Therefore the conditions of the Arzela-Ascoli theorem
hold [24] and the sequence is precompact.

Corollary II.6. ¹he u-limit set of the system M1 is not empty.

De5nition II.7. Pick a solution F (t)"+G(t), I (t), of M1, and from it

a converging subsequence. Denote as F "+g, u, the point of conver-

gence. ¸et F* (t)"+G* (t), I* (t), be the solution of M1 starting from the
initial condition F "+g, u,.

We remark that since F is in the u-limit set, which is positively

invariant [25], the solution F* (t) is de"ned for every positive real t and
F* can be estended as a solution for negative real t so that F* is in the

u-limit set for every t. Also, we will refer to the result on the approxi-
mation of the solutions to the limit solution (Saperstone IV,4) (25) as:

Lemma II.8. ¹here exists a sequence +¹ ,PR U G(t#¹ )PG* (t),

L L
I (t#¹ )PI* (t) as ¹ PR.
L L
Now we have assembled all the necessary machinery and we
proceed to prove stability.

Proposition II.9. lim sup G(t)6b /b . Moreover, if for some t 70

R   
G (t )6b /b , then G (t)(b /b ∀t't .
     
Proof. Majoring the solution for G we have


dG(t) R
6!b G (t)#b N G (t)6e\@R G(0)#b e\@R\Q ds
dt   

b
"e\@R G(0)#  (1!e\@R),
b

so that
b
lim sup G(t) 6  .
R b

Mathematical modelling of IVGTT 163

Similarly, applying this inequality for t't , the second statement follows.

We note in particular that if G (0) 6 b /b , G (t)(b /b ∀ t'0.
   
b
Corollary II.10. G* (t)6  ∀t.
b

Proof. Immediate by Proposition II.9 and Lemma II.8.
Having thus obtained a majoring bound for G*, for every couple
(G*, I* ) in the u-limit set, we use it to obtain a majoring bound for I*:
b b
Let x"   .
b b
 
Proposition II.11. I*(t)6x. Moreover, if for some t 70/I (t )6x and
 
b
G (t )6 , then I(t)(x ∀t't .
 b 

Proof.



dI* b R b b
"!b I*(t)#  G* (s) ds6!b I* (t)#  b  ∀ t 3 1.
dt  b !b  b  b
  
Therefore, for arbitrary t 6t, knowing that I* is bounded,



R b b
I* (t)6e\@ R\R I* (t )# e\@R\Q   ds,
 b
t
 
or
I* (t)6e\@R\R I*(t )#x (1!e\@R\R).

Now let t P!R to obtain the "rst result. For the second statement,

use the maxima given in the inequality
I (t)6e\@R\R I(t )#x (1!e\@R\R)

obtained in the same way.
A majoring bound on I* for each (G*, I*) in the u-limit set allows
us to derive a minoring bound on G*:
b
Proposition II.12. G* (t)7  . Moreover, if for some t 70 and
b #b x 
 
for all t't ,

b
G (t )7  and I(t)6x,
 b #b x
 
b
then G(t)'  ∀ t't .
b #b x 
 
164 A. De Gaetano, O. Arino

Proof.

dG*
"!b G* (t)!b G* (t) I* (t)#b
dt   

7!b G* (t)!b x G* (t)#b ∀t31.

  
Knowing that G* is bounded, and by the same reasoning as for
Proposition II.11, the "rst statement follows. The second statement
follows similarly to II.10 and II.11.

b b
Having minored G* we can now minor I*. Let p"   .
bb #b b b
    
Proposition II.13. I*(t)7p . Moreover, if for some t 70

I (t )7b p and G (t )7b p
   
then
I (t)'b p ∀ t't .
 
Proof.


dI* b R
"!b I*(t)#  G* (s) ds
dt  b t!b

dI* b
N 7!b I* (t)# b b p
dt  b  

and as for II.11 and for II.12, the stated inequality follows. The second
statement follows in the same fashion as II.11.
Having a lower bound for I*, we can "nd an upper bound for G*.
b
Let c"  .
b #b b p
  
Proposition II.14. G*(t)6c. Moreover, if for some t 60 and

∀ t't G(t ) 6c and I (t)'b p, then G (t)(c ∀ t't .
   
Proof. dG*/dt"!b G* (t)!b I* (t ) G* (t)#b 6!(b #b b p)
     
G* (t)#b , and as for II.11, the "rst statement follows. The statement

follows in the same fashion as II.12.

At this point, we can see the regularity of the general process: if

G* is bounded above by some a, then I* is bounded above by a/b ,

Mathematical modelling of IVGTT 165

then G* is bounded below by b /(b #b a/b ), then I* is bounded

   
below by b /(b (b #b a/b )), then G* is bounded above by
    
b /(b #b b /(b (b #b a/b ))), and so on. Similarly, if G is bounded
        
above by some a, and I (0) is bounded.
We can build two sequences, m and M , made up of the successive
H
lower and upper bounds to G*.
Proposition II.15. ¸et +m ,, +M , be the sequences respectively of
H H
lower and upper bounds to G* (t) obtained by repeated application
of Propositions II.11 through II.14, with M "b /b , m "
   
b /(b #b b b /(b b )). ¹hen
      
b
i) M "h (M ), m "h (m ), with h (x)"  ;
H> H H> H b b
b #  
 b b #b x
  
ii) h is monotonically increasing in x;
iii) h (M )6h (M ), h (m )7h (m );
H> H H> H
iv) +M ,, +m , are bounded.
H H
Proof.
i) by direct computation;
ii) obvious by inspection;
iii) obvious because of (ii) and since M (M and m 'm by direct
   
computation;
iv) since h monotonically increasing and M 'm , then M 'm ∀ j
  H
by induction. From this fact and from iii, all terms in the two
sequences +M , and +m , are bounded by M above and by
H H 
m below.

Now we have two monotonic sequences on R which are bounded,
hence converge. A point of convergence for any one of the two se-
quences must be a "xed point of h, and we study therefore the number
of "xed points of h.
Proposition II.16. h has only one positive ,xed point.
Proof. Let x be a "xed point of h. Then x"h (x) or
b
x"  N b b x#bb x!b b b "0 .
b b       
b #  
 b b #b x
  
Keeping in mind that all coe$cients are positive, the discriminant is
positive hence all roots are real; one root is negative, the other positive.
Hence the proposition.
166 A. De Gaetano, O. Arino

Since h only has one positive "xed point, and since both +M , and
H
+m , converge to a "xed point of h, then there must be only one value
H
to which both +M , and +m , converge, hence they converge to the
H H
same value. Since G* is bounded by the two sequences, it follows that
G* is constant, and all possible points in the u-limit set share the same
constant G component. Inserting this one value for G* (t) in the
R
original system, we obtain only one possible value for I*. We conclude
that the u-limit set consists exactly of one point. Since the equilibrium
point (G , I ) belongs to the u-limit set, it must be its only point.
@ @
Therefore, any solution to the original system converges to (G , I ).
@ @
Now, in order to demonstrate stability we need that solutions for
G and I remain within respective arbitrary intervals [G , G ], [I ,
  
I ], when starting from within suitable intervals contained in them.

But for each couple of intervals [G , G ], [I , I ], we can "nd
   
a couple of subintervals

 
b b
[m , M ],  m ,  M
H H b H b H
 
for some j such that
b b
m 'G , M (G ,  m 'I ,  M (I ,
 H  b H  b H 
 
by repeated applications of statements II.11 through II.14, knowing
that the two sequences +M , and +m , converge to the same point.
H b H b
Once G (t ) 3 [m , M ] and I(t ) 3 [ b m , b M ] , the solution remains
 H H  H H
within those intervals for every t't , again because of Propositions

II.11}II.14. Hence the system is stable and, because of global conver-
gence, asymptotically stable.

Acknowledgements. Thanks go to Prof. A.V. Greco and to Dr. G. Mingrone (Division

of Metabolic Diseases, Università Cattolica, Rome) for having kindly provided the
experimental data. The authors wish to thank Prof. Karl Hadeler and an anonymous
referee of the Journal of Mathematical Biology for the very thorough and valuable
review of the "rst draft of the manuscript.

References

1. Ackerman, E., L. C. Gatewood, J. W. Rosevear, and G. D. Molnar. Model studies

of blood-glucose regulation. Bull. Math. Biophys. 27: Suppl:21}Suppl:37, 1965
2. Bergman, R. N., Y. Z. Ider, C. R. Bowden, and C. Cobelli. Quantitative estimation
of insulin sensitivity. Am. J. Physiol. 236: E667}E677, 1979
3. Bolie, V. W. Coe$cients of normal blood glucose regulation. J. Appl. Physiol. 16:
783}788, 1961
Mathematical modelling of IVGTT 167

4. Buchanan, T. A., B. E. Metzger, N. Freinkel, and R. N. Bergman. Insulin sensitivity

and B-cell responsiveness to glucose during late pregnancy in lean and moderately
obese women with normal glucose tolerance or mild gestational diabetes. Am. J.
Obstet. Gynecol. 162: 1008}1014, 1990
5. Castagneto, M., A. De Gaetano, G. Mingrone, R. M. Tacchino, G. Nanni,
E. Capristo, G. Benedetti, P. A. Tataranni, and A. V. Greco. Normalization of
insulin sensitivity in the obese patient after stable weight reduction with bilio-
pancreatic diversion. Obesity Surgery 4: 161}168, 1994
6. Ceresa, F., F. Ghemi, P. F. Martini, P. Martino, G. Segre, and A. Vitelli. Control of
blood glucose in normal and in diabetic subjects. Studies by compartmental
analysis and digital computer technics. Diabetes 17: 570}578, 1968
7. De Gaetano, A., G. Mingrone, M. Castagneto, P. A. Tataranni, and A. V. Greco.
NONMEM provides improved parameter estimation for the minimal model of
glucose kinetics. Am. J. Physiol. 271 (Endocr.Metab.34): E932}E937, 1996
8. Defronzo, R. A. Insulin resistance, hyperinsulinemia, and coronary artery disease:
a complex metabolic web. J. Cardiovasc. Pharmacol. 20 Suppl 11: S1}16, 1992
9. Defronzo, R. A. and E. Ferrannini. Insulin resistance. A multifaceted syndrome
responsible for NIDDM, obesity, hypertension, dyslipidemia, and atherosclerotic
cardiovascular disease. Diabetes Care 14: 173}194, 1991
10. Defronzo, R. A., J. D. Tobin, and R. Andres. Glucose clamp technique: a method for
quantifying insulin secretion and resistance. Am. J. Physiol. 237: E214}E223, 1979
11. Frontoni, S., L. Ohman, J. R. Haywood, R. A. Defronzo, and L. Rossetti. In vivo
insulin action in genetic models of hypertension. Am. J. Physiol. 262: E191}E196,
1992
12. Gatewood, L. C., E. Ackerman, J. W. Rosevear, G. D. Molnar, and T. W. Burns.
Tests of a mathematical model of the blood-glucose regulatory system. Comput.
Biomed. Res. 2: 1}14, 1968
13. Hale, J. K., S. M. Verduyn Lunel. Introduction to Functional Di!erential
Equations, Springer Verlag, New York, 1993
14. Jansson, L., L. Lindskog, and N. E. Norden. Diagnostic value of the oral glucose
tolerance test evaluated with a mathematical model. Comput. Biomed. Res. 13:
512}521, 1980
15. MacDonald, N. Biological delay systems: linear stability theory. Cambridge:
Cambridge University Press, 1989,
16. Nomura, M., M. Shichiri, R. Kawamori, Y. Yamasaki, N. Iwama, and H. Abe.
A mathematical insulin-secretion model and its validation in isolated rat pancre-
atic islets perifusion. Comput. Biomed. Res. 17: 570}579, 1984
17. O'Connor, M. D., H. Landahl, and G. M. Grodsky. Comparison of storage- and
signal-limited models of pancreatic insulin secretion. Am. J. Physiol. 238:
R378}R389, 1980
18. Osei, K., D. A. Cottrell, and M. M. Orabella. Insulin sensitivity, glucose e!ec-
tiveness, and body fat distribution pattern in nondiabetic o!spring of patients with
NIDDM. Diabetes Care 14: 890}896, 1991
19. Pacini, G. and R. N. Bergman. MINMOD: a computer program to calculate insulin
sensitivity and pancreatic responsivity from the frequently sampled intravenous
glucose tolerance test. Comput. Methods Programs Biomed. 23: 113}122, 1986
20. Page, R., M. Boolell, A. Kalfas, S. Sawyer, R. Pestell, G. Ward, and F. Alford.
Insulin secretion, insulin sensitivity and glucose-mediated glucose disposal in
Cushing's disease: a minimal model analysis. Clin. Endocrinol. Oxf. 35: 509}517,
1991
 168 A. De Gaetano, O. Arino

21. Piccardo, M. G., G. Pacini, M. Rosa, and R. Vichi. Insulin resistance in myotonic
dystrophy. Enzyme 45: 14}22, 1991
22. Press, W. H., B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling. Numerical
recipes in C. The art of scienti"c computing. New York: Cambridge University
Press, 1988,
23. Proietto, J. Estimation of glucose kinetics following an oral glucose load. Methods
and applications. Horm. Metab. Res. Suppl. 24: 25}30, 1990
24. Rudin, W. Principles of mathematical analysis. New York: McGraw-Hill, 1976
25. Saperstone, S. H. Semidynamical systems in in"nite dimensional spaces. New
York: Springer Verlag, 1981
26. Sluiter, W. J., D. W. Erkelens, P. Terpstra, W. D. Reitsma, and H. Doorenbos.
Glucose tolerance and insulin release, a mathematical approach. II. Approxima-
tion of the peripheral insulin resistance after oral glucose loading. Diabetes 25:
245}249, 1976
27. Subba Rao, G., J. S. Bajaj, and J. Subba Rao. A mathematical model for
insulin kinetics. II. Extension of the model to include response to oral glucose
administration and application to insulin-dependent diabetes mellitus (IDDM).
J. Theor. Biol. 142: 473}483, 1990
28. To!olo, G., R. N. Bergman, D. T. Finegood, C. R. Bowden, and C. Cobelli.
Quantitative estimation of beta cell sensitivity to glucose in the intact organism:
a minimal model of insulin kinetics in the dog. Diabetes 29: 979}990, 1980
29. Turner, R. C., A. S. Rudenski, D. R. Matthews, J. C. Levy, S. P. O'Rahilly, and
J. P. Hosker. Application of structural model of glucose-insulin relations to assess
beta-cell function and insulin sensitivity. Horm. Metab. Res. Suppl. 24: 66}71, 1990
30. Yang, Y. J., J. H. Youn, and R. N. Bergman. Modi"ed protocols improve insulin
sensitivity estimation using the minimal model. Am. J. Physiol. 253: E595}E602,
1987

View publication stats