# Fall 2010 Modeling Diabetes Math 636 Diabetes is a disease that is characterized

August 30, 2017

Question
Fall 2010

Modeling Diabetes

Math 636

Diabetes is a disease that is characterized by excessive glucose in the blood stream. Currently, there is an epidemic of diabetes that has resulted from unhealthy lifestyles, which are
dramatically diÆerent from how humans survived when they evolved from small nomadic huntergatherer societies, and food was di±cult to ﬁnd. There are two forms of diabetes, Type 1, often
called juvenile diabetes, and Type 2, often referred to as adult onset diabetes (which now occurs
in children as young as 5). For our studies here, we will concentrate on Type 1 diabetes, which
is an autoimmune disease, and represents only 10% of all cases of diabetes. Type 1 diabetes
is a hereditary disease, which occurs in about 4-20 per 100,000 people with peak occurrence
around 14 years of age.
The classical symptoms of diabetes are increased hunger (polyphagia), increased thirst (polydipsia), and frequent urination (polyuria) [3]. These symptoms tend to develop rapidly (weeks
or months) in Type 1 diabetes. The frequent urination results from the poor water reabsorption in the kidney because of the osmotic imbalance of glucose in the blood. The dehydration
from frequent and dilute urine causes the symptom of being thirsty. Other symptoms include
blurred vision caused by glucose absorption in the lenses of the eyes, fatigue, weight loss, and
poor wound healing. Type 1 diabetes may result is diabetic ketoacidosis, where the urine smells
of acetone.
Diabetes increases the risk of heart disease, especially because of atherosclerosis from low
insulin. Over time, diabetes can damage blood vessels and nerves, which increases the chance
of foot injury and decreases the body’s ability to ﬁght infection. These problems can become
su±ciently severe as to require amputations. The osmotic imbalances over time result is kidney
damage (nephropathy) and nerve damage (neuropathy). The increased pressure in the eye or
on the optic nerve can lead to blindness (retinopathy). Thus, diabetes when untreated has very
serious consequences.
Glucose Metabolism
We begin with a short discussion of glucose metabolism. We ingest food to obtain energy to
sustain our bodies. The carbohydrates in the food are broken down into simple sugars, which
are absorbed into the blood. This raises the concentration of glucose in the blood, where cells
can access it for metabolism into energy.
However, when glucose levels get too high, then pressures (osmotic?) increase that can
cause problems in the tissues. For normal subjects, the rise in glucose concentration causes
the Ø-cells in the pancreas to release insulin into the blood (along with a number of other
hormones). Insulin aÆects glucose concentration in several ways, including the facilitation of
glucose transport across cell membranes, especially in skeletal muscles, and conversion of glucose
to glycogen in the liver, which provides a good storage of glucose for future consumption. Thus,
increasing the concentration of insulin results in blood glucose concentration decreasing. This
negative feedback system helps the body tightly regulate glucose levels to maintain a balance.
There are other signiﬁcant hormones that are involved in the regulation of blood glucose.
In response to high energy demands, epinephrine (adrenalin) is released to break down the
glycogen and produce glucose. This hormone works opposite insulin to increase blood glucose.
The glucocorticoids help metabolize carbohydrates, especially in the liver, and also help increase
blood glucose concentrations. Growth hormone can block the eÆects of insulin by reducing
the liver uptake of glucose and decrease muscle sensitivity to insulin. There are many other
hormones that regulate glucose levels in the blood, creating a complex regulatory system that
is crucial to maintenance of blood glucose for energy to all cells in the body.

Type 1 diabetes occurs when someone who is genetically predisposed to the disease incurs
some unknown environmental assault that initiates the auto-immune system to attack their
own Ø-cells. When the Ø-cells are destroyed, they boost the immune system to further attack
more Ø-cells, leaving the body without the ability to produce insulin. This severely limits its
ability to regulate glucose and results in the onset of diabetes. Because of the immune response,
the body cannot regenerate new Ø-cells nor can transplants succeed.
Glucose Tolerance Test
Our modeling eÆorts begin with a simple model developed by Ackerman et al [1,2] in the
1960s, which is based on the Glucose Tolerance Test (GTT). After a 12 hour fast, subject is
given a large amount of glucose (1.75 mg of glucose/kg of body weight). This sugar is rapidly
ingested, then blood sugar is monitored for the next 4-6 hours. The glucose concentration in
the blood of the subject is measured and these data are ﬁt to a model.
As noted above, the glucose regulatory system in the body is very complex. However, we
want to develop a simple model with a few parameters that can be ﬁt to the data generated
by the GTT. For simplicity, we create a system of diÆerential equations that follow the blood
concentrations of glucose (G(t)) and insulin (I(t)) though the later can be thought of as the
complex soup of hormones in the body regulating blood glucose. The general model is written:
dG
dt
dI
dt

= f1 (G, I) + J(t),
= f2 (G, I),

where J(t) is the ingested source of glucose. (This is an external control function represented
glucose sources coming from food ingested.) For the GTT, we can think of J(t) as a ±-function,
since we give an initial large quantity of glucose after a fast, then have no other glucose inputted
to the system. More complex models that have been developed and are continually being
improved often examine ideal ways to match J(t) with diÆerent foods.
The above model is very general, so gives little insight into the dynamics of glucose and
insulin. We assume that the body wants to maintain a homeostasis for glucose concentrations
in the blood. Also, we are assuming that J(t) is acting like a ±-function, so only aÆects the
initial conditions and is eÆectively zero away from t = 0. The homeostasis assumption means
that we want to consider a local perturbation of the dynamical system away from equilibrium.
Thus, create the perturbation variables,
g(t) = G(t) ° G0

and

i(t) = I(t) ° I0 ,

where G0 and I0 are the equilibrium values for blood glucose and insulin concentrations, respectively. Thus,
f1 (G0 , I0 ) = f2 (G0 , I0 ) = 0.
We expand the general model to linear terms with these deﬁnitions, yielding the linearized
perturbation model given by:
dg
dt
di
dt

=
=

@f1 (G0 , I0 )
@f1 (G0 , I0 )
g+
i,
dg
di
@f2 (G0 , I0 )
@f2 (G0 , I0 )
g+
i,
dg
di

where g(t) and i(t) now represent the linearized perturbed variables.

For the next step in our analysis, we examine the partial derivatives of the functions, f1 and
f2 , with our understanding of the physiology of glucose and insulin. An increase in glucose in
the blood stimulates tissue uptake of glucose and glycogen storage in the liver. Also, increases
in insulin facilitate the uptake of glucose in tissues and the liver. Hence, it is clear that
@f1 (G0 , I0 )
= °m1 < 0
dg

@f1 (G0 , I0 )
= °m2 < 0.
di

and

However, increases in blood glucose result in the release of insulin, while increases in insulin
only result increased metabolism of excess insulin. These physiological facts imply that
@f2 (G0 , I0 )
= m4 > 0
dg

@f2 (G0 , I0 )
= °m3 < 0.
di

and

With these deﬁnitions, the linearized system can be written:
µ ∂

g
˙
˙
i

=

µ

°m1
m4

°m2
°m3

∂µ ∂

g
i

,

where g = dg/dt and similarly for i(t).
˙
The characteristic equation for this linear system is given by
Ø
Ø °m ° ∏
°m2
Ø
1
det Ø
Ø
m4
°m3 ° ∏

Ø
Ø
Ø
Ø = ∏2 + (m1 + m2 )∏ + m1 m3 + m2 m4 = 0.
Ø

By our deﬁnitions above, this characteristic equation only has positive coe±cients. From basic
diÆerential equations (think spring mass system), this implies that the solutions ∏ are either
complex with negative real part or both eigenvalues are negative reals. Both situations give a
stable equilibrium as we would expected from this self-regulatory system.
We are only measuring the blood glucose level in the GTT, so we only need the linearized
solution for g(t). We expect the underdamped situation with complex eigenvalues. Physiologically, you can think of the your body’s response to a “sugar high” (maximum of blood glucose),
which is followed after an hour or two by a “sugar low” (minimum of blood glucose below
equilibrium) that encourages more eating. It follows that the general solution is given by
g(t) = e°Æt (c1 cos(!t) + c2 sin(!t)),
where
Æ=

m1 + m3
2

and

!=

If we take
c1 = A cos(!±)

1q
4(m1 m3 + m2 m4 ) ° (m1 + m3 )2 .
2
and

c2 = A sin(!±),

then we can approximate the blood glucose level by
G(t) = G0 + Ae°Æt cos(!(t ° ±)).

(1)

This solution has ﬁve unknown parameters to be ﬁt to the data. The parameter G0 represents
the equilibrium blood sugar level, Æ measures the ability of the system to return to equilibrium
state after being perturbed, and ! gives a frequency response to perturbations. One might
expect that measuring Æ should be the primary measure of whether someone was diabetic, as
people with diabetes should not be able to return rapidly to normal equilibrium levels. However,
the parameter Æ was found to have large errors from the many subjects tested by Ackerman

et al [1,2]. A more robust measure was the natural frequency of the system, !0 . (Recall the
natural frequency from forced damped oscillators in elementary diÆerential equations.) We
deﬁne

2
!0 = ! 2 + Æ2
and
T0 =
,
!0
where T0 is the natural period of the system. The natural period turns out to be a good
predictor of diabetes from Ackerman et al’s experiments [1,2]. In particular, they found that if
T0 < 4, then a person was generally normal, while if T0 > 4, then the person is likely to have
diabetes. Physiologically, you might relate to this by the idea that normally we get hungry
every 3-4 hours.
Examples
We examine the theory described above with a normal and a diabetic subject given the
GTT. The table below gives the data collected on two subjects.
t (hr)
0
0.5
0.75
1
1.5
2
2.5
3
4
6

Subject A
70
150
165
145
90
75
65
75
80
75

Subject B
100
185
210
220
195
175
105
100
85
90

Table 1: Data from the Glucose Tolerance Test. Subject A is a normal subject, while Subject
B has diabetes.

The data in the table are ﬁt to the model Eqn. (1). A least squares best ﬁt is performed
using either Excel or MatLab. Below is a table of the best ﬁtting parameters for each of the
subjects, including the least sum of square errors. The MatLab programs and Excel spreadsheet
are put together in a GTT.zip ﬁle for anyone interested.
Parameter
G0
Æ
A
!
±
LSSE

Subject A
79.1814
0.9927
171.5467
1.8127
0.90056
225.6757

Subject B
95.2125
0.6335
263.1521
1.0304
1.51604
718.6180

Table 2: Best Fitting Parameters to GTT Model. Subject A is a normal subject, while Subject
B has diabetes.

The models and the data are graphed in the ﬁgure below. One can readily see that the best
parameter ﬁt does very well matching the model to the data. From the deﬁnitions of !0 and
T0 above, we ﬁnd that for Subject A,
!0 = 2.0668

and

T0 = 3.0401,

so according to the criterion by Ackerman et al [1,2], this subject is clearly normal. For
Subject B, we ﬁnd
!0 = 1.2095
and
T0 = 5.1947,
so according to the criterion by Ackerman et al [1,2], this subject is clearly diabetic.
GTT Model
250

Glucose (mg/dl)

200

150

100

50

0
0

1

2

3
t (hr)

4

5

6

Diabetes in NOD Mice
Since diabetes is such a signiﬁcant human disease, an animal model has been developed to
gain a better understanding in the scientiﬁc community. One important animal with a diabetic
tendency is the non-obese diabetic (NOD) mouse. Type 1 diabetes arises in NOD mice when
T cells from the immune system become primed to speciﬁcally target and kill beta-cells. These
cytotoxic T cells belong to a class of lymphocytes displaying a surface marker called CD8
(denoted CD8+ T cells).
T Cell Immune Response
We start with a brief discussion of the T cell immune response, which is central to our
model for the onset of diabetes. The T cells mature in the thymus, and most T cells that
cross-react with self-proteins are destroyed to minimize autoimmune responses. Next the T
cells migrate to the lymph nodes, where they interact with antigen presenting cells (APC’s)
that display small fragments of proteins (about 9 amino acids) primarily from foreign sources,
such as viruses or bacteria. These antigens are held inside a cleft of a larger protein, the
major histocompatibility complex or MHC. Depending on the strength, duration, and number
of interactions that a particular T cell has determines whether it undergoes activation and
issues an immune response.
When a T cell with the appropriate speciﬁcity matches a foreign protein, then the activated
T cell proliferates. This activated T cell can rapidly reproduce, undergoing approximately 6

cell divisions, to create eÆector cells (also called cytotoxic T-lymphocytes or CTL’s), which
seek out and destroy target cells, which protects the host from this foreign invader. These
CTL’s can be dangerous in the body, so are short-lived. Alternately, the activated T cell can
issue a weaker response (mostly when there is not as much of the foreign protein present) and
produce long-lived memory cells, which have the same speciﬁcity, but wait until the stimulus is
encountered again in larger quantities to mount an immune response that destroys the target.
Autoimmunity for Type 1 Diabetes
An autoimmune disease, such as Type 1 diabetes, is caused by the immune cells attacking
the wrong targets. In particular, the T cells attack the Ø-cells in the pancreas that produce
insulin. Early in the development of NOD mice there is a wave of programmed cell death
(apoptosis) of pancreatic Ø-cells. It is conjectured that the clearance of the apoptotic cells is
reduced in these mice, which triggers an autoimmune response.
LYMPH NODE

PANCREAS
Apoptotic β cell

naive T cell

p−MHC

peptide p
activation

β Cell
apoptosis

Dendritic cell

f1

injury
A
E

Effector (CTL) T cells

Activated T cell

f2

M

Memory cells

Recent experiments suggest that a speciﬁc peptide fragment of the protein IGRP (glucose6-phosphatase catalytic subunit-related protein) may be the primary culprate in NOD mice.
Experiments have been designed to track levels of speciﬁc T cells circulating in the blood.
Speciﬁcally, some experiments of Trudeau et al.[5] track the levels of auto-reactive CD8+ T
cells in NOD mice for several weeks. Surprisingly, the levels of the CTL’s show dramatic
ﬂuctuations over time as seen in Fig 1.

Figure 1: Periodic waves of circulating T-cells occur in mice prone to diabetes (NOD mice) in
the weeks before the onset of the disease. Dark line, circles: T-cell level. Grey line, squares:
percentage of the animals that became diabetic.

Not all NOD mice develop diabetes, but the presence of cyclic T cell waves for individual
animals provide a predictor for the animal becoming diabetic. The data for each one of the
mice were aligned to the time of onset of high-blood sugar symptoms. The the pooled data
show three peaks in the level of T cells. The amplitude of the oscillations increases with time,
and we observe a slight increase in the period of oscillation. Below we present and analyze a
model that could shed some light on this process.
Model for Diabetes in NOD Mice
Our full model examines ﬁve state variables. The ﬁrst variable is the activated T cells,
denoted A(t). These T cells can become killer T cells (CTL’s) or eÆector cells, denoted E(t), or
memory cells, denoted M (t). The eÆector cells speciﬁcally attack the Ø-cells in the pancreas,
so we track the fraction of Ø-cells that remain, B(t). The destruction of the Ø-cells leads to
the production of the speciﬁc antigen peptide, p(t), which feedback and aÆects the number of
activated T cells. A schematic for this model is presented below.
f1
f1

M
f

2

A
1−f

2

B

E

p

This mathematical model assumes that on the time-scale of interest, a system of ordinary
diÆerential equations for a well-mixed compartment model can be applied. More details on
speciﬁc assumptions are given in MahaÆy and Keshet [4]. The detailed full model satisﬁes the
ﬁrst order system of diÆerential equations given by:
dA
dt
dM
dt
dE
dt
dp
dt
dB
dt

= (æ + ÆM )f1 (p) ° (Ø + ±A )A ° ≤A2 ,
= Ø2m1 f2 (p)A ° f1 (p)ÆM ° ±M M,
= Ø2m2 (1 ° f2 (p))A ° ±E E,

(2)

= REB ° ±p p,
= °∑EB,

with nonlinear feedback functions
f1 (p) =
f2 (p) =

pn
,
n
k1 + pn
ak2 m
.
m
k2 + pm

The nonlinear feedback function, f1 (p), represents the activation of the T cells by p(t) and
is given the saturation/switching form noted above that includes the Hill coe±cient n and
the kinetic constant, k1 . The other nonlinear feedback function, f2 (p), is a negative feedback
function that represents how many activated T cells become memory cells. This function has

the classic Michaelis-Menten (inhibition) form with a Hill coe±cient of m. Below we show the
form of these functions.
memory

activation

f 2 (p)

f1 (p)

k2

k1

p

The diÆerential equation for activated T cells in (2) begins with the production of activated T
cells from naive T cells (a constant æ) and memory cells (ÆM ). The activated T cells are lost by
becoming eÆector and memory T cells (rate Ø), decaying ( rate ±A ), and competing with others
(≤A2 ). The equation for the memory cells includes the production and ampliﬁcation (Ø2m1 f2 (p))
from activated T cells and the loss due to producing more activated T cells (f1 (p)ÆM ) and linear
decay (rate ±M ). A similar equation describes the dynamics of the eÆector cells with production
and ampliﬁcation (Ø2m2 (1 ° f2 (p))) from activated T cells and linear decay (rate ±E ). The last
two diÆerential equations describe the dynamics of how the eÆector T cells, E, destroy Ø-cells,
B, producing the protein, p, that activates T cells via the feedback function, f1 (p). This model
assumes that the Ø-cells are not replaced and that the peptide is rapidly removed (rate ±p ).
Model Reduction
The full model is a highly nonlinear model consisting of ﬁve diÆerential equations and 17
parameters. Experimental evidence limits the choice of parameters and gives insight into ways
of how this system can be reduced. Biological constraints simplify the analysis to a signiﬁcantly
smaller region of the parameter space. However, a ﬁve dimensional system is too complex for
detailed analysis.
It is important to note that the last two of the diÆerential equations of System (2) operate
on diÆerent time scales than the ﬁrst three diÆerential equations. The loss of Ø-cells is a long
term process lasting weeks, so the dynamics of B can be considered more like a slow moving
parameter. At the other end of the time scale is the kinetics for the peptide p, which is very rapid
(on the scale of hours). This latter information allows a quasi-steady state (QSS) assumption.
The QSS assumption on a variable is that the dynamics of this variable are su±ciently fast that
it remains eÆectively in equilibrium relative to the other variables. This allows the diÆerential
equation to be expressed as an algebraic equation. From the equation for p in System (2), the
QSS yields:
dp
RB
= 0,
so

E.
dt
±p
If B is considered a slow moving parameter, then the dynamics reduces to the 3-dimensional
system:
dA
= (æ + ÆM )f1 (p) ° (Ø + ±A )A ° ≤A2 = F1 (A, M, E),
dt
dM
= Ø2m1 f2 (p)A ° f1 (p)ÆM ° ±M M = F2 (A, M, E),
(3)
dt
dE
= Ø2m2 (1 ° f2 (p))A ° ±E E = F3 (A, E),
dt

where p =

RB
±p E.

Equilibrium and Linear Analysis
The equilibria for (3) are found by solving the three highly nonlinear equations:
F1 (Ae , Me , Ee ) = 0,

F2 (Ae , Me , Ee ) = 0,

F3 (Ae , Ee ) = 0.

From the form of the equations and the functions f1 and f2 , it is easy to see that one equilibrium
is the disease-free or trivial equilibrium:
(Ae , Me , Ee ) = (0, 0, 0).
From the nature of the nonlinearities, there may be from zero to four other equilibria. However,
there are relatively stringent biological constraints on the parameters, and in the range of
biologically relevant parameters there are two additional equilibria. These can be readily found
numerically, but have no analytic solution.
A linear analysis of this model is performed in the usual manner, ﬁnding the Jacobian matrix
and evaluating it near the equilibria.
0 @F1 (A,M,E)
@A
B
J(A, M, E) = @ @F2 (A,M,E)
@A
@F3 (A,E)
@A

@F1 (A,M,E)
@M
@F2 (A,M,E)
@M

0

@F1 (A,M,E)
@E
@F2 (A,M,E)
@E
@F3 (A,E)
@E

1
C
A

Provided n > 1 (which is expected, since f1 (p) is a type of switch), then the linearization about
the disease-free state gives the characteristic equation:
(∏ + Ø + ±A )(∏ + ±M )(∏ + ±E ) = 0.
Thus, the eigenvalues for the disease-free state are all negative, which makes this equilibrium a
stable node. One would expect the disease-free state to be stable as any small perturbation of
the equilibrium representing a minor assault on the body should result in the body returning
to the disease-free state.
The two non-zero equilibria result in a much more complicated characteristic equation.
However, solutions are readily found numerically. There exists one positive equilibrium that
corresponds to a state with all immune cell levels remaining elevated. In that state, eÆector
T cells are continuously killing Ø-cells, and this corresponds to an autoimmune attack, which
eventually results in diabetes. This equilibrium has various stability properties that depend on
the parameters and merit further detailed discussion below. The third equilibrium is a saddle
with a two-dimensional stable manifold, which for some parameters separates the “healthy” and
diseased equilibria. For these parameters, stimuli that fall on the wrong side of this separatrix
will be attracted to the diseased equilibrium. For other parameter values, the unstable manifold
of the diseased state connects to the stable manifold of the saddle point. In this case, almost
all positive initial conditions asymptotically, approach the “healthy” state.
Before continuing the analysis of this 3D model, we examine some solution trajectories of (3)
with parameters in the range predicted for NOD mice. The graph of Fig. (2) on the left shows all
three equilibria with the saddle node very close to the M -axis. One trajectory from the saddle
node (red) along one branch of its unstable manifold connects the saddle node equilibrium with
the disease-free equilibrium at the origin. The other trajectory (green) emanating in the other
direction along the saddle node’s unstable manifold wanders until it eventually approaches a
limit cycle about the diseased equilibrium. The graph of Fig. (2) on the right shows a close up

Diabetes Model

0.1
1

0.08
0.8

0.06
0.04

E

0.6

0.02
0.4

3

0.2

0
0.2

*

0.15

2.5
2

0
1

0.8
0.6

0.1

0.4

1.5
0.8

0.05

1

0.6

0.2

0.5

0.4
0.2
0

0

0

0

A

M

Figure 2: 3D Phase portraits for reduced model.
phase portrait near the diseased equilibrium (asterisk) with the solution trajectory approaching
a stable limit cycle.
Bifurcation Analysis
One conjecture for the onset of diabetes is that peptide fragments from the apoptosis of
Ø-cells in the pancreas are not cleared su±ciently fast, which allows an autoimmune response.
It follows that the bifurcation analysis should center on the parameter ±p . Recall that this
parameter appears in the QSS approximation, p = RB E. Note also that we ignored the
±p
diÆerential equation for the loss of Ø-cells, B, which also appears in this QSS approximation.
In fact, the loss of Ø-cells is equivalent to increasing the peptide clearance parameter, ±p . The
3D reduced model (3) was analyzed, and bifurcation diagrams were composed with the AUTO
feature of XPP, freely available software written by G Bard Ermentrout1 .
Y1

Y1

2

3.5

1.8
3
1.6
2.5

1.4
1.2

2

1
1.5

0.8
0.6

1

0.4
0.5
0.2
0

0
0

0.5

1

1.5

2
a15

2.5

3

3.5

4

2

4

(a)

6

8

10

12

14

16

18

a15

(b)

Figure 3: Bifurcation diagram for the peptide decay rate, a15 = ±p . The vertical axis is A in
units of 103 cells. (a) A portion of the diagram, enlarged, shows the typical bifurcation: A
Hopf bifurcation occurs at a15 = 0.5707 spawning a stable limit cycle. A homoclinic bifurcation
occurs at a15 = 2.268. (b) Further bifurcations on an expanded scale: another Hopf bifurcation
(to an unstable limit cycle) occurs at a15 = 4.063. This limit cycle vanishes at a15 = 20.28.
The diagram given in Figure 3(a) shows the basic bifurcation behaviour of the model (and
uses the default parameters values given in Tables ?? and ??. Moving across this diagram from
1

XPP is freely available at www.math.pitt.edu/~bard/xpp/xpp.html

left to right along the horizontal axis represents increasing values of the peptide decay rate ±p ,
or equivalently, a decreasing level of beta cells, B….

Get a 30 % discount on an order above \$ 100
Use the following coupon code:
RESEARCH
Positive SSL