September 25, 2017, Monday
University of Colorado at Boulder Search A to Z Campus Map University of Colorado at Boulder CU 
Search Links

MBW:Extensions to a 2-Delay Glucose-Insulin Regulatory Model

From MathBio

Jump to: navigation, search

The following wiki is a reproduction of the results produced by Jiaxu Li, Yang Kuang and Clinton Mason in 2006 [1]. Additionally, we seek to extend the original analysis in order to examine problems associated with dieting, meal type and diabetes. For a complete review of the original paper by Li et al. please visit MBW: Glucose-Insulin Regulatory System.

Executive summary

In this wiki, we replicate and extend a model for the ultradian glucose-insulin oscillations observed in the human body for a variety of health conditions. We first reproduce the results from a paper by Li, Kuang and Mason (2006), and then extend this model to account for various meal types and glucose-insulin regulatory system behavior under both type I and type II diabetes. The human glucose-insulin regulatory system is an especially important to model, since the CDC estimates that the number of cases of diabetes in the United States has increased from 12.0 million in 2000 to 25.8 million in 2011. In 2000, diabetes affected 4.4% of the population and now impacts 8.3% of the population. On an international scale, diabetes will significantly impact the health of many people. In India alone diabetes will kill over 1 million people next year and presents a significant danger to individuals escaping poverty[2]. A better understanding of how medications and diet affect glucose and insulin levels can lead to better treatment decisions, helping prevent serious blood-glucose-related conditions such as hypoglycemia, neuropathy, blindness and even death, from occurring.

We employ a variety of mathematical techniques in extending this model and analyzing results. The model is based on two delay differential equations (DDEs) that describe glucose and insulin levels in the body. There are two delays built into the model, one that describes hepatic glucose production and one that models \beta -cell insulin secretion. The differential system is solved using MATLAB’s numerical solver DDE23. As part of our model analysis, we examine the location and significance of any bifurcation points in the model’s behavior. We propose some corrections to the paper by Li et al. and we examine how changes in the initial history of the delay differential equations affect model behavior. We find that the model responds well to extensions that simulate various dietary and diabetic conditions and that some basic statistical analysis provides justification for current trends in research.

Biological Context

While a complete biological context for the mathematical model is given in the review wiki found at MBW: Glucose-Insulin Regulatory System, we briefly review the biological framing of the glucose-insulin regulatory system and type I and II diabetes. In order to function properly the human body aims to maintain a stable level of plasma-glucose concentration that enables cells to function appropriately. The glucose-insulin regulatory system is designed to promote appropriate levels of glucose and responds to both increases in decreases in glucose concentration. As the regulatory system responds there are two types of oscillations that are produced, neither of which is biologically well understood. The first type of oscillations are rapid and occur on a time scale ranging from 5 - 15 minutes, capturing the body's acute response to large changes in glucose and insulin concentration. The second type of oscillations, ultradian oscillations, occur on a time scale ranging from 50 - 150 minutes. These ultradian oscillations capture the body's large-scale response to glucose and insulin. Ultradian oscillations are the focus of the model derived by Li et al., which accurately portrays the long-term behavior of the glucose-insulin regulatory system.

Ideally, a typical healthy human body maintains a plasma glucose concentration of 70 – 110 mg/dl. When the plasma glucose level falls bellow 70 mg/dl, the body enters into a state of hypoglycemia. In order to correct this the \alpha -cells in the pancreas are stimulated, secreting glucagon. This glucagon travels to the liver where it is converted into glucose and made available to the body, thus raising the low glucose concentration. It takes time for the liver to convert glucagon into glucose and this time delay is built into the mathematical model as \tau _{2}. Prolonged states of hypoglycemia are incredibly dangerous as they cut off the body's source of energy. This becomes dangerous as the brain is insulin-independent and cannot function without glucose. Diabetic comas and death occur when blood sugar becomes too low. When the plasma glucose concentration rises above 110 mg/dl, the body enters into a state of hyperglycemia. In response to hyperglycemia, the \beta -cells in the pancreas secrete insulin. It takes time for the secreted insulin to become "remote" insulin after being produced inside the \beta -cells. This time delay is built into the mathematical model as \tau _{1}. Once insulin is available it does three things: (1) it helps insulin-dependent cells consume the available glucose, (2) it stops hepatic glucose production and (3) it degrades and is cleared by the body. In this way the body is able to regulate elevated glucose concentrations. This behavior is best captured in Figure 1.


Figure 1- Two time delay glucose–insulin regulatory system model. The dash-dot-dot lines indicate that insulin inhibits hepatic glucose production with time delay; the dash-dot lines indicate insulin secretion from the \beta -cells stimulated by elevated glucose concentration level and the short dashed line indicates the insulin caused acceleration of glucose utilization in cells with time delay; the dashed lines indicate low glucose concentration level triggering a-cells in the pancreas to release glucagon.

Diabetes Mellitus develops when the body's regulatory system is unable to appropriately balance the body's glucose and insulin concentrations. Type I diabetes occurs when the pancreas loses the ability to secrete the appropriate amount of insulin to process the available glucose. Type I diabetics require supplemental insulin to manage the glucose levels. Type I diabetics frequently can enter into states of hypoglycemia if they inject too much insulin without eating. Type II diabetes is known as 'insulin resistant diabetes' in that the insulin-dependent cells in the body become resistant to insulin. More specifically, insulin-dependent cells, consisting of most cells in the body besides brain and nervous cells, require insulin to be able to utilize glucose. When these cells become insulin resistant they are not triggered as well by insulin to consume glucose. This results in a state of chronic hyperglycemia which can be treated through weight loss, dieting and oral medications. Type I and II diabetes are both considered in the mathematical analysis below.

Existing Model: Reproduced Results

For a discussion of existing models that preceded the one developed by Li et al., please see MBW:Glucose-Insulin Regulatory System. Below we briefly develop the model given by Li et al.

The following conservation model is the guiding principle behind the proposed model. Namely, that any change in glucose and insulin can be represented by the following formulas:

{\tfrac  {dG(t)}{dt}}= glucose production - glucose utilization

{\tfrac  {dI(t)}{dt}}= insulin production - insulin clearance

Utilizing previous mathematical work, Li et al. proposed the coupled delay differential equations seen below.

G'(t)=G_{{in}}-f_{2}(G(t))-f_{3}(G(t))f_{4}(I(t))+f_{5}(I(t-\tau _{2})); I'(t)=f_{1}(G(t-\tau _{1}))-d_{i}(I(t))

An important strength of this model is that each term in the differential system has some biological significance; previous models had to include artificial state variables to account for the delays inherent in the regulatory system. These additional terms did not have a direct physiological analog to produce results that matched those observed empirically. The equations for functions f_{1} through f_{5} are given below and the justification of their mathematical forms is described in MBW:Glucose-Insulin Regulatory System.


Their forms of f_{1} through f_{5} are shown below:

The following table lists the parameter values used in these functions:


Before reproducing the results from the paper by Li et al., however, one must be extremely careful in how the functional terms are treated. It should be noted that functions require G(t) to be in mg and I to be in mU, giving Failed to parse (lexing error): G’(t)

in mg/min and Failed to parse (lexing error): I’(t)
in mU/min.  In plotting the functions for G(t) and I(t), however, Li et al. convert the units of G and I from mg and mU to mg/dl and \mu U/ml.  After difficulty in reproducing the exact results, we consulted the authors and concluded that in converting from strict amounts of glucose and insulin into concentrations the authors divided by the total blood volume. Using their parameter values and solutions, we found that Li et al. most likely assumed a total blood volume of 10 liters. This assumption is unrealistic as it is double the typical blood volume observed in humans. However, in order to be consistent with Li et al. we maintain a 10 liter assumption. The solutions do not change in type for a different amount of blood volume, but produce different glucose concentrations for the same parameter values. Thus, the appropriate conversion can be completed by dividing the output G from the first DDE by 100 and the output I from the second DDE by 10.  However, G_{{in}} is given in units of mg/(dl*min), so this must be converted into mg before solving the first DDE; this conversion can be done by multiplying the given G_{{in}} by 100.  With these conversions in place, the system indeed solves to give the exact results observed by Li et al.  Below we show the reproduced results for the ultradian oscillations in a healthy human given constant enteral nutrition, with \tau _{1}=7 min, \tau _{2}=12 min, G_{{in}}=1.08 mg/dl/min and d_{i}=.06: 

Normal solution.jpg

We also reproduced the bifurcation plots from the paper by Li et al:

Note that on p732 of Li et al., there is a typo in the bifurcation point in the right-hand plot of Figure 9; the paper cites a bifurcation in the first time delay \tau _{1} at a value of 2.8, but the actual bifurcation point is at \tau _{1} 3.8 (this can even be verified in the paper itself; the figure shows that the bifurcation clearly occurs at some value closer to 4 than to 3). The bifurcation analysis was important in understanding appropriate ranges for the parameters in the delay equations. Because there are indeed utradian oscillations in the glucose-insulin regulatory system we would expect parameter values to occur in ranges that produce stable oscillations. The results given here are especially consistent for \tau _{1}, \tau _{2} and G_{{in}} with observations from experiments and glucose tolerance tests. The bifurcations plots were reproduced by fixing three of the parameters and varying one of the parameters. For each value of the varied parameter the long term solution behavior was analyzed for stable steady states. The code that produced these plots and all of the following results can be found attached at the bottom of this document. The implications and importance of the bifurcation analysis is more thoroughly discussed in the review wiki of the paper by Li et al. With these model results now reproduced, we will now describe our extensions of this model.

Model Extensions

The model proposed by Li et al. produced the results seen in the previous sections of this wiki. The model provided insight into the ultradian oscillations that have been well observed within the body. In addition to being a more simple and accurate model than those previous posed, the two-delay model allowed for the previously discussed bifurcation analysis of the four parameters present in the model: G_{{in}}, d_{i}, \tau _{1} and \tau _{2}. The bifurcation analysis provided ranges for each parameter that produced the physically observed phenomena. The model developed by Li et al. provided an excellent foundation to understand the glucose-insulin regulatory system. We attempt to build upon that foundation by varying some of the assumptions of Li et al., while also exploring type I and type II diabetes. These extensions provided useful insight into the structure of the model, the importance of dieting, and the mechanisms of diabetes.

Non-Constant History

Figure 5 - Varying History

The delay differential equation structure proposed by Li et al. requires an initial history for both glucose and insulin. The analysis originally performed by the authors assumed that the glucose and insulin history was constant. Initially, this assumption seemed worth relaxing. Because the main goal of the two-delay model was to reproduce the ultradian stable oscillations present in the glucose-insulin regulatory system, Li et al. began their modeling by assuming an overnight fast. It would seem reasonable to think that the initial history should contain the oscillations similar to those produced by the model. We began to extend their model by using various histories, some with oscillations and some without. It was quickly seen that the initial history had a negligible impact on the overall solution. Due to the very stable nature of the delay differential equations being studied, they produced very similar solutions for all initial histories. After some amount of time, every history converged to the same range of oscillations. The solutions are shifted in time as it takes awhile for the solution to respond from large perturbations, but any significant difference in the solutions was usually settled after the first 500 minutes. This can clearly be seen in Figure 5, where we examined the solutions for a hyperglycemic history, hypoglycemic history and oscillatory history. This result justifies the assumption Li et al. made in assuming a constant history. For this system of two-delay equations, the history has very little impact. Furthermore, this small extension informed us that the behavior of this system is mostly determined by the parameter values. This result is important because it shows that a healthy glucose-insulin regulatory system is stable for large perturbations in the glucose and insulin concentrations. This makes physical sense as a non-diabetic person is able to acutely consume large amounts of glucose without suffering any long term side effects. More specifically, even with large changes in plasma glucose concentration a healthy regulatory system will return to ultradian oscillations after some amount of time.

Meal Types

Figure 6 - Glucose Arrival for Typical Meal
Figure 7 - Glucose Arrival for Protein Meal
Figure 8 - Glucose Arrival for Three Regular Meals
Figure 9 - Three Normal Meals
Figure 10 - Three Protein Meals
Figure 11 - Six Normal Meals
Figure 12 - Six Protein Meals

Issues physically present with the glucose-insulin regulatory system often manifest themselves as either type I or type II diabetes. Living well with diabetes requires very careful and diligent management of diet and meals. Here we examine the impact of various meal forms on the plasma glucose and insulin concentrations. Originally, Li et al. assumed that the glucose being infused into the system, G_{{in}}, was constant. While the assumption was useful for the previously shown bifurcation analysis, we extend the work of Li et al. by adding in a meal structure. Islam, Leech, Lin and Chrostowski[3] modeled the intake of glucose following a typical meal. They used a system of four coupled ordinary differential equations to model glucose appearance following a meal. Because the focus of our work was primarily in how the two delay differential equation model would handle and account for meals we did not reproduce the exact solution given by Islam et al. Rather we fit their solution with a scaled gamma distribution with \alpha =4 and \beta =2. A typical meal takes the form shown in Figure 6. It is important to note that there is a large spike in glucose arrival that decays fairly quickly.

Over the past twenty years there have been a large number of people who adhere to protein diets or diets that limit the amount of simple carbohydrates available in a meal. Combining the ideas of Islam et al. and information from a review found on[4] we produced a glucose arrival curve for a protein-based meal that took the form in Figure 7, a scaled gamma distribution with \alpha =10 and \beta =2. It is important to note that the peak glucose arrival rate is much lower for the protein meal than the typical meal and peaks later than the typical meal.

The first change we make to the model proposed by Li et al. is to change G_{{in}} from a constant to a function that assumes normal meal intake at 7:00 am, 12:00 pm and 5:30 pm. For one day the glucose intake rate is seen in Figure 8, differing from the original work where G_{{in}} was constant. Simulating the glucose-insulin regulatory system for the different G_{{in}} produced the result seen in Figure 9.

The changes to the model show three large spikes in the glucose level following the meal intake. These spikes occur while the system is adapting to the perturbed glucose insulin rate. Temporarily the glucose reaches unsafe levels of hyperglycemia but the system rebounds down to typical levels. This change to Li et al.’s model produces a general idea of how the system responds to meals, but it is not all together correct as the body is better able to handle large intakes of glucose without creating such dangerous levels of hyperglycemia. We would expect the peak levels of glucose to occur below 200 mg/dl. However, the purpose of this model is to capture the ultradian oscillations that occur throughout the day and not the rapid oscillations that would more appropriately model these types of spikes. Yet, this model does allow us to to compare different types of meals. By comparing Figure 9 to Figure 10 it can be seen that eating three protein meals a day is in fact "easier" on the glucose-insulin regulatory system. The peak glucose concentration is lower and the system is not forced to recover from as large of oscillations. However, it can also be noted that it induces a higher level of insulin for longer. Similarly, comparing six smaller sized meals, in Figure 11, a day to three larger meals a day shows that it produces a more stable glucose concentration that is easier for the regulatory system to process. Figure 12 shows the solution curve for six smaller protein meals a day. Some current trends in fitness and dieting suggest that eating smaller meals, equally spaced, with fewer carbohydrates puts less demands on the pancreas and liver. While it is far from conclusive, the extensions discussed here suggest that dieting can have a large impact on the glucose regulatory system and that smaller meals more equally spaced with more protein might help the body maintain its target glucose concentration.

Parameter Correlation Factors

We next attempted to develop an extension to the model that would allow us to capture the behavior of diabetes within the system. Both types of diabetes will be discussed in detail in the next subsections of this wiki, but it is first important to better understand the structure of the model. We initially struggled to adapt to the original model to produce the elevated plasma glucose factor found in type II diabetes. To better understand the structure of the model and get a sense of which parameters were most important, we calculated correlation factors for each parameter as it pertained to the maximum steady state of glucose concentration achieved. This was done by creating a 10,000 Latin Hypercube Sample of the four parameters present in the model and calculating the maximum glucose concentration for each grouping of parameters. The correlation factors were calculated using the mathematical formula below.

correlation(x,y)={\frac  {cov(x,y)}{\sigma _{x}\sigma _{y}}}

We found two interesting results:

  • That the highest stable glucose concentration maintained by the parameter ranges set forth by Li et. al was 167 mg/dl.
  • The following correlation factors were calculated:
    • \tau _{1}=-0.3478, \tau _{2}=-0.3467, G_{{in}}=0.5103 and d_{i}=0.6339.

Because we were unable to produce glucose levels present in typical case of Type II diabetes, we expected the statistical analysis of the parameters to reveal important trends. The negative correlation values of the delays were expected because as the delays decrease the amplitude of the oscillations present in the regulatory system decrease because the body's response is more instantaneous. The value of G_{{in}} clearly impacts the peak glucose level as was seen in the analysis of different meals. However, we did not think that it was as important for our analysis of type II diabetes because it is not always the case that type II diabetes intake more sugar - they simply have a hard time processing the sugar present. This left us with the highest correlation factor insulin clearance or degradation, d_{i}. Because it was the largest correlation factor in magnitude and made the most sense to modify, we decided to expand the range of d_{i} from .01-.12 1/min to .01-5 1/min. We were unsure how much to expand the range of insulin degradation and so started with 5/min. We found the following results:

  • For the new range of parameters we found the maximum peak glucose level to be 551 mg / dl.
  • The following correlation factors were calculated:
    • \tau _{1}=-.0254, \tau _{2}=-.0116, G_{{in}}=.3595 and d_{i}=.8532.

The maximum blood sugar found was slightly higher than we expected but has been found to occur in most type II diabetic patients. It was fascinating to see how dependent the peak blood glucose sugar was on insulin degradation - a result discussed later in this wiki. A correlation factor of .8532 is incredibly high, implying a nearly linear relationship between high plasma glucose concentration and insulin degradation. This result implies that there should be a very high correlation between type II diabetes and high insulin degradation. The results of the statistical analysis were consistent with results found from the examination of different histories. The stability of this system to large perturbations suggested that the most significant changes in solution would occur from a change in parameters. The results of the statistical analysis are discussed further in the analysis of each type of diabetes as well as the recent developments section of this wiki.

Type I Diabetes

Figure 12 - Model Results for Type I Diabetes

Type I diabetes is characterized by the body's inability for the body to produce insulin on its own. This introduces the necessity of injecting external insulin into the bloodstream to regulate glucose levels over the course of the day. With new technology, the regulation of insulin levels for individuals with type I diabetes is becoming increasingly simpler; now, those with type I often use insulin pumps that remain attached all day to secrete a steady supply of insulin into their bloodstream. This steady stream is referred to as ‘basal insulin replacement’; it corresponds to the constant working of the \beta -cells in a healthy individual during the day. Before eating a meal, individuals with type I diabetes must inject a bolus of insulin to help deal with the coming glucose spike caused by the meal. About half of a diabetic’s insulin intake goes to basal intake, and the other half is reserved for bolus injection[5]. Depending on where in the body the bolus is injected (the abdomen, arm, and thigh are some of the most common locations), the insulin level in the blood stream peaks within 60 and 120 minutes after injection.

In our model, we assume that bolus insulin is injected 10 minutes before a meal, and that the insulin level in the blood spikes 120 minutes after injection. We approximate the level of insulin in the blood after a bolus injection with a normal curve, with the area under the curve summing to the total amount of insulin injected in that bolus. We set the basal intake rate at 1 mU/min during the day and at .4 mU/min from midnight to 6:00 am. For an individual with type I diabetes eating three meals a day, our model results are included in Figure 13.

Note that with external insulin, the body’s glucose levels are kept within reasonable range, though this range is wider than that of a healthy person by about 120 mg/ml, with this difference occurring on the upper end of the glucose spikes (the range for a healthy person eating three meals is about 60 mg/ml to 260 mg/ml, whereas for a type I diabetic according to this model the range becomes about 60 mg/ml to 380 mg/ml). It is interesting to note that the diabetic also injects less insulin than is secreted by the \beta -cells in a healthy individual; the tradeoff is a less efficient control of glucose levels, but it also makes sense that a diabetic might require less insulin, since the diabetic can preempt a meal by ‘intelligently’ inject the right amount of insulin, whereas the body can only react to already increased glucose levels, resulting in some inefficiency.

Type I diabetes can only be treated successfully with external dosages of insulin, but type II diabetes is less severe, and can often be controlled by diet, medication, or both. An evaluation of the treatment methods for type II diabetes is the subject of the next section.

Type II Diabetes

Figure 14 - Type II diabetes modeled by elevated d_{i}
Figure 15 - Type II diabetes modeled by increased resistance
Figure 16 - Treatment of type II diabetes with insulin secretion drugs
Figure 17 - Treatment of type II diabetes by suppressing hepatic production

Type II diabetes is a chronic condition that typically has late onset in patients. Type II diabetes presents itself through elevated blood sugars caused primarily by increased insulin resistance in the cells. Because of the importance of parameters associated with the stable blood glucose concentration we take two approaches to adapting the model to type II diabetes. We first utilize the information from the latin hypercube simulation to specify parameters that led to a chronic state of oscillatory hyperglycemia. This was accomplished by using an increased rate of insulin degradation, which was inspired by the correlation factors calculated earlier in the wiki and justified in the conclusions and recent development section. The result is seen in Figure 14. The second approach is to multiply the product of f_{3}f_{4} by an efficiency coefficient in the range of zero to one. From the mathematical section of this wiki, it can be recalled that the product of f_{3}f_{4} is the consumption of glucose due to insulin dependent cells. Insulin resistance impacts these cells function, specifically their glucose consumption. The result of this is seen in Figure 15, where a stable non-oscillatory hyperglycemic solution is created. Both of the graphs represent ways of modeling elevated blood glucose concentration. It is important to note that in Figure 15 the efficiency coefficient used was .6 corresponding to a 94% insulin resistance. The difficulty that we had reproducing a situation using type II diabetes with the range of parameters suggests that higher insulin degradation rates may be the best way to capture behavior of type II diabetes. The graph of Figure 14 shows a stable oscillatory solution occurring at chronically elevated glucose levels. The parameters used to create Figures 14 and 15 are given below:

  • Figure 14: \tau _{1}=7.55 min, \tau _{2}=38.88 min, G_{{in}}=1.03 mg/dl/min and d_{i}=.9167 1/min
  • Figure 15: \tau _{1}=7 min, \tau _{2}=12 min, G_{{in}}=1.08, d_{i}=.06 1/min and efficiency constant .06.

In addition to weight loss and dieting programs, type II diabetes is treated with oral medication. The benefit of dieting can be seen in the graphs from meal section of this wiki. Changing the structure of meal creates lower spikes in blood glucose concentration which can be significant for a regulatory system that is not functioning properly. There are five main oral medications used to treat type II diabetes: sulfonylureas, meglitinides, thizolidinediones and a-glucosidase. Sulfonylureas and meglitinides attack type II diabetes be increase insulin secretion. By multiplying the insulin secretion function f_{1} by various multiples we were able to see a significant drop in glucose concentration. For the exact situation present in Figure 14 we were able to bring the glucose down to normal levels by multiplying f_{1} by 5 as seen in Figure 16. While this doesn't fully capture the specific nature or efficacy of treatment it does show that increasing insulin secretion would indeed help a type II diabetic. Efficacy of all treatments is discussed later and subject of current research for mathematical biologists. Metformin suppresses hepatic glucose production limiting the amount of glucose produced internally. Metformin in combination with diligent dieting can significantly lower the amount of glucose infused into the plasma. By supressing hepatic glucose production to 15 percent by multiplying f_{5} by .15 we were able to impact the glucose concentration seen in Figure 14 to what is shown in Figure 17. Again, this is not designed to be an accurate representation of treatment, but rather show that treatment is indeed effective and could impact type II diabetes. Thizolideniones helps prevent insulin resistance on a cellular level essentially increasing the efficiency factor we introduced, results are not shown for that. Lastly, a-glucosidase prevents complex carbohydrates from being broken down in the small intestine, lowering the glucose intake rate. All of these treatments have been implemented and studied both mathematically and clinically by other researchers[6]. Our worked primarily was to show that we could produce a situation that effectively captured the overall behavior of type II diabetes and show that treatment would indeed work.

Conclusions and Recent Developments

The work presented in this wiki is of the upmost importance due to the increasing prevalence of both type I and II diabetes in the United States and the world. Through reproducing all the results from Li et al. and making the described extensions we have developed a thorough understanding for the ultradian oscillations that occur within the glucose-insulin regulatory system and the two-delay model that captures the nature of the stable oscillations biologically observed. The bifurcation analysis provided insight into the various ranges of parameters that are physically realistic as they produce stable oscillations. This result was important throughout the literature describing the glucose-insulin regulatory system. By accounting for non-constant histories we saw that the history was not significant in the two-delay model proposed by Li et al. and that the system of equations was stable to large perturbations. Modeling different meal types and structures and both types of diabetes gave insight into the behavior of the regulatory systems and the types of solutions produced. The most important result derived on this wiki is the large correlation factor between elevated glucose concentrations and insulin degradation. After arriving at this result, we found current literature that also suggests the importance of insulin degradation in type II diabetics. Specifically there was a medical study examining insulin degradation in rats and diabetic patients[7]. Moreover, there was another paper by Wang, Li, and Kuang[8] that changed the functional form of insulin degradation to accurately measure the efficacy of oral treatment for type II diabetics. Wang et al. changed the functional form of d_{i} because of its "importance". However, we were unable to find a paper that displayed conclusive evidence for the importance of d_{i}. It seems that most people are attempting to better quantify and study the biological and mathematical significance of d_{i}, making us very pleased that our statistical analysis, though simple, showed a very high correlation factor. There is further work to be done on modeling the entire glucose-insulin regulatory system, diabetes and treatment options. As the disease continues to increase in prevalence, research will continue to occur in this field.


  1. Li, Jiaxu, Yang Kuang and Clinton C. Mason, 2006, Modeling the glucose–insulin regulatory system and ultradian insulin secretory oscillations with two explicit time delays, Journal of Theoretical Biology vol. 242
  2. Jason Gale. Nov. 7, 2010, India’s Diabetes Epidemic Cuts Down Millions Who Escape Poverty, Bloomberg Markets Magazine
  3. Islam, Md. Shafiqul, James Leech, Charles C.Y. Lin and Lukas Chrostowski, 2007, Peak Blood Glucose Prediction Algorithm Following a Meal Intake, University of British Columbia
  4. Insulin and Glucose, MedBio Info,, accessed April 2011
  5. University of California, San Francisco, 2010, Calculating Insulin Dose, Diabetes Education Online: Diabetes Teaching Center at the University of California, San Francisco
  6. Joseph L. Evans and Robert J. Rushakoff, 2002, Oral pharmacological agents for type 2 diabetes: sulfonylureas, meglitinides, metformin, thiazolidinediones, a-glucosidase inhibitors, and emerging approaches, Directory of Continuing Medical Education ch14
  7. Hossein Fakhrai-Rad, Andrej Nikoshkov, Ashraf Kamel, Maria Fernström, Juleen R. Zierath, Svante Norgren, Holger Luthman, and Joakim Galli, 2011, Insulin-degrading enzyme identified as a candidate diabetes susceptibility gene in GK rats, Human Molecular Genetics vol.9, issue 14
  8. Haiyan Wang, Jiaxu Li, and Yang Kuang, 2011, Enhanced modelling of the glucose–insulin system and its applications in insulin therapies, J. of Biological Dynamics vol.3, issue 1, p22-38



%Stephen Kissler and Cody Cichowitz
%Appm 4390
%Code that runs numerical integration of two delay system

%Important parameters:
    %I = Insulin Level
    %G = Glucose Level

%Necessary functions:

%Key for 'parameters':
    %parameters(1) = Vg
    %parameters(2) = Rm
    %parameters(3) = a1
    %parameters(4) = C1
    %parameters(5) = Ub
    %parameters(6) = C2
    %parameters(7) = C3
    %parameters(8) = Vp
    %parameters(9) = Vi
    %parameters(10)= E
    %parameters(11)= U0
    %parameters(12)= Um
    %parameters(13)= beta
    %parameters(14)= C4
    %parameters(15)= Rg
    %parameters(16)= alpha
    %parameters(17)= C5
    %parameters(18)= tp
    %parameters(19)= ti
    %parameters(20)= td

%Use this one for analyzing treatment and history changes
function [] =  ToyPlot(tau1,tau2,G,d,efficent,treatment,ghxt,ihxt) 

%Use this one for bifurcation analysis prints steady states
% function[ming,maxg,mini,maxi] = ToyPlot(tau1,tau2,G,d) 

%Only use the variables you need-------
global tot_v
global Gin
global di
global solution
global e
global tx
global ghx
global ihx
ghx = ghxt;
ihx = ihxt;
tx = treatment;
e = efficent;
Gin = G;
di = d;
tot_v = 10; %liters

tau = [tau1, tau2]; %delays

%Solve delay differential equations:
sol = dde23(@ddexIGde,tau,@ddexIGhist,[0, 60*24]);

%glucose and insulin solutions
glucose_sol = sol.y(1,:);
insulin_sol = sol.y(2,:);
measure = length(glucose_sol);

%calculating steady states
ming = min( glucose_sol(1,(floor(.8*measure): measure)));
maxg = max( glucose_sol(1,(floor(.8*measure): measure)));
mini = min( insulin_sol(1,(floor(.8*measure): measure)));
maxi = max( insulin_sol(1,(floor(.8*measure): measure)));

%Plotting ---------------------
title('Glucose Insulin Plot')
ylabel('Glucose (mg / dl) / Insulin Concentration (\muU / ml)')

%Define history functions:

function[s] = ddexIGhist(t)
%global variables 
global tot_v
global solution
global ghx
global ihx

s = [ghx,ihx]; %for when history is input
% s = [220,45]; %hyperglycemic history
% s = [55, 15]; %hypo glycemic history

%Oscillatory history 
% s(1) = 22.34/2*cos(2*pi*t/110.92) +94.18 ;
% s(2) = 5.2240/2*sin(2*pi*t/110.88)+10.178;

%Define differential equations:
function[dydt] = ddexIGde(t,y,Z)
global tot_v
parameters = zeros(20,1);
%Define function paramters; store in an array (note key at top):
parameters(1) = 10;
parameters(2) = 210;
parameters(3) = 300;
parameters(4) = 2000;
parameters(5) = 72;
parameters(6) = 144;
parameters(7) = 1000;
parameters(8) = 3;
parameters(9) = 11;
parameters(10)= .2;
parameters(11)= 40;
parameters(12)= 940;
parameters(13)= 1.77;
parameters(14)= 80;
parameters(15)= 180;
parameters(16)= .29;
parameters(17)= 26;
parameters(18)= 6;
parameters(19)= 100;
parameters(20)= 36;

%Define other paramters:
global di
global Gin
global e
global tx

ylag1 = Z(:,1);
ylag2 = Z(:,2);

%complete equations. Get rid of e (efficency used to model insulin resistance),
%tx (used to model treatment) and ifunc(t) (used for Type 1 diabetes) as necessary
dydt = [gflux(t,Gin) - ftwo(y(1)*10*tot_v,parameters)/(tot_v*10)...
        - fthree(y(1)*10*tot_v,parameters)*1*e*ffour(y(2)*tot_v,parameters)/(tot_v*10)...
            + tx*1*ffive(ylag2(2)*tot_v,parameters)/(tot_v*10) ;
      ifunc(t) + fone(ylag1(1)*10*tot_v,parameters)/tot_v - di*y(2)];

Auxiliary functions

%Stephen Kissler and Cody Cichowitz

%This is a subfunction to 'ToyPlot.m'.  This function governs insulin
%secretion, and depends on past glucose levels.

function[out] = fone(G,par)
out = par(2) / (1 + exp((par(4) - G/par(1))/par(3)));

%Stephen Kissler and Cody Cichowitz

%This is a subfunction to 'ToyPlot.m'.  This function governs glucose use
%by the brain and nerve cells (insulin-independent users), which is a
%function of current glucose level only.

function[out] = ftwo(G,par)
out = par(5)*(1-exp(-G/(par(6)*par(1))));
%Stephen Kissler and Cody Cichowitz

%This is a subfunction to 'ToyPlot.m'.  This function governs part of the
%glucose use by muscle, fat, and 'other' cells (insulin-dependent users).
%As such, it is a function of glucose level only, to be paired with ffour,
%which is a function of insuln levels only.

function[out] = fthree(G,par)
out = G/(par(7)*par(1));
%Stephen Kissler and Cody Cichowitz

%This is a subfunction to 'ToyPlot.m'.  This function governs part of the
%glucose use by muscle, fat, and 'other' cells (insulin-dependent users).
%As such, it is a function of insulin level only, to be paired with fthree,
%which is a function of glucose levels only.

function[out] = ffour(I,par)
out = par(11) + (par(12) - par(11))/(1+exp(-par(13)*log(I/par(14)*(1/par(9)...
    + 1/(par(10)*par(19))))));
%Stephen Kissler and Cody Cichowitz

%This is a subfunction to 'ToyPlot.m'.  This function governs the hepatic
%release of glucose.  It depends on past insulin levels.

function[out] = ffive(I,par)
out = par(15)/(1+exp(par(16)*(I/par(8) - par(17))));


%Cody Cichowitz and Stephen Kissler
%Pick appropriate meal type

function out = gflux(t, gin)

out = gin; %constant

 % % % % Three normal meals a day----------------------------
% if t < 420
%     out = 0;
% elseif t>=1860
%     out = fitgamma(t-1860,4,2); 
% elseif t >= 1050
%     out = fitgamma(t-1050,4,2); 
% elseif t >= 720
%     out = fitgamma(t-720,4,2);
% elseif t >= 420
%     out = fitgamma(t-420,4,2);
% end
% % % ----------------------------------------------------

% %Protien for 3 meals.....
% if t < 420
%     out = 0;
% elseif t >= 2550
%     out = fitgamma(t-2550,10,2)+fitgamma(t-2160,10,2)+fitgamma(t-1860,10,2)...
%       +fitgamma(t-1050,10,2)+fitgamma(t-720,10,2)+fitgamma(t-420,10,2);
% elseif t >=2160
%     out = fitgamma(t-2160,10,2)+fitgamma(t-1860,10,2)+fitgamma(t-1050,10,2)...
%       +fitgamma(t-720,10,2)+fitgamma(t-420,10,2);
% elseif t>=1860
%         out =
%         fitgamma(t-1860,10,2)+fitgamma(t-1050,10,2)+fitgamma(t-720,10,2)+fitgamma(t-420,10,2);
% elseif t >= 1050
%     out =
%     fitgamma(t-1050,10,2)+fitgamma(t-720,10,2)+fitgamma(t-420,10,2); 
% elseif t >= 720
%     out = fitgamma(t-720,10,2)+fitgamma(t-420,10,2); 
% elseif t >= 420
%     out = fitgamma(t-420,10,2);
% %----------------------------------------------------

%Six smaller meals a day carbs-----------------------------
% if t < 420
%     out = 0 ;%- 10*1000/60/10/10;
% elseif t>= 1260
%     out = .5*fitgamma(t-1260,4,2);%- 10*1000/60/10/10;
% elseif t >= 1050
%     out = .5*fitgamma(t-1050,4,2) ;%-50*1000/60/10/10;
% elseif t >= 900
%     out = .5*fitgamma(t-900,4,2); %-50*1000/60/10/10;
% elseif t >= 720
%     out = .5*fitgamma(t-720,4,2) ;%-50*1000/60/10/10;
% elseif t >= 600
%     out = .5*fitgamma(t-600,4,2) ;%-50*1000/60/10/10;
% elseif t >= 420
%     out = .5*fitgamma(t-420,4,2) ;%-50*1000/60/10/10;
% end

% % Six smaller meals a day protien-----------------------------
% if t < 420
%     out = 0;% - 10*1000/60/10/10;
% elseif t>= 1260
%     out = .5*fitgamma(t-1260,10,2)+.5*fitgamma(t-1050,10,2)
%     +.5*fitgamma(t-900,10,2)+.5*fitgamma(t-720,10,2) +.5*fitgamma(t-600,10,2)...
%       +.5*fitgamma(t-420,10,2);
% elseif t >= 1050
%     out = .5*fitgamma(t-1050,10,2)
%     +.5*fitgamma(t-900,10,2)+.5*fitgamma(t-720,10,2) +.5*fitgamma(t-600,10,2)+.5*fitgamma(t-420,10,2);
% elseif t >= 900
%     out = .5*fitgamma(t-900,10,2)+.5*fitgamma(t-720,10,2)...
%     +.5*fitgamma(t-600,10,2)+.5*fitgamma(t-420,10,2);
% elseif t >= 720
%     out = .5*fitgamma(t-720,10,2)...
%     +.5*fitgamma(t-600,10,2)+.5*fitgamma(t-420,10,2) ;
% elseif t >= 600
%     out = .5*fitgamma(t-600,10,2)+.5*fitgamma(t-420,10,2) ;%-50*1000/60/10/10;
% elseif t >= 420
%     out = .5*fitgamma(t-420,10,2) ;%-50*1000/60/10/10;
% end

%Cody Cichowitz and Stephen Kissler
%function for scaled gamma distributions for protient and carb meals

function out = fitgamma(x,alpha,beta)

x = x/60;
out= ((beta^alpha)/gamma(alpha)) .* x.^(alpha-1) .* exp(-beta*x);
out = 78/.4481*out;
out = out* .3;


%Cody Cichowitz and Stephen Kissler
%Runs through biffurcation solutions for all parameters

clear all

tau2 = 12;
Gin = 1.08;
di = .06;
for z = 0.01 :.8 : 20; %loop through tau1
    j = j + 1;
    [ a, b, c, d] = ToyPlot(z,tau2,Gin,di);
    ming(j) = a;
    maxg(j) = b;
    mini(j) = c;
    maxi(j) = d;
    tau1(j) = z;
%Format and Plot result
tau = [tau1,tau1];
g = [ming,maxg];
i = [mini,maxi];
xlabel('\tau_1 (min)')
ylabel('Glucose (mg / dl) and Insulin (\muU / ml )')
title('Bifurcation Analysis for \tau_1')

clear all
tau1 = 7;
gin = 1.08;
di = .06;

for z = 0.1 :1 : 40; %loop through tau2
    j = j + 1;
    [ a, b, c, d] = ToyPlot(tau1,z,gin,di);
    ming(j) = a;
    maxg(j) = b;
    mini(j) = c;
    maxi(j) = d;
    tau2(j) = z;
%format and plot solution
tau = [tau2,tau2];
g = [ming,maxg];
i = [mini,maxi];
xlabel('\tau_2 (min)')
ylabel('Glucose (mg / dl) and Insulin (\muU / ml )')
title('Bifurcation Analysis for \tau_2')

clear all
tau1 = 7;
tau2 = 12;
di = .06;
for z = 0 :.05 : 1.5;
    j = j + 1;
    [ a, b, c, d] = ToyPlot(tau1,tau2,z,di);
    ming(j) = a;
    maxg(j) = b;
    mini(j) = c;
    maxi(j) = d;
    gin(j) = z;
gin = [gin,gin];
g = [ming,maxg];
i = [mini,maxi];
xlabel('G_{in} (mg / dl / min)')
ylabel('Glucose (mg / dl) and Insulin (\muU / ml )')
title('Bifurcation Analysis for G_{in}')

clear all
tau1 = 7;
tau2 = 12;
gin = 1.08;

for z = .01:.005 : .12;
    j = j + 1;
    [ a, b, c, d] = ToyPlot(tau1,tau2,gin,z);
    ming(j) = a;
    maxg(j) = b;
    mini(j) = c;
    maxi(j) = d;
    di(j) = z;
di1 = [di,di];
g = [ming,maxg];
i = [mini,maxi];
xlabel('d_i (1 / min)')
ylabel('Glucose (mg / dl) and Insulin (\muU / ml )')
title('Bifurcation Analysis for d_i')

Correlation Factors

%Cody Cichowitz and Stephen Kissler
%Appm 4390

%function that outputs correlation factors, need date in.
function [cor_tau1,cor_tau2,cor_gin,cor_di] = cor_factors(data)

z = [data(:,2:5),data(:,6)];

%Calculate correlations factors...
cor_tau1 = tt(1,2)/(sqrt(tt(1,1))*sqrt(tt(2,2)));
tt1 = cov(z(:,2),z(:,5));
cor_tau2 = tt1(1,2)/(sqrt(tt1(1,1))*sqrt(tt1(2,2)));
tt2 = cov(z(:,3),z(:,5));
cor_gin = tt2(1,2)/(sqrt(tt2(1,1))*sqrt(tt2(2,2)));
tt3 = cov(z(:,4),z(:,5));
cor_di = tt3(1,2)/(sqrt(tt3(1,1))*sqrt(tt3(2,2)));

%Cody Cichowitz and Stephen Kissler
%Appm 4390
%function that creates data to be put in cor_factors.m

function [solution] = stochastic()
clear all
close all

%creates Latin Hyper cube over parameters
inputs(:,1) = 20*inputs(:,1); %tau1
inputs(:,2) = 40*inputs(:,2); %tau2
inputs(:,3) = 1.5*inputs(:,3); %g_in
inputs(:,4) = ((.12-.01)*inputs(:,4)+.01); %di

solution = zeros(10000,9);

%creates solution
for k = 1 : length(inputs(:,1))
    solution(k,:) = [k, inputs(k,:),ming,maxg,mini,maxi];

Type I diabetes

%Stephen Kissler and Cody Cichowitz
%APPM 4390

%This function determines insulin input in the case of type I diabetes.
%There are spikes in insuiln that occur ten minutes before each meal; these
%spikes are defined by a normal curve included in 'normalfunc.m'.
%Otherwise, insulin levels are kept at a constant basal level, reflecting
%the action of an insulin pump.

function out = ifunc(t)

if t > 2510 && t < 2630
    out = normalfunc(t-2510); %Insulin spike, day 2 dinner
elseif t > 2150 && t < 2270
    out = normalfunc(t-2150); %Insulin spike, day 2 lunch 
elseif t > 1790 && t < 1910
    out = normalfunc(t-1790); %Insulin spike, day 2 breakfast
elseif t > 1070 && t < 1190
    out = normalfunc(t-1070); %Insulin spike, day 1 dinner
elseif t > 710 && t < 830
    out = normalfunc(t-710); %Insulin spike, day 1 lunch
elseif t > 350 && t < 470
    out = normalfunc(t-350); %Insulin spike, day 1 breakfast
elseif t > 350 && t < 1190
    out = 1; %Daytime insulin level, day 1
elseif t > 1790 && t < 2630
    out = 1; %Daytime insulin level, day 2
elseif t > 1190 && t < 1400
    out = 1.3; %Insulin level from last meal to midnight, day 1
elseif t > 2630 && t < 2880
    out = 1; %Insulin level from last meal to midnight, day 2
    out = .4; %Insulin level from midnight to morning
%Stephen Kissler and Cody Cichowitz
%APPM 4390

%This function is used in ifunc.m, and defines a bolus of insulin
%corresponding to an insuiln shot that would precede a meal for a type I

function[out] = normalfunc(x)

%This function gives a normal curve for the uptake of insulin after a
mu = 120;
sigma = 35;
TotAmnt = 2000;
out = TotAmnt/sqrt(2*pi*sigma^2) * exp(-(x-mu).^2./(2*sigma^2));