Fall 2010Modeling DiabetesMath 636Diabetes 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 aredramatically 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, oftencalled juvenile diabetes, and Type 2, often referred to as adult onset diabetes (which now occursin children as young as 5). For our studies here, we will concentrate on Type 1 diabetes, whichis an autoimmune disease, and represents only 10% of all cases of diabetes. Type 1 diabetesis a hereditary disease, which occurs in about 4-20 per 100,000 people with peak occurrencearound 14 years of age. The classical symptoms of diabetes are increased hunger (polyphagia), increased thirst (polydipsia), and frequent urination (polyuria) . These symptoms tend to develop rapidly (weeksor 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 dehydrationfrom frequent and dilute urine causes the symptom of being thirsty. Other symptoms includeblurred vision caused by glucose absorption in the lenses of the eyes, fatigue, weight loss, andpoor wound healing. Type 1 diabetes may result is diabetic ketoacidosis, where the urine smellsof acetone. Diabetes increases the risk of heart disease, especially because of atherosclerosis from lowinsulin. Over time, diabetes can damage blood vessels and nerves, which increases the chanceof foot injury and decreases the body’s ability to ﬁght infection. These problems can becomesu±ciently severe as to require amputations. The osmotic imbalances over time result is kidneydamage (nephropathy) and nerve damage (neuropathy). The increased pressure in the eye oron the optic nerve can lead to blindness (retinopathy). Thus, diabetes when untreated has veryserious consequences. Glucose MetabolismWe begin with a short discussion of glucose metabolism. We ingest food to obtain energy tosustain our bodies. The carbohydrates in the food are broken down into simple sugars, whichare absorbed into the blood. This raises the concentration of glucose in the blood, where cellscan access it for metabolism into energy. However, when glucose levels get too high, then pressures (osmotic?) increase that cancause problems in the tissues. For normal subjects, the rise in glucose concentration causesthe Ø-cells in the pancreas to release insulin into the blood (along with a number of otherhormones). Insulin aÆects glucose concentration in several ways, including the facilitation ofglucose transport across cell membranes, especially in skeletal muscles, and conversion of glucoseto 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. Thisnegative 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 theglycogen and produce glucose. This hormone works opposite insulin to increase blood glucose. The glucocorticoids help metabolize carbohydrates, especially in the liver, and also help increaseblood glucose concentrations. Growth hormone can block the eÆects of insulin by reducingthe liver uptake of glucose and decrease muscle sensitivity to insulin. There are many otherhormones that regulate glucose levels in the blood, creating a complex regulatory system thatis 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 incurssome unknown environmental assault that initiates the auto-immune system to attack theirown Ø-cells. When the Ø-cells are destroyed, they boost the immune system to further attackmore Ø-cells, leaving the body without the ability to produce insulin. This severely limits itsability 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 TestOur modeling eÆorts begin with a simple model developed by Ackerman et al [1,2] in the1960s, which is based on the Glucose Tolerance Test (GTT). After a 12 hour fast, subject isgiven a large amount of glucose (1. 75 mg of glucose/kg of body weight). This sugar is rapidlyingested, then blood sugar is monitored for the next 4-6 hours. The glucose concentration inthe 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, wewant to develop a simple model with a few parameters that can be ﬁt to the data generatedby the GTT. For simplicity, we create a system of diÆerential equations that follow the bloodconcentrations of glucose (G(t)) and insulin (I(t)) though the later can be thought of as thecomplex soup of hormones in the body regulating blood glucose. The general model is written:dGdtdIdt= f1 (G, I) + J(t),= f2 (G, I),where J(t) is the ingested source of glucose. (This is an external control function representedglucose 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 inputtedto the system. More complex models that have been developed and are continually beingimproved 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 andinsulin. We assume that the body wants to maintain a homeostasis for glucose concentrationsin the blood. Also, we are assuming that J(t) is acting like a ±-function, so only aÆects theinitial conditions and is eÆectively zero away from t = 0. The homeostasis assumption meansthat we want to consider a local perturbation of the dynamical system away from equilibrium. Thus, create the perturbation variables,g(t) = G(t) ° G0andi(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 linearizedperturbation model given by:dgdtdidt==@f1 (G0 , I0 )@f1 (G0 , I0 )g+i,dgdi@f2 (G0 , I0 )@f2 (G0 , I0 )g+i,dgdiwhere 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 andf2 , with our understanding of the physiology of glucose and insulin. An increase in glucose inthe blood stimulates tissue uptake of glucose and glycogen storage in the liver. Also, increasesin insulin facilitate the uptake of glucose in tissues and the liver. Hence, it is clear that@f1 (G0 , I0 )= °m1 <, 0dg@f1 (G0 , I0 )= °m2 <, 0. diandHowever, increases in blood glucose result in the release of insulin, while increases in insulinonly result increased metabolism of excess insulin. These physiological facts imply that@f2 (G0 , I0 )= m4 >, 0dg@f2 (G0 , I0 )= °m3 <, 0. diandWith these deﬁnitions, the linearized system can be written:µ ∂g˙˙i=µ°m1m4°m2°m3∂µ ∂gi,where g = dg/dt and similarly for i(t). ˙The characteristic equation for this linear system is given byØØ °m ° ∏°m2Ø1det ØØm4°m3 ° ∏ØØØØ = ∏2 + (m1 + m2 )∏ + m1 m3 + m2 m4 = 0. ØBy our deﬁnitions above, this characteristic equation only has positive coe±cients. From basicdiÆerential equations (think spring mass system), this implies that the solutions ∏ are eithercomplex with negative real part or both eigenvalues are negative reals. Both situations give astable 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 linearizedsolution 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 belowequilibrium) that encourages more eating. It follows that the general solution is given byg(t) = e°Æt (c1 cos(!t) + c2 sin(!t)),whereÆ=m1 + m32and!=If we takec1 = A cos(!±)1q4(m1 m3 + m2 m4 ) ° (m1 + m3 )2 . 2andc2 = A sin(!±),then we can approximate the blood glucose level byG(t) = G0 + Ae°Æt cos(!(t ° ±)). (1)This solution has ﬁve unknown parameters to be ﬁt to the data. The parameter G0 representsthe equilibrium blood sugar level, Æ measures the ability of the system to return to equilibriumstate after being perturbed, and ! gives a frequency response to perturbations. One mightexpect that measuring Æ should be the primary measure of whether someone was diabetic, aspeople 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 Ackermanet al [1,2]. A more robust measure was the natural frequency of the system, !0 . (Recall thenatural frequency from forced damped oscillators in elementary diÆerential equations. ) Wedeﬁne2º2!0 = ! 2 + Æ2andT0 =,!0where T0 is the natural period of the system. The natural period turns out to be a goodpredictor of diabetes from Ackerman et al’s experiments [1,2]. In particular, they found that ifT0 <, 4, then a person was generally normal, while if T0 >, 4, then the person is likely to havediabetes. Physiologically, you might relate to this by the idea that normally we get hungryevery 3-4 hours. ExamplesWe examine the theory described above with a normal and a diabetic subject given theGTT. The table below gives the data collected on two subjects. t (hr)00. 50. 7511. 522. 5346Subject A70150165145907565758075Subject B1001852102201951751051008590Table 1: Data from the Glucose Tolerance Test. Subject A is a normal subject, while SubjectB has diabetes. The data in the table are ﬁt to the model Eqn. (1). A least squares best ﬁt is performedusing either Excel or MatLab. Below is a table of the best ﬁtting parameters for each of thesubjects, including the least sum of square errors. The MatLab programs and Excel spreadsheetare put together in a GTT. zip ﬁle for anyone interested. ParameterG0ÆA!±LSSESubject A79. 18140. 9927171. 54671. 81270. 90056225. 6757Subject B95. 21250. 6335263. 15211. 03041. 51604718. 6180Table 2: Best Fitting Parameters to GTT Model. Subject A is a normal subject, while SubjectB has diabetes. The models and the data are graphed in the ﬁgure below. One can readily see that the bestparameter ﬁt does very well matching the model to the data. From the deﬁnitions of !0 andT0 above, we ﬁnd that for Subject A,!0 = 2. 0668andT0 = 3. 0401,so according to the criterion by Ackerman et al [1,2], this subject is clearly normal. ForSubject B, we ﬁnd!0 = 1. 2095andT0 = 5. 1947,so according to the criterion by Ackerman et al [1,2], this subject is clearly diabetic. GTT Model250Glucose (mg/dl)2001501005000123t (hr)456Diabetes in NOD MiceSince diabetes is such a signiﬁcant human disease, an animal model has been developed togain a better understanding in the scientiﬁc community. One important animal with a diabetictendency is the non-obese diabetic (NOD) mouse. Type 1 diabetes arises in NOD mice whenT cells from the immune system become primed to speciﬁcally target and kill beta-cells. Thesecytotoxic T cells belong to a class of lymphocytes displaying a surface marker called CD8(denoted CD8+ T cells). T Cell Immune ResponseWe start with a brief discussion of the T cell immune response, which is central to ourmodel for the onset of diabetes. The T cells mature in the thymus, and most T cells thatcross-react with self-proteins are destroyed to minimize autoimmune responses. Next the Tcells 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, themajor histocompatibility complex or MHC. Depending on the strength, duration, and numberof interactions that a particular T cell has determines whether it undergoes activation andissues an immune response. When a T cell with the appropriate speciﬁcity matches a foreign protein, then the activatedT cell proliferates. This activated T cell can rapidly reproduce, undergoing approximately 6cell divisions, to create eÆector cells (also called cytotoxic T-lymphocytes or CTL’s), whichseek out and destroy target cells, which protects the host from this foreign invader. TheseCTL’s can be dangerous in the body, so are short-lived. Alternately, the activated T cell canissue a weaker response (mostly when there is not as much of the foreign protein present) andproduce long-lived memory cells, which have the same speciﬁcity, but wait until the stimulus isencountered again in larger quantities to mount an immune response that destroys the target. Autoimmunity for Type 1 DiabetesAn autoimmune disease, such as Type 1 diabetes, is caused by the immune cells attackingthe wrong targets. In particular, the T cells attack the Ø-cells in the pancreas that produceinsulin. 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 isreduced in these mice, which triggers an autoimmune response. LYMPH NODEPANCREASApoptotic β cellnaive T cellp−MHCpeptide pactivationβ CellapoptosisDendritic cellf1injuryAEEffector (CTL) T cellsActivated T cellf2MMemory cellsRecent 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.  track the levels of auto-reactive CD8+ Tcells 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) inthe 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 individualanimals provide a predictor for the animal becoming diabetic. The data for each one of themice were aligned to the time of onset of high-blood sugar symptoms. The the pooled datashow 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 amodel that could shed some light on this process. Model for Diabetes in NOD MiceOur 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), ormemory 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 tothe production of the speciﬁc antigen peptide, p(t), which feedback and aÆects the number ofactivated T cells. A schematic for this model is presented below. f1f1Mf2A1−f2BEpThis mathematical model assumes that on the time-scale of interest, a system of ordinarydiÆerential equations for a well-mixed compartment model can be applied. More details onspeciﬁc assumptions are given in MahaÆy and Keshet . The detailed full model satisﬁes theﬁrst order system of diÆerential equations given by:dAdtdMdtdEdtdpdtdBdt= (æ + Æ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 functionsf1 (p) =f2 (p) =pn,nk1 + pnak2 m. mk2 + pmThe nonlinear feedback function, f1 (p), represents the activation of the T cells by p(t) andis given the saturation/switching form noted above that includes the Hill coe±cient n andthe kinetic constant, k1 . The other nonlinear feedback function, f2 (p), is a negative feedbackfunction that represents how many activated T cells become memory cells. This function hasthe classic Michaelis-Menten (inhibition) form with a Hill coe±cient of m. Below we show theform of these functions. memoryactivationf 2 (p)f1 (p)k2k1pThe diÆerential equation for activated T cells in (2) begins with the production of activated Tcells from naive T cells (a constant æ) and memory cells (ÆM ). The activated T cells are lost bybecoming 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 lineardecay (rate ±M ). A similar equation describes the dynamics of the eÆector cells with productionand ampliﬁcation (Ø2m2 (1 ° f2 (p))) from activated T cells and linear decay (rate ±E ). The lasttwo 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 modelassumes that the Ø-cells are not replaced and that the peptide is rapidly removed (rate ±p ). Model ReductionThe full model is a highly nonlinear model consisting of ﬁve diÆerential equations and 17parameters. Experimental evidence limits the choice of parameters and gives insight into waysof how this system can be reduced. Biological constraints simplify the analysis to a signiﬁcantlysmaller region of the parameter space. However, a ﬁve dimensional system is too complex fordetailed analysis. It is important to note that the last two of the diÆerential equations of System (2) operateon diÆerent time scales than the ﬁrst three diÆerential equations. The loss of Ø-cells is a longterm process lasting weeks, so the dynamics of B can be considered more like a slow movingparameter. 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 thatit remains eÆectively in equilibrium relative to the other variables. This allows the diÆerentialequation to be expressed as an algebraic equation. From the equation for p in System (2), theQSS yields:dpRB= 0,sopºE. dt±pIf B is considered a slow moving parameter, then the dynamics reduces to the 3-dimensionalsystem:dA= (æ + ÆM )f1 (p) ° (Ø + ±A )A ° ≤A2 = F1 (A, M, E),dtdM= Ø2m1 f2 (p)A ° f1 (p)ÆM ° ±M M = F2 (A, M, E),(3)dtdE= Ø2m2 (1 ° f2 (p))A ° ±E E = F3 (A, E),dtwhere p =RB±p E. Equilibrium and Linear AnalysisThe 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 equilibriumis 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 ofbiologically relevant parameters there are two additional equilibria. These can be readily foundnumerically, but have no analytic solution. A linear analysis of this model is performed in the usual manner, ﬁnding the Jacobian matrixand evaluating it near the equilibria. 0 @F1 (A,M,E)@ABJ(A, M, E) = @ @F2 (A,M,E)@A@F3 (A,E)@A@F1 (A,M,E)@M@F2 (A,M,E)@M0@F1 (A,M,E)@E@F2 (A,M,E)@E@F3 (A,E)@E1CAProvided n >, 1 (which is expected, since f1 (p) is a type of switch), then the linearization aboutthe 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 astable node. One would expect the disease-free state to be stable as any small perturbation ofthe equilibrium representing a minor assault on the body should result in the body returningto 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 thatcorresponds to a state with all immune cell levels remaining elevated. In that state, eÆectorT cells are continuously killing Ø-cells, and this corresponds to an autoimmune attack, whicheventually results in diabetes. This equilibrium has various stability properties that depend onthe parameters and merit further detailed discussion below. The third equilibrium is a saddlewith a two-dimensional stable manifold, which for some parameters separates the “healthy” anddiseased equilibria. For these parameters, stimuli that fall on the wrong side of this separatrixwill be attracted to the diseased equilibrium. For other parameter values, the unstable manifoldof the diseased state connects to the stable manifold of the saddle point. In this case, almostall 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 allthree equilibria with the saddle node very close to the M -axis. One trajectory from the saddlenode (red) along one branch of its unstable manifold connects the saddle node equilibrium withthe disease-free equilibrium at the origin. The other trajectory (green) emanating in the otherdirection along the saddle node’s unstable manifold wanders until it eventually approaches alimit cycle about the diseased equilibrium. The graph of Fig. (2) on the right shows a close upDiabetes Model0. 110. 080. 80. 060. 04E0. 60. 020. 430. 200. 2*0. 152. 52010. 80. 60. 10. 41. 50. 80. 0510. 60. 20. 50. 40. 20000AMFigure 2: 3D Phase portraits for reduced model. phase portrait near the diseased equilibrium (asterisk) with the solution trajectory approachinga stable limit cycle. Bifurcation AnalysisOne 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 thisparameter appears in the QSS approximation, p = RB E. Note also that we ignored the±pdiÆ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 . The3D reduced model (3) was analyzed, and bifurcation diagrams were composed with the AUTOfeature of XPP, freely available software written by G Bard Ermentrout1 . Y1Y123. 51. 831. 62. 51. 41. 2211. 50. 80. 610. 40. 50. 20000. 511. 52a152. 533. 5424(a)681012141618a15(b)Figure 3: Bifurcation diagram for the peptide decay rate, a15 = ±p . The vertical axis is A inunits of 103 cells. (a) A portion of the diagram, enlarged, shows the typical bifurcation: AHopf bifurcation occurs at a15 = 0. 5707 spawning a stable limit cycle. A homoclinic bifurcationoccurs 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 (anduses the default parameters values given in Tables ?? and ??. Moving across this diagram from1XPP is freely available at www. math. pitt. edu/~bard/xpp/xpp. htmlleft to right along the horizontal axis represents increasing values of the peptide decay rate ±p ,or equivalently, a decreasing level of beta cells, B. . . .
Originally posted 2018-07-07 18:53:17. Republished by Blog Post Promoter