|Home | About | Journals | Submit | Contact Us | Français|
In insulin pump therapy, optimization of bolus and basal insulin dose settings is a challenge. We introduce a new algorithm that provides individualized basal rates and new carbohydrate ratio and correction factor recommendations. The algorithm utilizes a mathematical model of blood glucose (BG) as a function of carbohydrate intake and delivered insulin, which includes individualized parameters derived from sensor BG and insulin delivery data downloaded from a patient’s pump.
A mathematical model of BG as a function of carbohydrate intake and delivered insulin was developed. The model includes fixed parameters and several individualized parameters derived from the subject’s BG measurements and pump data. Performance of the new algorithm was assessed using n = 4 diabetic canine experiments over a 32 h duration. In addition, 10 in silico adults from the University of Virginia/Padova type 1 diabetes mellitus metabolic simulator were tested.
The percentage of time in glucose range 80–180 mg/dl was 86%, 85%, 61%, and 30% using model-based therapy and [78%, 100%] (brackets denote multiple experiments conducted under the same therapy and animal model), [75%, 67%], 47%, and 86% for the control experiments for dogs 1 to 4, respectively. The BG measurements obtained in the simulation using our individualized algorithm were in 61–231 mg/dl min–max envelope, whereas use of the simulator’s default treatment resulted in BG measurements 90–210 mg/dl min–max envelope.
The study results demonstrate the potential of this method, which could serve as a platform for improving, facilitating, and standardizing insulin pump therapy based on a single download of data.
Studies indicate that adult subjects with type 1 diabetes using a sensor-augmented pump (SAP) have lower hemoglobin A1c values than subjects undergoing multiple daily injection (MDI) therapy. This is achieved without an increase in the number of hypoglycemic events.1 In the landmark STAR 3 study,2 the advantages of SAP therapy hemoglobin A1c reduction were demonstrated over the control group receiving MDIs. Subjects in the SAP group were more likely to meet age-specific hemoglobin A1c targets and had lower area under the glucose concentration-time curve values for hyperglycemia with no increased risk of hypoglycemia. Glucose variability improved in the SAP group compared with the MDI group. However, the application of evidence from research to individualized care still remains a clinical challenge.3
The basis of SAP therapy is to emulate the pancreatic beta cells’ functionality of insulin secretion. A small amount of basal insulin is constantly secreted by the pancreas to maintain euglycemia. However, during meals, additional insulin is secreted to compensate for the increase in blood glucose (BG). Similarly, with SAP therapy, basal insulin is continuously infused with boluses provided to compensate for the amount of food intake. Boluses can be delivered instantaneously or can be programmed to be delivered over a prolonged period of time. Optimal management of pump therapy patients requires knowledge of insulin pharmacodynamics and pharmacokinetics in order to accurately estimate insulin-to-carbohydrate ratio (I:C) and correction factor (CF), information that is not available for the majority of the SAP users.4
The ultimate goal of diabetes management is the development of the artificial pancreas with fully automated insulin delivery controlled by corresponding algorithms. In attempts to reach this goal, a few promising algorithms have already been published,5–13 with some of the algorithms showing significant potential in clinical trials.14–19 However, there is still a significant way to go before any of these algorithms could be commercialized where the general diabetes population could benefit. Safety still remains the main obstacle for performing fully closed-loop insulin delivery, and this is mainly attributed to sensor inaccuracy and failure. However, currently available SAP therapy can be further improved by optimization based on personal data derived from SAP.
Currently, manufacturers of SAP systems do not provide estimates of either the basal insulin profiles or the I:C and CF. Physicians set up a basal pattern for each patient, and some use data from continuous glucose monitoring (CGM) devices. Consequently, the pattern and doses are adjusted to further improve the patient’s glucose profile. Traditionally, initial basal patterns for first-time pump users have been generated by applying the 1700/1900 “rule” or variation thereof.
Another common practice is to set the basal rates lower than needed and to compensate for the missing insulin through insulin boluses throughout the day. Health care professionals use a variety of heuristics in order to estimate a subject’s insulin therapy. For example, if a patient with a high basal/bolus insulin ratio demonstrates a tendency toward hyperglycemia, this would indicate that boluses were frequently missed, whereas frequent hypoglycemia may indicate a need for basal reevaluation.3
Typically, the clinicians’ first step is to estimate total daily insulin (TDI) requirements, which can be evaluated as a function of the subject’s body weight (BW).
A few publications have attempted to address the basal insulin adjustment issue. Palerm and coauthors21 introduced a recursive method, Wang and coauthors20 simulated an adaptive technique, and Miller and coauthors22 showed the ability of a fuzzy logic learning algorithm to adapt basal insulin rates. All these algorithms can potentially succeed in adjusting insulin therapy. However, the main inherent drawback is the recursive need for multiple adjustments that may burden patients and physicians.
We believe in causality in glycemic response to insulin and meal intake. This causality can be translated into a mathematical model that is able to predict the BG response adequately enough to deduce conclusions about the insulin therapy. If a model could generate glucose profiles that were predictive for a specific individual, it might ultimately be used to adjust therapy for that same individual.23 Herein we introduce an algorithm that uses a subject’s past insulin and glucose data to individualize (calibrate) a mathematical model that is then used to optimize the subject’s basal profile, I:C, and CF.
The algorithm first uses suitable data sets to identify the insulin and meal parameters of the mathematical model. Later, based on the mathematical model predictions, the subject’s optimal insulin therapy is estimated.
In general, physiological mathematical models can be categorized as empirical, semi-empirical, or fundamental. The first two approaches are mainly used in the attempt to model the human glucose-to-insulin relationships. Empirical models capture the system behavior through input–output data; Finan and coauthors,24 for example, used autoregressive models.25 Semi-empirical models are based on simplification of the human glucose metabolism. One ofthe best known human glucose models is the minimal model of Bergman and coauthors,26 which simplifies the human glucose system into a minimal number of compartments to model insulin sensitivity and glucose effectiveness. More fundamental models will require a thorough understanding of glucose metabolism on subcellular, cellular, and tissue levels that are currently unavailable.
In this article, we first introduce the mathematical model formulation followed by the optimization method used for evaluating insulin therapy. Next, we present the experimental and simulation results that were conducted in order to evaluate our method.
A mathematical model that enables individualization of a subject’s unique BG dynamics was formulated. The model describes the BG dynamics as a function of insulin and meal intake using three main compartments (Figure 1). It includes several unique features such as a biphasic meal model that is used to estimate the postprandial BG response that is superimposed on the BG response to insulin. The presented model was chosen for following reasons: it is linear, is physiologically based, and includes only parameters that have direct connection/understanding to measurable data (BG and insulin delivered).
The insulin pharmacokinetics is estimated by the Steil and coauthors27 study, and the pharmacodynamics are computed per subject.
It was observed that changes in glycemia after ingestion of a mixed meal remains above baseline for a longer period than predicted by the 120 min glycemic index.28 Therefore, the meal influence on BG is modeled as a BG biphasic time-delayed meal response. This is based on the assumption that food is digested in two delayed phases, which is manifested as a faster response near the meal intake followed by a second delayed response.
The following equation describes the subcutaneous insulin concentration as a function of the administrated insulin:
where ÎD is the delivered insulin represented as a pulse between time t1 and t2 (min), Îin (U/h) is the insulin intake between time t1 and t2 (min), and s represents the complex argument of the Laplace transform.
The pharmacokinetic model (Figure 1) that describes the plasma insulin’s (Ip) dynamics as a function of subcutaneous insulin and Ip initial conditions is described by Equation (3). This pharmacokinetic model was evaluated using in vivo data by Steil and coauthors:27
where ÎP is Ip in deviation form and ÎP0 and dIP0 are ÎP and derivative initial conditions, respectively.
All the insulin states are formulated in deviation form from given insulin value I0. In deviation form, I0 is subtracted from all the insulin states, and therefore, zero insulin in deviation form is equal to the absolute value of I0:
where x represents the different insulin states D, in, or P.
The following equation expresses the glucose kinetic model (Figure 1) in deviation form:
where Ĝ, KI, Ĝ0, dG0, τ1, and τ2 are BG in deviation form, the insulin gain, the BG initial conditions in deviation form, the BG derivative initial conditions, and two time constants, respectively.
The following two equations describe the contribution of the biphasic meal response to the BG output:
where Min is the meal intake in grams of carbohydrate per hour between time t1 and t2 (min),
where GM (mg/dl) is the BG increment due to meal consumption, Dmeal1,2 (min) are time delays between the time of the meal intake and the time each meal hump starts to appear in the BG, KM1,2 are the meal response gains, and tM1,2 are time constants.
Final BG is calculated by
where Gout,0 is the fasting BG (mg/dl) at a given insulin rate of I0. The mathematical model predicts that, if the given insulin rate I0 is delivered over a long fasting period, the resulting patient BG will be fasting BG.
The mathematical model includes nine parameters that are estimated per subject (Table 1).
Through the 1800 “rule,”29,30 the maximum drop in mg/dl for a unit bolus of rapid-acting insulin can be estimated as a function of the TDI. Using Equation (3), we can estimate that KI (Table 1) is approximately three times the drop of BG as evaluated by the 1800 “rule”:
where DpU1800 is the drop in mg/dl per unit of insulin as it is predicted by the 1800 “rule.”
I0 and Gout,0 are estimated by data collected over a few fasting periods (normally night time). Next, we chose a few meal responses to estimate the remaining eight parameters: τ1,2, KM1,2, τM1,2, and Dmeal1,2 [Equations (2)–(8)]. For each data set, in addition to the eight parameters, we also estimate Ĝ0 and dG0.
The parametric estimation is conducted in two steps. First, a rough estimation is accomplished by using genetic algorithm.31 On the second step, the solution obtained by the genetic algorithm is used as an initial guess to a local solver (Matlab, Fminsearch).
The parametric optimization cost function is based on the root mean square error between the CGM values and model predicted BG:
where CGMk is the CGM at time instant k and BGmodel,k is the model BG prediction at sampling time k.
The parametric optimization is conducted under the following constrains:
Parameter ranges were selected in order to be realistic. The time constants represent the time it takes the system’s step response to reach ~63% of its final (asymptotic) value. The initial conditions are also bounded by assuming that a CGM has maximum absolute relative distance of 14%:
The derivatives of the initial conditions are bounded to be less than a change of 5 mg/dl/min, which is an estimation of the maximum physiological BG change:
We estimate the optimal basal that will bring a subject’s BG to 120 mg/dl (IBasal) by using the insulin gain (KI), the estimation of the Gout,0, and its corresponding basal value (I0):
The meal insulin therapy is estimated by the minimization formulation shown in Equation (15). The minimization is using the model-predicted meal response to estimate the I:C, the CF, and the 6 h postprandial square bolus:
where SqB1,2,3 are three values that produce the 6 h postprandial bolus when multiplied by the meal (g·CHO/h) and added to the IBasal. The maximum basal, basalmax, is set by default to 4 U/h.
The final insulin therapy is a superposition of the IBasal, the insulin bolus at the time of the meal intake, and the postprandial insulin bolus as illustrated in Figure 2.
Six diabetic canines were chosen in order to estimate model parameters (two dogs were excluded in the early stage of the experiment). The dogs’ BG levels were measured over three 24 h periods during nominal insulin pump therapy that was set in advance by the veterinary facility personal and were used as the control group. The overnight insulin in the first control period was set to the nominal basal rate and to ±0.1 U/h of the nominal basal between 10:00 pm and 7:00 am in the other two control experiments. In the first control experiment, where no insulin therapy change was implemented at the overnight period, the blood samples were taken every hour. For the rest of the experiment (including the model-based pump therapy), blood samples were taken every 30 min. The four dogs’ glucose levels were controlled using the model-based pump therapy for 32 h from 11:00 pm until 7:00 am two days after.
Ten in silico adult subjects based on the Food and Drug Administration (FDA)-accepted University of Virginia (UVa)/Padova type 1 diabetes mellitus metabolic simulator32 were used to assess the optimization algorithm using insulin intake and CGM data incorporating actual characteristics of the Enlite™ glucose sensor,33 the latest in Medtronic glucose sensing technology, integrated with the glucose simulator. The identification data were generated by using a scenario of 60, 60, and 80 g carbohydrate meals given at 7:00 am, 12:00 pm, and 6:00 pm, respectively. Each subject’s simulation started at midnight of day 1 and ended at 8:00 am of day 2 (overall, 32 h of simulation).
In the next step, I:C, CF, SqB, and IBasal were adjusted by using the solely calibrated mathematical model, and only the final results were tested with the UVa/Padova simulator using three meals in the 24 h length scenario of 60, 60, and 80 g carbohydrate for meals provided at 7:00 am, 12:00 pm, and 6:00 pm, respectively.
Table 2 summarizes experimental results for 24 h duration studies for four dogs from 6:00 am to 6:00 am. It can be seen that none of the canines under the model-based therapy reached BG values lower than 60 mg/dl. The percentage of time glucose was in the 80–180 mg/dl range was 86%, 85%, 61%, and 30% using model-based therapy and [78%, 100%] (brackets denote multiple experiments conducted under the same therapy and animal model), [75%, 67%], 47%, and 86% for the control experiments for dogs 1 to 4, respectively. The BG at 6:00 am of the second day were 110, 116, 166, and 170 using the model-based therapy and [88, 84], [101, 82], 74, and 100 for the control days for dogs 1 to 4, respectively. The TDI requirements were ~49, 49, 34, and 27 U for the model-based therapy and ~49, 47, 26, and 24 U using the nominal therapy for dogs 1 to 4, respectively. The basal/bolus rations were ~51%, 49%, 40%, and 50% applying the model based-therapy and 65%, 61%, 63% and 57% using the nominal therapy for dogs 1 to 4, respectively. This percentage was calculated taking into account the postprandial square bolus as meal bolus after the overnight basal was subtracted.
Table 3 summarizes the mathematical model parameters that were identified based on the experimental data. It can be seen that dog 1 is less sensitive to insulin, while dog 2 is relatively insulin resistant and also has the slowest pharmacodynamics—both τ1 and τ2 are 150 min (large values of time constants indicate on potentially slow pharmacokinetics). Dogs 3 and 4 have relatively higher sensitivity to insulin and faster pharmacodynamics—τ2 approximately 0.1 min. All dog identifications included a second meal hump delay response between 223 and 300 min, which indicates a very slow digestion for these specific dogs; however, this is specific to these dogs and does not apply to other dogs or to humans.
Figures 3 to 6 depict the 24 h experimental results for dogs 1–4, respectively. All dogs were provided with 13.2 oz. of canned dog food mixed with dog kibble (the amount of kibble is set per dog). The BG, the insulin delivered (in logarithmic scale), and the meal intake are shown in subplots (A) to (C), respectively. The brown, blue, and black curves are the model-based therapy experiment, the model-based therapy predicted outcome, and the control experimental days, respectively. The grey curves indicate outliers that were taken from the final statistics.
Figure 7 shows the control variability grid analysis34 plot that compares the model-based therapy and the simulator’s default treatment algorithm, blue circles and pink rectangles, respectively, using a 10 adult population of the UVa/Padova simulator with a meal scenario of 60, 60, and 80 g carbohydrate meals given at 7:00 am, 12:00 pm, and 6:00 pm, respectively. Both therapies for all virtual adult patients resulted in zone A and zone B.
Table 4 shows a comparison between the I:C and basal rates obtained by the model-based therapy compared with the simulator’s default treatment algorithm I:C and basal rates.
Six diabetic dogs were chosen in order to estimate model-based therapy algorithm parameters. Two dogs were later excluded from the experiment. One was emitting its meals, and therefore it was hard to estimate its meal intake, and the second was constantly removing the insulin pump catheter. All the dogs in at least one of the control experiments and in some of the model-based therapy experiments showed a significant rise in BG in the afternoon that sometimes lasted until midnight. This nonconsistent rise in BG can be a result of the dogs being taken for exercise, and although closely supervised, they might have eaten their own feces or temporarily disconnected the insulin pump catheter. We omitted control days 1, 1, 2, and 2 for dogs 1 to 4, respectively. In retrospect, since we already estimated the dogs’ CF, we should have given the dogs a correction bolus when we noticed the unexplained rise in BG.
We repeated the model-based therapy for dog 2 because the personnel who were conducting the experiment decided to suspend insulin delivery for 10 min when BG reached 61 mg/dl. The experiment was resumed after 10 min and carried over to the second day’s breakfast, where we did not observe a fast postprandial drop in BG. However, we decided to repeat the experiment, and when the postprandial response was very good without any dangerous drops in BG, however, we experienced a sudden elevation in BG in the afternoon after the dog was taken to exercise, as was observed with some of the control experiment days.
In retrospect, we probably should not have used the data sets that included the afternoon BG elevation in the identification procedure. This is especially noticeable for dogs 3 and 4, where the model predicts relatively high delayed meal responses in the afternoon as a result of using two of three control experiment days for training the model. Interestingly, the model-delayed meal hump for dog 4 predicts the experiment results for the model-based experiment in very good agreement. Three of the four experiments days with dog 4 ended with an afternoon BG rise, while one of the control days resulted in a relatively low BG profile over the most of the experiment, which made it hard to understand what the outliers are.
We demonstrated that, even with a less accurate mathematical model, the model-based therapy resulted in greater percentages of time spent at glucose levels of 80 to 180 mg/dl that are equally comparable to the nominal therapy for dog 1 and better for dogs 2 and 3. Dog 4 ended in an elevated BG, but it reached a lower maximum value (269 mg/dl) than two of the control periods that were omitted (391 and 314 mg/dl). The relatively elevated BG when the model-based therapy was applied for dog 4 may also arise from setting the basal rate at a too low level as a result of incorrectly evaluating the dog insulin sensitivity and BG.
We target our insulin therapy to 120 mg/dl, and the BG measured at 6:00 am after 31 h of model-based insulin therapy were 110, [116 123], 166, and 170 for dogs 1–4, respectively. However, dog 3 did reach the value of 126 mg/dl at 11:30 pm and stayed in the 131–117 mg/dl range until 2:00 am, when it started drifting to higher BG values. This study was a proof of concept, and many of the problems in the reported results were solved in the algorithm version used at the day of the publication.
The model-based therapy resulted in unchanged TDI requirements for dogs 1 and 2 when compared with the nominal therapy. For dogs 3 and 4, the TDI requirements were increased by 31% and 12%, respectively. The model-based therapy, in comparison with the nominal therapy, delivered more insulin to all dogs to compensate for the meals and less insulin for the basal.
The algorithm was tested initially on the simulator in order to assess its ability in optimal conditions. The indication that it works in silico and under ideal conditions is the first step in testing it in more realistic conditions. The UVa/Padova simulator is designed with no meal intraindividual variability, and hence the use for multiple meals is nonrealistic.
The algorithm was tested with 10 in silico adult subjects using the UVa/Padova FDA-accepted simulator. The population response BG values using the algorithm andthe simulator’s default treatment algorithm were in the 71–236 mg/dl range and in the 90–210 mg/dl min–max envelope, with a median of 121 and 127 mg/dl and interquartile of 26 and 23 mg/dl, respectively. The percentage of time at glucose levels of 70 and 180 mg/dlwere medians of 100% and 98% and interquartiles of 2% and 6% for the simulator’s default treatment and the model-based algorithm, respectively.
In order to test the robustness of the system to meal uncertainties, the simulations were repeated with 20% of overestimations and underestimations of the meal carbohydrate count. The population response BG values using 20% overestimations and underestimations of the meal carbohydrates were in the 57–229 mg/dl range and in the 86–242 mg/dl min–max envelope, with medians of 120 and 128 mg/dl and interquartiles of 23 and 32 mg/dl, respectively. The percentage of time at glucose levels of 70 and 180 mg/dl were medians of 100% and 98% and interquartiles of 4% and 10% for the 20% overestimations and underestimations of the meal carbohydrate count, respectively. These results indicate of the robustness of the method to meal uncertainties, at least in the simulated studies.
Continuous glucose monitoring data often include calibration errors and noisy measurements; however, the method we introduced is based on offline calculations that enable us to reduce calibration errors significantly by retrospective calibration as well as to reduce noisy data through filtering and by selecting less noisy data sets.
Our model-based algorithm was able to accurately predict and adjust the insulin therapy of all 10 in silico adult subjects by maintaining the BG levels at a safe glycemic zone. Only CGM with the Enlite sensor noise collected from the simulator was used for the offline mathematical model calibration and calculation of the insulin therapy, which is equivalent to using a single sensor download of data of a subject wearing a SAP to adjust his/her insulin therapy.
We introduced a model-based insulin therapy algorithm that can serve as a platform to optimize basal insulin profiles, I:C, and CF based on data derived from a patient’s SAP. The model is designed for routine BG control and will not handle sudden BG excursions from infection, steroid intake, or starting a diet.
It should be noted that we compared our algorithm with two relatively successful insulin therapies; the dog insulin therapy was set by professional laboratory technicians who had a relatively long time to set and reset the dog insulin therapy, and the simulator’s default treatment was optimized by the inventors of the simulator, assuming they had access to all the equations in the simulator’s model. Nevertheless, our model-based therapy was comparable to them.
We assume that, currently, a large number of insulin-dependent diabetes patients are given less than optimal insulin treatment, but even standardizing the SAP therapy has a merit per se. We can envision that, in the near future, diabetes patients around the world will be able to get effective insulin therapy by using a short period of SAP data. These recommendations can later be modified by the patient and his physician to further improve insulin therapy. This algorithm can be used in subjects who utilize SAP therapy regularly, as well as in SAP naïve subjects who will be put on a pump for a dedicated time in order to set their insulin therapy.
As it was manifested by the divergence in the experiment results, animal models provide a richer environment to test algorithms than computer simulations. The unexplained rise in BG in the afternoon occurred only in some of the experiments, with the same dog model having exactly the same food and exactly the same insulin treatment. This is part of the challenge when attempting to regulate BG in reality.
Having accurate insulin pump settings is necessary for effective BG regulation. However, the subject’s awareness of possible changes in BG is still a crucial factor in regulating his/her BG, and this need has great influence on the individual’s lifestyle. In the near future we can expect that closed-loop algorithms that regulate insulin basal rates will be approved. This model or other model-based therapy will complement these algorithms.
This work was funded by Medtronic Minimed Inc.
All authors are full-time employees of Medtronic Minimed Inc.