
MBW:Extensions to a 2Delay GlucoseInsulin Regulatory ModelFrom MathBioThe 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: GlucoseInsulin Regulatory System. ContentsExecutive summaryIn this wiki, we replicate and extend a model for the ultradian glucoseinsulin 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 glucoseinsulin regulatory system behavior under both type I and type II diabetes. The human glucoseinsulin 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 bloodglucoserelated 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 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 ContextWhile a complete biological context for the mathematical model is given in the review wiki found at MBW: GlucoseInsulin Regulatory System, we briefly review the biological framing of the glucoseinsulin regulatory system and type I and II diabetes. In order to function properly the human body aims to maintain a stable level of plasmaglucose concentration that enables cells to function appropriately. The glucoseinsulin 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 largescale response to glucose and insulin. Ultradian oscillations are the focus of the model derived by Li et al., which accurately portrays the longterm behavior of the glucoseinsulin 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 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 . Prolonged states of hypoglycemia are incredibly dangerous as they cut off the body's source of energy. This becomes dangerous as the brain is insulinindependent 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 cells in the pancreas secrete insulin. It takes time for the secreted insulin to become "remote" insulin after being produced inside the cells. This time delay is built into the mathematical model as . Once insulin is available it does three things: (1) it helps insulindependent 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 dashdotdot lines indicate that insulin inhibits hepatic glucose production with time delay; the dashdot lines indicate insulin secretion from the 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 acells 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 insulindependent cells in the body become resistant to insulin. More specifically, insulindependent 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 ResultsFor a discussion of existing models that preceded the one developed by Li et al., please see MBW:GlucoseInsulin 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: glucose production  glucose utilization insulin production  insulin clearance Utilizing previous mathematical work, Li et al. proposed the coupled delay differential equations seen below.
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 through are given below and the justification of their mathematical forms is described in MBW:GlucoseInsulin Regulatory System. Their forms of through 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 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 and , however, Li et al. convert the units of and from mg and mU to mg/dl and 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, 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 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 min, min, mg/dl/min and : 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 righthand plot of Figure 9; the paper cites a bifurcation in the first time delay at a value of 2.8, but the actual bifurcation point is at 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 glucoseinsulin regulatory system we would expect parameter values to occur in ranges that produce stable oscillations. The results given here are especially consistent for , and 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 ExtensionsThe 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 twodelay model allowed for the previously discussed bifurcation analysis of the four parameters present in the model: , , and . 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 glucoseinsulin 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. NonConstant HistoryThe 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 twodelay model was to reproduce the ultradian stable oscillations present in the glucoseinsulin 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 twodelay 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 glucoseinsulin regulatory system is stable for large perturbations in the glucose and insulin concentrations. This makes physical sense as a nondiabetic 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
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 medbio.info.com^{[4]} we produced a glucose arrival curve for a proteinbased meal that took the form in Figure 7, a scaled gamma distribution with and . 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 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 was constant. Simulating the glucoseinsulin regulatory system for the different 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 glucoseinsulin 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 FactorsWe 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.
We found two interesting results:
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 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, . Because it was the largest correlation factor in magnitude and made the most sense to modify, we decided to expand the range of from . 1/min to 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:
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 DiabetesType 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 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 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 DiabetesType 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 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 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 nonoscillatory 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:
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 aglucosidase. Sulfonylureas and meglitinides attack type II diabetes be increase insulin secretion. By multiplying the insulin secretion function 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 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 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, aglucosidase 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 DevelopmentsThe 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 glucoseinsulin regulatory system and the twodelay 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 glucoseinsulin regulatory system. By accounting for nonconstant histories we saw that the history was not significant in the twodelay 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 because of its "importance". However, we were unable to find a paper that displayed conclusive evidence for the importance of . It seems that most people are attempting to better quantify and study the biological and mathematical significance of , 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 glucoseinsulin regulatory system, diabetes and treatment options. As the disease continues to increase in prevalence, research will continue to occur in this field. Resources
CodeToyPlot.m%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: %fone(G) %ftwo(G) %fthree(G) %ffour(I) %ffive(I) %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  figure plot(sol.x,glucose,'r',sol.x,insulin,'b') title('Glucose Insulin Plot') xlabel('time(min)') ylabel('Glucose (mg / dl) / Insulin Concentration (\muU / ml)') legend('Glucose','Insulin') % % %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 %delays 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 (insulinindependent users), which is a %function of current glucose level only. function[out] = ftwo(G,par) out = par(5)*(1exp(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 (insulindependent 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 (insulindependent 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 %3.22.2011 %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)))); gflux.m%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(t1860,4,2); % elseif t >= 1050 % out = fitgamma(t1050,4,2); % elseif t >= 720 % out = fitgamma(t720,4,2); % elseif t >= 420 % out = fitgamma(t420,4,2); % end % % %  % %Protien for 3 meals..... % if t < 420 % out = 0; % elseif t >= 2550 % out = fitgamma(t2550,10,2)+fitgamma(t2160,10,2)+fitgamma(t1860,10,2)... % +fitgamma(t1050,10,2)+fitgamma(t720,10,2)+fitgamma(t420,10,2); % elseif t >=2160 % out = fitgamma(t2160,10,2)+fitgamma(t1860,10,2)+fitgamma(t1050,10,2)... % +fitgamma(t720,10,2)+fitgamma(t420,10,2); % elseif t>=1860 % out = % fitgamma(t1860,10,2)+fitgamma(t1050,10,2)+fitgamma(t720,10,2)+fitgamma(t420,10,2); % elseif t >= 1050 % out = % fitgamma(t1050,10,2)+fitgamma(t720,10,2)+fitgamma(t420,10,2); % elseif t >= 720 % out = fitgamma(t720,10,2)+fitgamma(t420,10,2); % elseif t >= 420 % out = fitgamma(t420,10,2); %end % % % % %Six smaller meals a day carbs % if t < 420 % out = 0 ;% 10*1000/60/10/10; % elseif t>= 1260 % out = .5*fitgamma(t1260,4,2);% 10*1000/60/10/10; % elseif t >= 1050 % out = .5*fitgamma(t1050,4,2) ;%50*1000/60/10/10; % elseif t >= 900 % out = .5*fitgamma(t900,4,2); %50*1000/60/10/10; % elseif t >= 720 % out = .5*fitgamma(t720,4,2) ;%50*1000/60/10/10; % elseif t >= 600 % out = .5*fitgamma(t600,4,2) ;%50*1000/60/10/10; % elseif t >= 420 % out = .5*fitgamma(t420,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(t1260,10,2)+.5*fitgamma(t1050,10,2) % +.5*fitgamma(t900,10,2)+.5*fitgamma(t720,10,2) +.5*fitgamma(t600,10,2)... % +.5*fitgamma(t420,10,2); % elseif t >= 1050 % out = .5*fitgamma(t1050,10,2) % +.5*fitgamma(t900,10,2)+.5*fitgamma(t720,10,2) +.5*fitgamma(t600,10,2)+.5*fitgamma(t420,10,2); % elseif t >= 900 % out = .5*fitgamma(t900,10,2)+.5*fitgamma(t720,10,2)... % +.5*fitgamma(t600,10,2)+.5*fitgamma(t420,10,2); % elseif t >= 720 % out = .5*fitgamma(t720,10,2)... % +.5*fitgamma(t600,10,2)+.5*fitgamma(t420,10,2) ; % elseif t >= 600 % out = .5*fitgamma(t600,10,2)+.5*fitgamma(t420,10,2) ;%50*1000/60/10/10; % elseif t >= 420 % out = .5*fitgamma(t420,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.^(alpha1) .* exp(beta*x); out = 78/.4481*out; out = out* .3;
bifurcation.m%Cody Cichowitz and Stephen Kissler %Runs through biffurcation solutions for all parameters %tau1 clc clear all tau2 = 12; Gin = 1.08; di = .06; j=0; 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; end %Format and Plot result tau = [tau1,tau1]; g = [ming,maxg]; i = [mini,maxi]; figure plot(tau,g,'ro',tau,i,'bo') xlabel('\tau_1 (min)') ylabel('Glucose (mg / dl) and Insulin (\muU / ml )') title('Bifurcation Analysis for \tau_1') legend('glucose','insulin') %tau2 clear all tau1 = 7; gin = 1.08; di = .06; j=0; 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; end %format and plot solution tau = [tau2,tau2]; g = [ming,maxg]; i = [mini,maxi]; figure plot(tau,g,'ro',tau,i,'bo') xlabel('\tau_2 (min)') ylabel('Glucose (mg / dl) and Insulin (\muU / ml )') title('Bifurcation Analysis for \tau_2') legend('glucose','insulin') %Gin clear all tau1 = 7; tau2 = 12; di = .06; j=0; 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; end gin = [gin,gin]; g = [ming,maxg]; i = [mini,maxi]; figure plot(gin,g,'ro',gin,i,'bo') xlabel('G_{in} (mg / dl / min)') ylabel('Glucose (mg / dl) and Insulin (\muU / ml )') title('Bifurcation Analysis for G_{in}') legend('glucose','insulin') %di clear all tau1 = 7; tau2 = 12; gin = 1.08; j=0; 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; end di1 = [di,di]; g = [ming,maxg]; i = [mini,maxi]; figure plot(di1,g,'ro',di1,i,'bo') xlabel('d_i (1 / min)') ylabel('Glucose (mg / dl) and Insulin (\muU / ml )') title('Bifurcation Analysis for d_i') legend('glucose','insulin') 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... tt=cov(z(:,1),z(:,5)); 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 clc close all %creates Latin Hyper cube over parameters inputs=lhsdesign(10000,4); 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)) [ming,maxg,mini,maxi]=ToyPlot(inputs(k,1),inputs(k,2),inputs(k,3),inputs(k,4)); solution(k,:) = [k, inputs(k,:),ming,maxg,mini,maxi]; end 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(t2510); %Insulin spike, day 2 dinner elseif t > 2150 && t < 2270 out = normalfunc(t2150); %Insulin spike, day 2 lunch elseif t > 1790 && t < 1910 out = normalfunc(t1790); %Insulin spike, day 2 breakfast elseif t > 1070 && t < 1190 out = normalfunc(t1070); %Insulin spike, day 1 dinner elseif t > 710 && t < 830 out = normalfunc(t710); %Insulin spike, day 1 lunch elseif t > 350 && t < 470 out = normalfunc(t350); %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 else out = .4; %Insulin level from midnight to morning end  %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 %diabetic. function[out] = normalfunc(x) %This function gives a normal curve for the uptake of insulin after a %bolus. mu = 120; sigma = 35; TotAmnt = 2000; out = TotAmnt/sqrt(2*pi*sigma^2) * exp((xmu).^2./(2*sigma^2)); <\pre> 