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

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º

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

pº

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….

**30 %**discount on an order above

**$ 100**

Use the following coupon code:

RESEARCH