- Email: [email protected]
Diabetic Blood Glucose Control via Optimization over Insulin and Glucose Doses
8th IFAC Symposium on Biological and Medical Systems The International Federation of Automatic Control August 29-31, 2012. Budapest, Hungary
Diabetic Blood Glucose Control via Optimization over Insulin and Glucose Doses ⋆ Meike Stemmann ∗ Rolf Johansson ∗∗ ∗
Lund University, Sweden (e-mail: [email protected]). Lund University, Sweden (e-mail: [email protected])
∗∗
Abstract: The proposed algorithm aims at helping diabetic patients which use multi-dose injection (MDI) therapy to determine the dose and time of insulin and glucose intakes to stabilize their blood glucose (BG) concentration. The objective is to spend most of the time in normal BG range (70-180 mg/dl). The control was done by minimizing the risk of hypo- and hyperglycemia over the doses and times of insulin and glucose intakes. Keywords: Biomedical Control, Diabetes, Optimization 1. INTRODUCTION
for CSII patients, MDI patients need single insulin advices just a couple of times per day.
Diabetes Mellitus Type I (DMTI) is a chronic disease characterized by the body’s inability to regulate its own blood glucose concentration. In Type 1 diabetes, this is caused by a lack of insulin secretion resulting from a failure of the β-cells in the pancreas, due to either injuries or diseases. Viral infections and autoimmune disorders have both been associated with the destruction of the β-cells (Guyton and Hall 2006). Hence, patients with DMTI have to administer external insulin to substitute the lack of insulin production by the β-cells. This is done either by multiple daily injections (MDI) or through continuous subcutaneous insulin infusion (CSII) using an insulin pump. MDI patients usually take a single dose of basal insulin on daily basis to supply the body with the basic insulin needed and several additional bolus insulin doses at times when the concentration of glucose in the blood is high, e.g., around meal times. A diabetic patient has to determine the correct doses of insulin needed to stabilize the blood glucose (BG) concentration, and has to solve an optimization problem every day.
To fit this regime, the control algorithm proposed in this paper aims at determining the doses and time points for the insulin and also glucose administration by solving a nonlinear optimization problem. While insulin treatment is important in most of the everyday situations, additional glucose intakes can be important to prevent hypoglycemia under special conditions like exercising or stress (Guyton and Hall 2006). Therefore, in addition to insulin advices, also glucose advices are determined by the proposed control algorithm.
To help the patient with this task, blood glucose prediction algorithms (St˚ ahl and Johansson 2009, Cescon et al. 2009) as well as many different control algorithms have been proposed. These control algorithms reach from proportionalintegral-derivative controller (PID), pole placement over adaptive and run-to-run methods to model predictive control (MPC) (Cobelli et al. 2009, Harvey et al. 2010, Parker et al. 2001). Many of the proposed controllers aim at having a continuous insulin signal, that could be used by an insulin pump to continuously infuse exogenous insulin into the diabetic patient.
2. AN OPTIMIZATION-BASED CONTROL ALGORITHM
However, in the scope of the European project DIAdvisor (Poulsen et al. 2010, DIAdvisor 2012) as a blood glucose prediction and treatment advisory system, the patient should be able to use either CSII or MDI treatment. Whereas continuous insulin administration can be useful ⋆ FP7 IST-216592 DIAdvisor
978-3-902823-10-6/12/$20.00 © 2012 IFAC
40
This paper is organized as follows. First, the overall control algorithm is presented in the beginning of Sec. 2. Then the subsections 2.1, 2.2 and 2.3 present the core part of the control algorithm, which is the optimization problem determining the insulin and glucose advices. Section 3.1 shows the simulation of the control algorithm using virtual patients and Sec. 3.2 shows the evaluation of the control algorithm with data measured on real patients.
The proposed control algorithm determines the doses of insulin and glucose together with the times they should be administered to the patient by solving an optimization problem. This optimization is included in an algorithm, which invokes the optimization when the blood glucose concentration leaves a certain range, see Fig. 1. When the blood glucose concentration y falls under 90 [mg/dL] and at least 120 [min] have passed since the last intake, or it falls under 80 [mg/dL] and at least 15 [min] have passed since the last intake, the optimization problem is solved to determine time and dose for an insulin or glucose intake. Similarly, when the BG concentration rises over 130 [mg/dL] and at least 120 [min] have passed since the last intake, the optimization problem is solved to determine the doses. 10.3182/20120829-3-HU-2029.00069
blood glucose concentration y [mg/dL]
8th IFAC Symposium on Biological and Medical Systems August 29-31, 2012. Budapest, Hungary
The parameters n1 to n6 were estimated individually for the patient to be controlled. This was done using nonlinear constrained optimization, solving the optimization problem (2). The measured blood glucose concentration is denoted by y and the desired reference BG concentration to be reached by yref . The parameters n1 and n4 are constrained to ensure the correct gain of the respective model.
130
90
minimize
||yg (t, tg ) + yi (t, ti ) + yref − y||2
80
subject to
n1 > 0 n4 < 0
n1 ...n6
Figure 2 shows yg (t, tg , ug ) and yi (t, ti , ui ) with ug = 10 g and ui = 1 unit at times tg = ti = 100 as inputs for a sample patient. When both insulin and glucose are taken at different times, those two functions are added to form the total change of blood glucose concentration.
15 120 time since last intake [min]
yg [mg/dL]
Fig. 1. Optimization-based control algorithm: The filled areas represent when the optimization problem is solved in the space of time since last intake vs. BG concentration. The blue, vertically striped area for hyperglycemia and the green, horizontally striped area for hypoglycemia.
yi [mg/dL]
The proposed control algorithm solves the optimization to determine insulin or glucose doses when the blood glucose concentration leaves the range between 90 and 130 [mg/dL], in contrast to the previously mentioned normal range being between 70 to 180 [mg/dL]. This is done because the goal is to keep the BG concentration inside the normal range for as long as possible (Kovatchev et al. 2005). This section is organized as follows. First, a model describing the dynamics of a diabetic patient is described in Sec. 2.1. This model relates the glucose or insulin intakes to the change of BG concentration and is needed for the optimization. Next, the cost function used for the optimization problem is presented in Sec. 2.2. The optimization problem to be solved to determine doses and times of insulin and glucose administration is formulated in Sec. 2.3. This optimization problem uses the cost function and the patient model mentioned above.
To describe the change in blood glucose concentration as a response to an intake of insulin and glucose, a nonlinear patient model is used (Trogmann et al. 2010b, Trogmann et al. 2010a). This model takes the doses and times of glucose or insulin intakes as inputs and has the change of blood glucose concentration as an output. Eq. (1) shows this model for insulin and glucose intakes, where yg (t, tg , ug ) is the change in blood glucose concentration as response to a glucose intake of size ug [g] at time tg [min] and yi (t, ti , ui ) the same for an insulin intake with the size ui [IU] at time ti [min]. e−n2 ·(t−tg ) · (t − tg )n3 · ug 1 + e−0.5·(t−tg ) e−n5 ·(t−ti ) · (t − ti )n6 · ui yi (t, ti , ui ) = n4 · 1 + e−0.5·(t−ti )
40
yg (t, 100, 10): Change in BG for 10 g glucose
30 20 10 0 0
2 0 −2 −4 −6 −8 0
500
yi (t, 100, 1):
500
1000
time [min]
1500
2000
Change in BG for 1 unit insulin
1000
time [min]
1500
2000
Fig. 2. Output of the nonlinear model for a sample patient. Upper panel: Change of BG as response to 10 g of Glucose intake at time tg = 100. Lower panel: Change of BG as response to 1 unit of insulin intake at time ti = 100.
2.1 Glucose and Insulin Dynamic Model
yg (t, tg , ug ) = n1 ·
(2)
2.2 Cost Function For a diabetic patient, the risk connected to hypoglycemia is much larger than the risk connected to hyperglycemia. The proposed cost function for the optimization takes this into account through its asymmetric shape, so that too low BG values imply more risk than too high BG values. In this way, the higher risk of low BG compared to high BG is expressed. The advantage of an asymmetic over a quadratic cost function was already presented in (Kirchsteiger and del Re 2009) and (Dua et al. 2009) for example. The asymmetric cost function used here is shown in Eq. (3) (Kovatchev et al. 2005).
(1) R(y) = 10−3 · (y − yref )2 · (c1 + c2 · 30 · exp(−0.05 · y)) (3)
41
8th IFAC Symposium on Biological and Medical Systems August 29-31, 2012. Budapest, Hungary
The BG concentration is denoted by y and the desired reference BG by yref . The parameters c1 and c2 allow to shape the slopes of the asymmetric cost function.
subject to
70 ≤ yP ≤ 450 ui < 20 ug < 80 The reason for the large upper bound is to have a feasible solution of the optimization problem after meals.
Figure 3 shows the cost function for different values of c1 and c2 . The black dashed line shows the default with c1 = 1 and c2 = 1. The effect on the cost function of increasing or decreasing these parameters is shown for some examples. However, not loose the asymmetric shape it should be kept c2 ≥ 0.4 and c1 ≤ 4. 300
The minimization in the optimization problem (5) is done over the doses and times of insulin and glucose as optimization variables using the Optimization Toolbox from Matlab (Mathworks 2011). 2.4 Data
c1 = 0.1, c2 = 1 c1 = 2, c2 = 1
250
The patient data used in Sec. 3.2 below was collected during an in-hospital period at the Centre Hospitalier Universitaire de Montpellier, in the framework of the European FP7-IST research project DIAdvisorT M , see (DIAdvisor 2012) and (Poulsen et al. 2010). The blood glucose concentration was measured with a continuous blood glucose measurement device (CGMS) (Dexcom 2012). Furthermore, for comparison the blood glucose was measured by the Yellow Springs Instruments (YSI) blood glucose analyser (YSI 2012) and through self-monitored of blood glucose (SMBG). Also the insulin doses taken by a patient were recorded, both for patients using an insulin pump and for patients using multi-dose injection (MDI). The example patient taken for evaluation in Sec. 3.2 is an MDI patient.
c1 = 4, c2 = 1 c1 = 1, c2 = 0.1
200
c1 = 1, c2 = 0.4 c1 = 1, c2 = 2
R
150
c1 = 1, c2 = 3 c1 = 1, c2 = 1
100 50 0 0
50
100
150
200
blood glucose concentration y [mg/dL]
250
Fig. 3. The asymmetric cost function for different values of the parameters ci , i = 1, 2.
3. SIMULATION 3.1 Simulation using a Virtual Patient
2.3 Optimization Problem
In order to simulate the controller in a closed-loop for testing and evaluation, an implementation of a virtual patient is used. The virtual patient is a nonlinear physiologically based model (Dalla Man et al. 2007).
The optimization problem needed to be solved for determination of the doses and times for insulin and glucose intake is given in (5). It uses Eq. (3) as the cost function to be minimized. The BG concentration y in the cost function is replaced with the total BG concentration in the future horizon Hp , called yP , using Eq. (1), as shown in Eq. (4). yP (ui , ug , ti , tg , t) = yg (t, tg , ug ) + yi (t, ti , ui ) + yˆ(t) (4) Here, yˆ(t) is the predicted BG concentration for the future horizon Hp assuming no future insulin or glucose intakes. In this way, yP (t) together with the cost function depend on the doses and times for insulin and glucose administration.
The simulation set-up is shown in Fig. 4. The virtual patient takes glucose and insulin as inputs and gives out the measured CGMS signal y(t). All these signals are used by a prediction algorithm to predict the CGMS over a future horizon. The predictor uses a linear state-space model of the patient, here a virtual patient, and a Kalman filter to calculate the predictions (Cescon et al. 2009). The CGMS predictions and the measured CGMS are used by the controller to determine the doses of insulin and glucose to be given to the virtual patient, as described in Sec. 2.
The idea is to minimize the risk connected to a certain amount of measured blood glucose concentration y over the doses and times of glucose and insulin intakes. Furthermore, the optimization problem (5) penalizes too large doses of insulin or glucose in order to keep them as small as possible. The amount of insulin ui and glucose ug to be given as well as the predicted total BG concentration y are constrained. minimize ui ,ug ,ti ,tg
Hp X
The simulation result for three example virtual patients are shown in Fig. 5 and Fig. 6. In the former, the first four panels show the results for virtual patient 1 and the last four panels for virtual patient 2. For each patient, the first panel shows the measured CGMS signal y at the output of the virtual patients and the 180 mg/dL limit of the normal blood glucose range. The insulin ui suggested by the controller is shown in the 2nd panel, while the glucose ug suggested by the controller is shown in the 3rd panels. The glucose intake from simulated meals M is shown the 4th panel. The red impulse in the 2ndpanel around the time 686 is an insulin disturbance, which is not suggested by the controller but acts on the virtual patient. This is introduced to test the controllers ability to deal with dropping blood glucose, since the virtual patient
10−3 (yP (ui , ug , ti , tg ) − yref )2 (5)
i=1 · c1 + c2 · 30 · e−0.05·yP (ui ,ug ,ti ,tg ) +c3 · |ui | + c4 · |ug |
42
8th IFAC Symposium on Biological and Medical Systems August 29-31, 2012. Budapest, Hungary
Predictor
y [mg/dL]
Virtual Patient
y(t)
400
ui (t) ug [g]
ug (t) Controller
200 100
00:00
M [g]
In addition to closed-loop simulations with a virtual patient, the proposed control algorithm was evaluated using data measured on a real patient. A patient-specific model was estimated according to section 2.1 from this data. Since the data is already measured, it is not possible to evaluate the control algorithm in a closed-loop manner. The evaluation should instead consider weather or not the extra insulin or glucose intakes suggested by the controller are reasonable. The CGMS data measured y by the Dexcom device is fed into the controller input. The advices of insulin ui and glucose ug suggested by the controller are then recorded and compared with the measured CGMS data and actual applied insulin ia and glucose intakes M , that the patient took during the in-hospital stay. Results from this comparison for a sample patient are shown in Fig. 7 and Fig. 8 for two different days. On the first day for example (Fig. 7), the BG concentration still increases after the second meal, so that an additional insulin intake could make sense. On the second day, shown in Fig. 8, the BG concentration decreases into hypoglycemic range, so that the suggested extra carbohydrates appear reasonable. 4. DISCUSSION The proposed control algorithm was designed specifically for MDI patients, which require that advices on insulin boluses are given only a few times per day. Although in many approaches, MPC is used to produce a continuous insulin signal (De Nicolao et al. 2011, Percival et al. 2011), it is possible to use the MPC scheme to produce the required insulin advices for MDI patients through approximation of the insulin signal from the MPC (Kirchsteiger an del Re 2009, Kirchsteiger et al. 2009). The
43
18:00
0:00
06:00
10:00
18:00
0:00
06:00
10:00
18:00
0:00
06:00
10:00
18:00
0:00
20
100 50 time [h]
y [mg/dL]
300
BG concentration of virtual patient 2
200 100 00:00
ui [units]
3.2 Open-Loop Simulation with Real Patient Data
10:00
40
00:00
ug [g]
After a meal, the controller applies insulin to the virtual patient, which brings the patient back into a normal blood glucose range. It counteracts low blood glucose values with an intake of glucose. The evaluation of the measured CGMS signal using the risk function by Kovatchev and the High Blood Glucose Index (HBGI) and Low Blood Glucose Index (LBGI) (Kovatchev et al. 2005), results in HBGI = {7.60, 5.63, 3.72} and LBGI = {1.71, 0.68, 2.51} for the virtual patients 1, 2 and 3, respectively.
06:00
10:00
18:00
0:00
06:00
10:00
18:00
0:00
06:00
10:00
18:00
0:00
06:00
10:00
18:00
0:00
50 00:00 100 50 00:00
M [g]
implementation does not model conditions like exercise or stress resulting in dropping blood glucose.
06:00
50
00:00
Fig. 4. The simulation set-up.
BG concentration of virtual patient 1
300
00:00
yˆ(t) ui [units]
M [g]
100 50 00:00
time [h]
Fig. 5. Simulation Results for virtual patients 1 and 2. First half of the panels: virtual patient 1; Second half of the panels: virtual patient 2. For each half, the first panel shows the measured CGMS signal y, the second panel the insulin ui suggested by the controller (blue) and an insulin disturbance (red), the third panel the glucose ug suggested by the controller and the fourth panel glucose intake from simulated meals M . proposed controller however aims at explicitly considering the amounts and time points for insulin and glucose administration as optimization variables instead of a discrete approximation of a continuous signal. Other control schemes used previously to determine a protocol for insulin injections around meal times make use of run-to-run or fuzzy control algorithms (Campos-Delegado et al. 2006, Campos-Cornejo et al. 2010, Owens et al. 2006). Mostly those control algorithms require a preset daily plan for meals, which can be limiting for the patient. The controller proposed here however determines the advices of insulin and glucose only from predictions and measurements of the blood glucose concentrations, which gives greater freedom
300
BG concentration of virtual patient 3 300
100 06:00
10:00
18:00
0:00
10:00
18:00
ia [units]
12 : 00 06:00
0:00
50 00:00
06:00
10:00
18:00
0:00
100 50 00:00
06:00
10:00
18:00
0:00
time [h]
y [mg/dL]
300
SMBG
Upper bound
ia [units] M [g]
18 : 00
00 : 00
06 : 00
18 : 00
00 : 00
06 : 00
18 : 00
00 : 00
06 : 00
60
ug [g]
7.5 5 2.5 12 : 00
18 : 00
00 : 00
06 : 00
18 : 00
00 : 00
06 : 00
18 : 00
00 : 00
06 : 00
00 : 00
06 : 00
15 5 12 : 00 60 20 12 : 00
12 : 00 60 20 12 : 00
18 : 00
time [h]
18 : 00
00 : 00
06 : 00
18 : 00
00 : 00
06 : 00
In the case of the real patient data test, the controller seems to give good agreement in addition to the already taken insulin and glucose amounts. However, the performance of the controller depends strongly on the quality of the blood glucose predictions. To determine the doses of insulin and glucose advices, the controller uses blood glucose predictions calculated by a predictor that is not part of the control algorithm. If this predictor delivers unreliable BG predictions to the controller, for example due to bad quality of measured data, the controller will not be able to produce a reliable advice.
60 20 12 : 00
06 : 00
For the virtual patients used as an example here, the control algorithm can bring the blood glucose concentration back into the normal range between 70 mg/dL and 180 mg/dL, both after meals and after an insulin disturbance by suggesting additional glucose and insulin bolus intakes. The virtual patients used for the simulations does not include the effect of, for example, exercise and stress, which both result in dropping blood glucose concentration. Hence, an insulin disturbance was introduced to be able to test the controllers ability to cope with dropping blood glucose.
15
20 12 : 00
00 : 00
to a diabetic patient. Furthermore, an asymmetric cost function is used to determine the MDI advices.
100
5 12 : 00
18 : 00
Lower bound
Lower bound
200
12 : 00
ui [units]
BGTarget YSI
Upper bound
Fig. 8. Simulation Result for real data controller test, day 2. Upper panel: measured BG concentration; 2nd panel: actual insulin taken; 3rd panel: actual meals taken; 4th panel: insulin advice by controller; last panel: carbohydrate advice by the controller
Day 1 CGM
SMBG
7.5 5 2.5
ug [g]
Fig. 6. Simulation Result for virtual patient 3. The first panel shows the measured CGMS signal y, the second panel the insulin ui suggested by the controller (blue) and an insulin disturbance (red), the third panel the glucose ug suggested by the controller and the fourth panel glucose intake from simulated meals M .
BGTarget YSI
100
50 00:00
CGM
200
M [g]
ui [units] ug [g]
y [mg/dL]
200
00:00
ug [g]
day 2
ui [units]
y [mg/dL]
8th IFAC Symposium on Biological and Medical Systems August 29-31, 2012. Budapest, Hungary
time [h]
Fig. 7. Simulation Result for real data controller test, day 1. Upper panel: measured BG concentration; 2nd panel: actual insulin taken; 3rd panel: actual meals taken; 4th panel: insulin advice by controller; last panel: carbohydrate advice by the controller
The proposed controller was tested on three virtual patients and one two-day data set of a real patient. Conclusions are drawn from just a few results. Although this might seem to be limiting generalization, it nonetheless
44
8th IFAC Symposium on Biological and Medical Systems August 29-31, 2012. Budapest, Hungary
shows that the concept of the presented control algorithm is promising. 5. CONCLUSION The proposed controller presented in this paper is based on an algorithm including an optimization problem, which determines the doses and times for insulin and glucose intakes for MDI patients. By using an asymmetric cost function, the optimization minimizes the risk connected to the blood glucose concentration of a patient. In the examples, it can be seen that it is possible to use this control algorithm to produce advices of insulin and glucose intakes for MDI patients, where the advices are boluses and not continuous signals. ACKNOWLEDGEMENTS This research was financially supported by DIAdvisor, an integrated project funded under the European Union’s Seventh Framework Programme (Ref. FP7 IST-216592 DIAdvisor). REFERENCES Campos-Cornejo, F., Campos-Delgado, D., EspinozaTrejo, D., Zisser, H., Jovanovic, L., Doyle, F., and Dassau, E. (2010). An advisory protocol for rapidand slow-acting insulin therapy based on a run-to-run methodology. Diabetes Technology and Therapeutics, 12(v), 555–565. Campos-Delgado, D., Hernandez-Ordonez, M., Femat, R., and Gordillo-Moscoso, A. (2006). Fuzzy-based controller for glucose regulation in type-1 diabetic patients by subcutaneous route. Biomedical Engineering, IEEE Transactions on, 53(11), 2201 –2210. Cescon, M., St˚ ahl, F., and Johansson, R. (2009). Subspacebased model identification of diabetic blood glucose dynamics. In Proc. 15th IFAC Symposium on System Identification (SYSID 2009), July 6 - 8, 2009, SaintMalo, France. Cobelli, C., Dalla Man, C., Sparacino, G., De Nicolao, G., and Kovatchev, B. (2009). Diabetes: Models, signals and control. Biomedical Engineering, IEEE Reviews in, 2, 54–96. Dalla Man, C., Rizza, R., and Cobelli, C. (2007). Meal simulation model of glucose-insulin system. Biomedical Engineering, IEEE Transactions on, 54(10), 1740–1749. De Nicolao, G., Magni, L., Dalla Man, C., and Cobelli, C. (2011). Modeling and control of diabetes: Towards the artificial pancreas. In Proceedings of the 18th IFAC World Congress, volume 18. Milano, Italy. Dexcom (2012). URL http://www.dexcom.com. DIAdvisor (2012). URL http://www.diadvisor.org. Dua, P., Doyle, F., and Pistikopoulos, E. (2009). Multiobjective blood glucose control for type 1 diabetes. In Medical and Biological Engineering and Computing, 343 –352. Guyton, A.C. and Hall, J.E. (2006). Textbook of medical physiology, chapter 78. Elsevier Saunders, 11 edition. Harvey, R., Wang, Y., Grosman, B., Percival, M., Bevier, W., Finan, D., Zisser, H., Seborg, D., Jovanovic, L., and Doyle, F.J.and Dassau, E. (2010). Quest for the artificial
45
pancreas: Combining technology with treatment. Engineering in Medicine and Biology Magazine, IEEE, 29, 53–62. Kirchsteiger, H. and del Re, L. (2009). Reduced hypoglycemia risk in insulin bolus therapy using asymmetric cost functions. In 7th Asian Control Conference (ASCC 2009), 751 –756. Hong Kong, China. Kirchsteiger, H., del Re, L., Renard, E., and Mayrhofer, M. (2009). Robustness properties of optimal insulin bolus administrations for type 1 diabetes. In American Control Conference (ACC 2009), volume 47, 2284 – 2289. St Louis, MO, USA. Kovatchev, B.P., Clarke, W.L., Breton, M., Brayman, K., and McCall, A. (2005). Quantifying temporal glucose variability in diabetes via continuous glucose monitoring: Mathematical methods and clinical applications. Diabetes Technology and Therapeutics, 7(6), 849–862. Mathworks (2011). Matlab optimization toolbox user’s guide. Owens, C., Zisser, H., Jovanovic, L., Srinivasan, B., Bonvin, D., and Doyle, J. (2006). Run-to-run control of blood glucose concentrations for people with type 1 diabetes mellitus. Biomedical Engineering, IEEE Transactions on, 53(6), 996 –1005. Parker, R., Doyle, F.I., and Peppas, N. (2001). The intravenous route to blood glucose control. Engineering in Medicine and Biology Magazine, IEEE, 20, 65–73. Percival, M., Wang, Y., Grosman, B., Dassau, E., Zisser, H., Jovanovic, L., and Doyle, F. (2011). Development of a multi-parametric model predictive control algorithm for insulin delivery in type 1 diabetes mellitus using clinical parameters. Journal of Process Control, 21(3), 391 – 404. Poulsen, J.U., Avogaro, A., Chauchard, F., Cobelli, C., Johansson, R., Nita, L., Pogose, M., del Re, L., Renard, E., Sampath, S., Saudek, F., Skillen, M., and Soendergaard, J. (2010). A diabetes management system empowering patients to reach optimised glucose control: From monitor to advisor. In Engineering in Medicine and Biology Society (EMBC2010), 2010 Annual International Conference of the IEEE, 5270 –5271. Buenos Aires, Argentina. Stahl, F. and Johansson, R. (2009). Diabetes mellitus modeling and short-term prediction based on blood glucose measurements. Mathematical Biosciences, 217(2), 101–117. Trogmann, H., Kirchsteiger, H., Castillo Estrada, G., and del Re, L. (2010a). Fast estimation of meal/insulin bolus effects in T1DM for in silico testing using hybrid approximation of physiological meal/insulin model. Diabetes, 59, A136. Trogmann, H., Kirchsteiger, H., and del Re, H. (2010b). Hybrid control of type 1 diabetes bolus therapy. In Decision and Control (CDC2010), 2010 49th IEEE Conference on, 4721 –4726. Atlanta, GA. YSI (2012). URL http://www.ysi.com.
Diabetic Blood Glucose Control via Optimization over Insulin and Glucose Doses ⋆ Meike Stemmann ∗ Rolf Johansson ∗∗ ∗
Lund University, Sweden (e-mail: [email protected]). Lund University, Sweden (e-mail: [email protected])
∗∗
Abstract: The proposed algorithm aims at helping diabetic patients which use multi-dose injection (MDI) therapy to determine the dose and time of insulin and glucose intakes to stabilize their blood glucose (BG) concentration. The objective is to spend most of the time in normal BG range (70-180 mg/dl). The control was done by minimizing the risk of hypo- and hyperglycemia over the doses and times of insulin and glucose intakes. Keywords: Biomedical Control, Diabetes, Optimization 1. INTRODUCTION
for CSII patients, MDI patients need single insulin advices just a couple of times per day.
Diabetes Mellitus Type I (DMTI) is a chronic disease characterized by the body’s inability to regulate its own blood glucose concentration. In Type 1 diabetes, this is caused by a lack of insulin secretion resulting from a failure of the β-cells in the pancreas, due to either injuries or diseases. Viral infections and autoimmune disorders have both been associated with the destruction of the β-cells (Guyton and Hall 2006). Hence, patients with DMTI have to administer external insulin to substitute the lack of insulin production by the β-cells. This is done either by multiple daily injections (MDI) or through continuous subcutaneous insulin infusion (CSII) using an insulin pump. MDI patients usually take a single dose of basal insulin on daily basis to supply the body with the basic insulin needed and several additional bolus insulin doses at times when the concentration of glucose in the blood is high, e.g., around meal times. A diabetic patient has to determine the correct doses of insulin needed to stabilize the blood glucose (BG) concentration, and has to solve an optimization problem every day.
To fit this regime, the control algorithm proposed in this paper aims at determining the doses and time points for the insulin and also glucose administration by solving a nonlinear optimization problem. While insulin treatment is important in most of the everyday situations, additional glucose intakes can be important to prevent hypoglycemia under special conditions like exercising or stress (Guyton and Hall 2006). Therefore, in addition to insulin advices, also glucose advices are determined by the proposed control algorithm.
To help the patient with this task, blood glucose prediction algorithms (St˚ ahl and Johansson 2009, Cescon et al. 2009) as well as many different control algorithms have been proposed. These control algorithms reach from proportionalintegral-derivative controller (PID), pole placement over adaptive and run-to-run methods to model predictive control (MPC) (Cobelli et al. 2009, Harvey et al. 2010, Parker et al. 2001). Many of the proposed controllers aim at having a continuous insulin signal, that could be used by an insulin pump to continuously infuse exogenous insulin into the diabetic patient.
2. AN OPTIMIZATION-BASED CONTROL ALGORITHM
However, in the scope of the European project DIAdvisor (Poulsen et al. 2010, DIAdvisor 2012) as a blood glucose prediction and treatment advisory system, the patient should be able to use either CSII or MDI treatment. Whereas continuous insulin administration can be useful ⋆ FP7 IST-216592 DIAdvisor
978-3-902823-10-6/12/$20.00 © 2012 IFAC
40
This paper is organized as follows. First, the overall control algorithm is presented in the beginning of Sec. 2. Then the subsections 2.1, 2.2 and 2.3 present the core part of the control algorithm, which is the optimization problem determining the insulin and glucose advices. Section 3.1 shows the simulation of the control algorithm using virtual patients and Sec. 3.2 shows the evaluation of the control algorithm with data measured on real patients.
The proposed control algorithm determines the doses of insulin and glucose together with the times they should be administered to the patient by solving an optimization problem. This optimization is included in an algorithm, which invokes the optimization when the blood glucose concentration leaves a certain range, see Fig. 1. When the blood glucose concentration y falls under 90 [mg/dL] and at least 120 [min] have passed since the last intake, or it falls under 80 [mg/dL] and at least 15 [min] have passed since the last intake, the optimization problem is solved to determine time and dose for an insulin or glucose intake. Similarly, when the BG concentration rises over 130 [mg/dL] and at least 120 [min] have passed since the last intake, the optimization problem is solved to determine the doses. 10.3182/20120829-3-HU-2029.00069
blood glucose concentration y [mg/dL]
8th IFAC Symposium on Biological and Medical Systems August 29-31, 2012. Budapest, Hungary
The parameters n1 to n6 were estimated individually for the patient to be controlled. This was done using nonlinear constrained optimization, solving the optimization problem (2). The measured blood glucose concentration is denoted by y and the desired reference BG concentration to be reached by yref . The parameters n1 and n4 are constrained to ensure the correct gain of the respective model.
130
90
minimize
||yg (t, tg ) + yi (t, ti ) + yref − y||2
80
subject to
n1 > 0 n4 < 0
n1 ...n6
Figure 2 shows yg (t, tg , ug ) and yi (t, ti , ui ) with ug = 10 g and ui = 1 unit at times tg = ti = 100 as inputs for a sample patient. When both insulin and glucose are taken at different times, those two functions are added to form the total change of blood glucose concentration.
15 120 time since last intake [min]
yg [mg/dL]
Fig. 1. Optimization-based control algorithm: The filled areas represent when the optimization problem is solved in the space of time since last intake vs. BG concentration. The blue, vertically striped area for hyperglycemia and the green, horizontally striped area for hypoglycemia.
yi [mg/dL]
The proposed control algorithm solves the optimization to determine insulin or glucose doses when the blood glucose concentration leaves the range between 90 and 130 [mg/dL], in contrast to the previously mentioned normal range being between 70 to 180 [mg/dL]. This is done because the goal is to keep the BG concentration inside the normal range for as long as possible (Kovatchev et al. 2005). This section is organized as follows. First, a model describing the dynamics of a diabetic patient is described in Sec. 2.1. This model relates the glucose or insulin intakes to the change of BG concentration and is needed for the optimization. Next, the cost function used for the optimization problem is presented in Sec. 2.2. The optimization problem to be solved to determine doses and times of insulin and glucose administration is formulated in Sec. 2.3. This optimization problem uses the cost function and the patient model mentioned above.
To describe the change in blood glucose concentration as a response to an intake of insulin and glucose, a nonlinear patient model is used (Trogmann et al. 2010b, Trogmann et al. 2010a). This model takes the doses and times of glucose or insulin intakes as inputs and has the change of blood glucose concentration as an output. Eq. (1) shows this model for insulin and glucose intakes, where yg (t, tg , ug ) is the change in blood glucose concentration as response to a glucose intake of size ug [g] at time tg [min] and yi (t, ti , ui ) the same for an insulin intake with the size ui [IU] at time ti [min]. e−n2 ·(t−tg ) · (t − tg )n3 · ug 1 + e−0.5·(t−tg ) e−n5 ·(t−ti ) · (t − ti )n6 · ui yi (t, ti , ui ) = n4 · 1 + e−0.5·(t−ti )
40
yg (t, 100, 10): Change in BG for 10 g glucose
30 20 10 0 0
2 0 −2 −4 −6 −8 0
500
yi (t, 100, 1):
500
1000
time [min]
1500
2000
Change in BG for 1 unit insulin
1000
time [min]
1500
2000
Fig. 2. Output of the nonlinear model for a sample patient. Upper panel: Change of BG as response to 10 g of Glucose intake at time tg = 100. Lower panel: Change of BG as response to 1 unit of insulin intake at time ti = 100.
2.1 Glucose and Insulin Dynamic Model
yg (t, tg , ug ) = n1 ·
(2)
2.2 Cost Function For a diabetic patient, the risk connected to hypoglycemia is much larger than the risk connected to hyperglycemia. The proposed cost function for the optimization takes this into account through its asymmetric shape, so that too low BG values imply more risk than too high BG values. In this way, the higher risk of low BG compared to high BG is expressed. The advantage of an asymmetic over a quadratic cost function was already presented in (Kirchsteiger and del Re 2009) and (Dua et al. 2009) for example. The asymmetric cost function used here is shown in Eq. (3) (Kovatchev et al. 2005).
(1) R(y) = 10−3 · (y − yref )2 · (c1 + c2 · 30 · exp(−0.05 · y)) (3)
41
8th IFAC Symposium on Biological and Medical Systems August 29-31, 2012. Budapest, Hungary
The BG concentration is denoted by y and the desired reference BG by yref . The parameters c1 and c2 allow to shape the slopes of the asymmetric cost function.
subject to
70 ≤ yP ≤ 450 ui < 20 ug < 80 The reason for the large upper bound is to have a feasible solution of the optimization problem after meals.
Figure 3 shows the cost function for different values of c1 and c2 . The black dashed line shows the default with c1 = 1 and c2 = 1. The effect on the cost function of increasing or decreasing these parameters is shown for some examples. However, not loose the asymmetric shape it should be kept c2 ≥ 0.4 and c1 ≤ 4. 300
The minimization in the optimization problem (5) is done over the doses and times of insulin and glucose as optimization variables using the Optimization Toolbox from Matlab (Mathworks 2011). 2.4 Data
c1 = 0.1, c2 = 1 c1 = 2, c2 = 1
250
The patient data used in Sec. 3.2 below was collected during an in-hospital period at the Centre Hospitalier Universitaire de Montpellier, in the framework of the European FP7-IST research project DIAdvisorT M , see (DIAdvisor 2012) and (Poulsen et al. 2010). The blood glucose concentration was measured with a continuous blood glucose measurement device (CGMS) (Dexcom 2012). Furthermore, for comparison the blood glucose was measured by the Yellow Springs Instruments (YSI) blood glucose analyser (YSI 2012) and through self-monitored of blood glucose (SMBG). Also the insulin doses taken by a patient were recorded, both for patients using an insulin pump and for patients using multi-dose injection (MDI). The example patient taken for evaluation in Sec. 3.2 is an MDI patient.
c1 = 4, c2 = 1 c1 = 1, c2 = 0.1
200
c1 = 1, c2 = 0.4 c1 = 1, c2 = 2
R
150
c1 = 1, c2 = 3 c1 = 1, c2 = 1
100 50 0 0
50
100
150
200
blood glucose concentration y [mg/dL]
250
Fig. 3. The asymmetric cost function for different values of the parameters ci , i = 1, 2.
3. SIMULATION 3.1 Simulation using a Virtual Patient
2.3 Optimization Problem
In order to simulate the controller in a closed-loop for testing and evaluation, an implementation of a virtual patient is used. The virtual patient is a nonlinear physiologically based model (Dalla Man et al. 2007).
The optimization problem needed to be solved for determination of the doses and times for insulin and glucose intake is given in (5). It uses Eq. (3) as the cost function to be minimized. The BG concentration y in the cost function is replaced with the total BG concentration in the future horizon Hp , called yP , using Eq. (1), as shown in Eq. (4). yP (ui , ug , ti , tg , t) = yg (t, tg , ug ) + yi (t, ti , ui ) + yˆ(t) (4) Here, yˆ(t) is the predicted BG concentration for the future horizon Hp assuming no future insulin or glucose intakes. In this way, yP (t) together with the cost function depend on the doses and times for insulin and glucose administration.
The simulation set-up is shown in Fig. 4. The virtual patient takes glucose and insulin as inputs and gives out the measured CGMS signal y(t). All these signals are used by a prediction algorithm to predict the CGMS over a future horizon. The predictor uses a linear state-space model of the patient, here a virtual patient, and a Kalman filter to calculate the predictions (Cescon et al. 2009). The CGMS predictions and the measured CGMS are used by the controller to determine the doses of insulin and glucose to be given to the virtual patient, as described in Sec. 2.
The idea is to minimize the risk connected to a certain amount of measured blood glucose concentration y over the doses and times of glucose and insulin intakes. Furthermore, the optimization problem (5) penalizes too large doses of insulin or glucose in order to keep them as small as possible. The amount of insulin ui and glucose ug to be given as well as the predicted total BG concentration y are constrained. minimize ui ,ug ,ti ,tg
Hp X
The simulation result for three example virtual patients are shown in Fig. 5 and Fig. 6. In the former, the first four panels show the results for virtual patient 1 and the last four panels for virtual patient 2. For each patient, the first panel shows the measured CGMS signal y at the output of the virtual patients and the 180 mg/dL limit of the normal blood glucose range. The insulin ui suggested by the controller is shown in the 2nd panel, while the glucose ug suggested by the controller is shown in the 3rd panels. The glucose intake from simulated meals M is shown the 4th panel. The red impulse in the 2ndpanel around the time 686 is an insulin disturbance, which is not suggested by the controller but acts on the virtual patient. This is introduced to test the controllers ability to deal with dropping blood glucose, since the virtual patient
10−3 (yP (ui , ug , ti , tg ) − yref )2 (5)
i=1 · c1 + c2 · 30 · e−0.05·yP (ui ,ug ,ti ,tg ) +c3 · |ui | + c4 · |ug |
42
8th IFAC Symposium on Biological and Medical Systems August 29-31, 2012. Budapest, Hungary
Predictor
y [mg/dL]
Virtual Patient
y(t)
400
ui (t) ug [g]
ug (t) Controller
200 100
00:00
M [g]
In addition to closed-loop simulations with a virtual patient, the proposed control algorithm was evaluated using data measured on a real patient. A patient-specific model was estimated according to section 2.1 from this data. Since the data is already measured, it is not possible to evaluate the control algorithm in a closed-loop manner. The evaluation should instead consider weather or not the extra insulin or glucose intakes suggested by the controller are reasonable. The CGMS data measured y by the Dexcom device is fed into the controller input. The advices of insulin ui and glucose ug suggested by the controller are then recorded and compared with the measured CGMS data and actual applied insulin ia and glucose intakes M , that the patient took during the in-hospital stay. Results from this comparison for a sample patient are shown in Fig. 7 and Fig. 8 for two different days. On the first day for example (Fig. 7), the BG concentration still increases after the second meal, so that an additional insulin intake could make sense. On the second day, shown in Fig. 8, the BG concentration decreases into hypoglycemic range, so that the suggested extra carbohydrates appear reasonable. 4. DISCUSSION The proposed control algorithm was designed specifically for MDI patients, which require that advices on insulin boluses are given only a few times per day. Although in many approaches, MPC is used to produce a continuous insulin signal (De Nicolao et al. 2011, Percival et al. 2011), it is possible to use the MPC scheme to produce the required insulin advices for MDI patients through approximation of the insulin signal from the MPC (Kirchsteiger an del Re 2009, Kirchsteiger et al. 2009). The
43
18:00
0:00
06:00
10:00
18:00
0:00
06:00
10:00
18:00
0:00
06:00
10:00
18:00
0:00
20
100 50 time [h]
y [mg/dL]
300
BG concentration of virtual patient 2
200 100 00:00
ui [units]
3.2 Open-Loop Simulation with Real Patient Data
10:00
40
00:00
ug [g]
After a meal, the controller applies insulin to the virtual patient, which brings the patient back into a normal blood glucose range. It counteracts low blood glucose values with an intake of glucose. The evaluation of the measured CGMS signal using the risk function by Kovatchev and the High Blood Glucose Index (HBGI) and Low Blood Glucose Index (LBGI) (Kovatchev et al. 2005), results in HBGI = {7.60, 5.63, 3.72} and LBGI = {1.71, 0.68, 2.51} for the virtual patients 1, 2 and 3, respectively.
06:00
10:00
18:00
0:00
06:00
10:00
18:00
0:00
06:00
10:00
18:00
0:00
06:00
10:00
18:00
0:00
50 00:00 100 50 00:00
M [g]
implementation does not model conditions like exercise or stress resulting in dropping blood glucose.
06:00
50
00:00
Fig. 4. The simulation set-up.
BG concentration of virtual patient 1
300
00:00
yˆ(t) ui [units]
M [g]
100 50 00:00
time [h]
Fig. 5. Simulation Results for virtual patients 1 and 2. First half of the panels: virtual patient 1; Second half of the panels: virtual patient 2. For each half, the first panel shows the measured CGMS signal y, the second panel the insulin ui suggested by the controller (blue) and an insulin disturbance (red), the third panel the glucose ug suggested by the controller and the fourth panel glucose intake from simulated meals M . proposed controller however aims at explicitly considering the amounts and time points for insulin and glucose administration as optimization variables instead of a discrete approximation of a continuous signal. Other control schemes used previously to determine a protocol for insulin injections around meal times make use of run-to-run or fuzzy control algorithms (Campos-Delegado et al. 2006, Campos-Cornejo et al. 2010, Owens et al. 2006). Mostly those control algorithms require a preset daily plan for meals, which can be limiting for the patient. The controller proposed here however determines the advices of insulin and glucose only from predictions and measurements of the blood glucose concentrations, which gives greater freedom
300
BG concentration of virtual patient 3 300
100 06:00
10:00
18:00
0:00
10:00
18:00
ia [units]
12 : 00 06:00
0:00
50 00:00
06:00
10:00
18:00
0:00
100 50 00:00
06:00
10:00
18:00
0:00
time [h]
y [mg/dL]
300
SMBG
Upper bound
ia [units] M [g]
18 : 00
00 : 00
06 : 00
18 : 00
00 : 00
06 : 00
18 : 00
00 : 00
06 : 00
60
ug [g]
7.5 5 2.5 12 : 00
18 : 00
00 : 00
06 : 00
18 : 00
00 : 00
06 : 00
18 : 00
00 : 00
06 : 00
00 : 00
06 : 00
15 5 12 : 00 60 20 12 : 00
12 : 00 60 20 12 : 00
18 : 00
time [h]
18 : 00
00 : 00
06 : 00
18 : 00
00 : 00
06 : 00
In the case of the real patient data test, the controller seems to give good agreement in addition to the already taken insulin and glucose amounts. However, the performance of the controller depends strongly on the quality of the blood glucose predictions. To determine the doses of insulin and glucose advices, the controller uses blood glucose predictions calculated by a predictor that is not part of the control algorithm. If this predictor delivers unreliable BG predictions to the controller, for example due to bad quality of measured data, the controller will not be able to produce a reliable advice.
60 20 12 : 00
06 : 00
For the virtual patients used as an example here, the control algorithm can bring the blood glucose concentration back into the normal range between 70 mg/dL and 180 mg/dL, both after meals and after an insulin disturbance by suggesting additional glucose and insulin bolus intakes. The virtual patients used for the simulations does not include the effect of, for example, exercise and stress, which both result in dropping blood glucose concentration. Hence, an insulin disturbance was introduced to be able to test the controllers ability to cope with dropping blood glucose.
15
20 12 : 00
00 : 00
to a diabetic patient. Furthermore, an asymmetric cost function is used to determine the MDI advices.
100
5 12 : 00
18 : 00
Lower bound
Lower bound
200
12 : 00
ui [units]
BGTarget YSI
Upper bound
Fig. 8. Simulation Result for real data controller test, day 2. Upper panel: measured BG concentration; 2nd panel: actual insulin taken; 3rd panel: actual meals taken; 4th panel: insulin advice by controller; last panel: carbohydrate advice by the controller
Day 1 CGM
SMBG
7.5 5 2.5
ug [g]
Fig. 6. Simulation Result for virtual patient 3. The first panel shows the measured CGMS signal y, the second panel the insulin ui suggested by the controller (blue) and an insulin disturbance (red), the third panel the glucose ug suggested by the controller and the fourth panel glucose intake from simulated meals M .
BGTarget YSI
100
50 00:00
CGM
200
M [g]
ui [units] ug [g]
y [mg/dL]
200
00:00
ug [g]
day 2
ui [units]
y [mg/dL]
8th IFAC Symposium on Biological and Medical Systems August 29-31, 2012. Budapest, Hungary
time [h]
Fig. 7. Simulation Result for real data controller test, day 1. Upper panel: measured BG concentration; 2nd panel: actual insulin taken; 3rd panel: actual meals taken; 4th panel: insulin advice by controller; last panel: carbohydrate advice by the controller
The proposed controller was tested on three virtual patients and one two-day data set of a real patient. Conclusions are drawn from just a few results. Although this might seem to be limiting generalization, it nonetheless
44
8th IFAC Symposium on Biological and Medical Systems August 29-31, 2012. Budapest, Hungary
shows that the concept of the presented control algorithm is promising. 5. CONCLUSION The proposed controller presented in this paper is based on an algorithm including an optimization problem, which determines the doses and times for insulin and glucose intakes for MDI patients. By using an asymmetric cost function, the optimization minimizes the risk connected to the blood glucose concentration of a patient. In the examples, it can be seen that it is possible to use this control algorithm to produce advices of insulin and glucose intakes for MDI patients, where the advices are boluses and not continuous signals. ACKNOWLEDGEMENTS This research was financially supported by DIAdvisor, an integrated project funded under the European Union’s Seventh Framework Programme (Ref. FP7 IST-216592 DIAdvisor). REFERENCES Campos-Cornejo, F., Campos-Delgado, D., EspinozaTrejo, D., Zisser, H., Jovanovic, L., Doyle, F., and Dassau, E. (2010). An advisory protocol for rapidand slow-acting insulin therapy based on a run-to-run methodology. Diabetes Technology and Therapeutics, 12(v), 555–565. Campos-Delgado, D., Hernandez-Ordonez, M., Femat, R., and Gordillo-Moscoso, A. (2006). Fuzzy-based controller for glucose regulation in type-1 diabetic patients by subcutaneous route. Biomedical Engineering, IEEE Transactions on, 53(11), 2201 –2210. Cescon, M., St˚ ahl, F., and Johansson, R. (2009). Subspacebased model identification of diabetic blood glucose dynamics. In Proc. 15th IFAC Symposium on System Identification (SYSID 2009), July 6 - 8, 2009, SaintMalo, France. Cobelli, C., Dalla Man, C., Sparacino, G., De Nicolao, G., and Kovatchev, B. (2009). Diabetes: Models, signals and control. Biomedical Engineering, IEEE Reviews in, 2, 54–96. Dalla Man, C., Rizza, R., and Cobelli, C. (2007). Meal simulation model of glucose-insulin system. Biomedical Engineering, IEEE Transactions on, 54(10), 1740–1749. De Nicolao, G., Magni, L., Dalla Man, C., and Cobelli, C. (2011). Modeling and control of diabetes: Towards the artificial pancreas. In Proceedings of the 18th IFAC World Congress, volume 18. Milano, Italy. Dexcom (2012). URL http://www.dexcom.com. DIAdvisor (2012). URL http://www.diadvisor.org. Dua, P., Doyle, F., and Pistikopoulos, E. (2009). Multiobjective blood glucose control for type 1 diabetes. In Medical and Biological Engineering and Computing, 343 –352. Guyton, A.C. and Hall, J.E. (2006). Textbook of medical physiology, chapter 78. Elsevier Saunders, 11 edition. Harvey, R., Wang, Y., Grosman, B., Percival, M., Bevier, W., Finan, D., Zisser, H., Seborg, D., Jovanovic, L., and Doyle, F.J.and Dassau, E. (2010). Quest for the artificial
45
pancreas: Combining technology with treatment. Engineering in Medicine and Biology Magazine, IEEE, 29, 53–62. Kirchsteiger, H. and del Re, L. (2009). Reduced hypoglycemia risk in insulin bolus therapy using asymmetric cost functions. In 7th Asian Control Conference (ASCC 2009), 751 –756. Hong Kong, China. Kirchsteiger, H., del Re, L., Renard, E., and Mayrhofer, M. (2009). Robustness properties of optimal insulin bolus administrations for type 1 diabetes. In American Control Conference (ACC 2009), volume 47, 2284 – 2289. St Louis, MO, USA. Kovatchev, B.P., Clarke, W.L., Breton, M., Brayman, K., and McCall, A. (2005). Quantifying temporal glucose variability in diabetes via continuous glucose monitoring: Mathematical methods and clinical applications. Diabetes Technology and Therapeutics, 7(6), 849–862. Mathworks (2011). Matlab optimization toolbox user’s guide. Owens, C., Zisser, H., Jovanovic, L., Srinivasan, B., Bonvin, D., and Doyle, J. (2006). Run-to-run control of blood glucose concentrations for people with type 1 diabetes mellitus. Biomedical Engineering, IEEE Transactions on, 53(6), 996 –1005. Parker, R., Doyle, F.I., and Peppas, N. (2001). The intravenous route to blood glucose control. Engineering in Medicine and Biology Magazine, IEEE, 20, 65–73. Percival, M., Wang, Y., Grosman, B., Dassau, E., Zisser, H., Jovanovic, L., and Doyle, F. (2011). Development of a multi-parametric model predictive control algorithm for insulin delivery in type 1 diabetes mellitus using clinical parameters. Journal of Process Control, 21(3), 391 – 404. Poulsen, J.U., Avogaro, A., Chauchard, F., Cobelli, C., Johansson, R., Nita, L., Pogose, M., del Re, L., Renard, E., Sampath, S., Saudek, F., Skillen, M., and Soendergaard, J. (2010). A diabetes management system empowering patients to reach optimised glucose control: From monitor to advisor. In Engineering in Medicine and Biology Society (EMBC2010), 2010 Annual International Conference of the IEEE, 5270 –5271. Buenos Aires, Argentina. Stahl, F. and Johansson, R. (2009). Diabetes mellitus modeling and short-term prediction based on blood glucose measurements. Mathematical Biosciences, 217(2), 101–117. Trogmann, H., Kirchsteiger, H., Castillo Estrada, G., and del Re, L. (2010a). Fast estimation of meal/insulin bolus effects in T1DM for in silico testing using hybrid approximation of physiological meal/insulin model. Diabetes, 59, A136. Trogmann, H., Kirchsteiger, H., and del Re, H. (2010b). Hybrid control of type 1 diabetes bolus therapy. In Decision and Control (CDC2010), 2010 49th IEEE Conference on, 4721 –4726. Atlanta, GA. YSI (2012). URL http://www.ysi.com.