The glucose-insulin model is defined [11Chase JG, Shaw GM, Lin J, et al. Adaptive bolus-based targeted glucose regulation of hyperglycaemia in critical care Med Eng Phys 2005; 27(1): 1-11.-13Wong XW, Singh-Levett I, Hollingsworth LJ, et al. A novel, model-based insulin and nutrition delivery controller for glycemic regulation in critically ill patients Diabetes Technol Ther 2006; 8(2): 174-90., 23Hann CE, Chase JG, Lin J, et al. Integral-based parameter identification for long-term dynamic verification of a glucose-insulin system model Comput Methods Programs Biomed 2005; 77(3): 259-70.]:

where G(t) is the plasma glucose concentration (mmol/L); G_E the equilibrium level of plasma glucose concentration (mmol/L); Q(t) the interstitial insulin; I(t) the concentration of the plasma insulin above basal level (mU/L); P(t) the exogenous glucose infusion rate (mmol/(L min)); u(t) the insulin infusion rate (mU/min); V the assumed insulin distribution volume (L); n the delay in interstitial transfer of insulin (min^-1); p_G the fractional clearance of plasma glucose at basal insulin (min^-1); S_I the time-varying insulin sensitivity (L/mU min); k the parameter controlling the effective half life of insulin (min^-1); and α_G the Michaelis-Menten parameter for glucose clearance saturation. For more details on the construction and physiological interpretation of the model Equations (1)-(3) see [11Chase JG, Shaw GM, Lin J, et al. Adaptive bolus-based targeted glucose regulation of hyperglycaemia in critical care Med Eng Phys 2005; 27(1): 1-11.-13Wong XW, Singh-Levett I, Hollingsworth LJ, et al. A novel, model-based insulin and nutrition delivery controller for glycemic regulation in critically ill patients Diabetes Technol Ther 2006; 8(2): 174-90., 23Hann CE, Chase JG, Lin J, et al. Integral-based parameter identification for long-term dynamic verification of a glucose-insulin system model Comput Methods Programs Biomed 2005; 77(3): 259-70.].

For the control of blood glucose G(t) in Equation (1), measurements are assumed to be taken every hour with a normally distributed absolute error of 7%, which is typical for a commercial glucometer [42Arkray GlucocardTM Test Strip 2 Data Sheet Japan Arkray Inc 2001.]. Model-based control of glucose typically starts by taking two measurements G₀ and G₆₀ at the times 0 and 60 minutes and computing the insulin sensitivity S_I from Equation (10). The goal is to determine the required insulin infusion u_I or bolus u_B in Equations (1) and (6) that will bring glucose down to a target value G_target in the next hour.

Let S_I,1 be the solution of Equation (10) that determines the insulin sensitivity in the first hour. Define S_I,2 as the insulin sensitivity in the second hour. In the ICU a patient’s condition can change rapidly as a result of a disease state or drug therapy. Therefore, S_I can change significantly over time [2Chase JG, Shaw GM, Wong XW, et al. Model-based glycaemic control in critical care - A review of the state of the possible Biomed Signal Process Control [review and tutorial] 2006; 1(1): 3-21., 23Hann CE, Chase JG, Lin J, et al. Integral-based parameter identification for long-term dynamic verification of a glucose-insulin system model Comput Methods Programs Biomed 2005; 77(3): 259-70.]. Given S_I,1 in the first hour, it is thus possible that S_I,2 in the second hour may have changed. An approximation to the insulin sensitivity S_I,2 in the second hour for predicting potential outcomes of an intervention at the end of hour 1, is defined S¯I,2=SI,1 . As long as the true SI,2 doesn’t change significantly from S_I,1, this value S¯I,2 can be used to determine the insulin control input u(t) in Equation (1) that will bring the glucose to the target value of G_target. Any significant changes will induce increasingly, unavoidable errors in the prediction.

First assume that u_B = 0 in Equation (6) and that only constant insulin infusion in the second hour is used. An example is given in Fig. (1), which includes a “true” glucose response to an infusion of u_I = 2 units over 1 hour with a nutritional input of P(t) = 0.03mmol/L/min. Insulin sensitivity S_I,1 = 0.0008 (L/mU min), S_I,2 = 0.001 (L/mU min), G_E=4.5 mmol and the rest of the parameters in Equation (1) are given in Equation (5). The “measurement” points G₀=8 and G₆₀=7.26 are denoted by crosses (+) in Fig. (1) and the target glucose G_target=5 mmol/L is denoted by a circle (o). No noise is added.

Fig. (1) Controlling glucose to a target value of Gtarget=5 mmol/L. The true SI in the first and second hours are defined as SI,1= 0.0008 and SI,2= 0.001 (L/mU min).

For simplicity, it is assumed that S_I is precisely known in the first hour. In practice, either the solution to Equation (10) or Equation (12) would approximate S_I . Assuming that S¯I,2=0.0008 in the second hour, the goal is to find the infusion u_I such that the numerical solution G(t) to Equation (1) with S_I=S_I,2, and initial conditions, {G(0)=G₆₀, Q(0)=Q₆₀, I(0)=I₆₀} satisfies G(60)=G_target. The values of Q₆₀ and I₆₀ are determined by the evaluating the numerical solution to Equations (2)-(3) at t = 60. Note that without loss of generality, the time at the beginning of drug intervention is assumed to be at 0 minutes and the target value is assumed to be at 60 minutes.

To determine u_I, Equation (1) is solved numerically for a range of infusion values u_I, and the resulting end glucose value is compared to the target value. The end glucose value is represented as a function G¯targetuI , and is defined:

The correct u_I is denoted u_I,target and is defined:

where G¯targetuI is defined in Equation (13). Define the points:

where u_I is treated as a variable on the y axis. Fig. (2) shows the points G¯target,1,uI,1,...,G¯target,7,uI,7 plotted as crosses (+). A cubic spline is then fitted to the data, which is shown as the solid line in Fig. (2). Evaluating the cubic spline at G_target in Fig. (2) allows a good approximation to u_I,target of Equation (14) with only 4 numerical solutions of Equations (1)-(3) required. This overall method for determining u_I,target of Equation (14) is summarized in Fig. (3).

Fig. (2) Plotting the points G¯target,1,uI,1,...,G¯target,7,uI,7 from Equation (16) and fitting a cubic spline to determine uI,target in Equation (14).

For this example, the target infusion was calculated to be u_I = u_I,target = 3.27U. The resulting glucose response with S_I held constant at the approximate value of 0.0008 L/mU min is denoted by a dashed line in Fig. (1). This approximate glucose response hits the target G_target as required. However, the “true” glucose response, which comes from using the output infusion u_I,target in Fig. (2) with the true value of S_I,2 = 0.001, slightly undershoots G_target in Fig. (1). The end result in this case is still accurate, with an error of 8%. In practice both noise, modelling error and natural variation in S_I can effect the accuracy of hitting the target glucose and is investigated in detail in the results.

The most common approach to parameter identification as discussed in the introduction are methods that rely on many forward simulations. A standard non-linear regression least squares (NRLS) gradient descent algorithm was tested rigorously in [23Hann CE, Chase JG, Lin J, et al. Integral-based parameter identification for long-term dynamic verification of a glucose-insulin system model Comput Methods Programs Biomed 2005; 77(3): 259-70.]. Assuming a reasonable starting guess, the NRLS method was thousands of times slower than the integral method of [23Hann CE, Chase JG, Lin J, et al. Integral-based parameter identification for long-term dynamic verification of a glucose-insulin system model Comput Methods Programs Biomed 2005; 77(3): 259-70.]. Furthermore, local minima’s were often found so that the best insulin sensitivity estimate S_I was not always found.

The problem of local minima’s in NRLS can always be corrected by starting at many starting points, like the method of Cobelli [19Cobelli C, Caumo A, Omenetto M. Minimal model SG overestimation and SI underestimation: improved accuracy by a Bayesian two-compartment model Am J Physiol 1999; 277(3 Pt 1): E481-8.]. However, this dramatically increases the number of forward simulations. For example in [23Hann CE, Chase JG, Lin J, et al. Integral-based parameter identification for long-term dynamic verification of a glucose-insulin system model Comput Methods Programs Biomed 2005; 77(3): 259-70.], the integral method was 1000 times fast than the NRLS algorithm which started from one initial guess. If 10-100 starting points were used for the NRLS algorithm, which is quite typical to ensure accuracy (e.g. Cobelli [19Cobelli C, Caumo A, Omenetto M. Minimal model SG overestimation and SI underestimation: improved accuracy by a Bayesian two-compartment model Am J Physiol 1999; 277(3 Pt 1): E481-8.]), the integral method would be 10000-100000 times faster than NRLS. The speed gain increases even further as the complexity of the model and number of fitted parameters increase, for example a cardiovascular model (e.g. [24Hann CE, Chase JG, Shaw GM. Integral-based identification of patient specific parameters for a minimal cardiac model Comput Methods Programs Biomed 2006; 81(2): 181-92., 26Starfinger C, Hann CE, Chase JG, et al. Model-based cardiac diagnosis of pulmonary embolism Comput Methods Programs Biomed 2007; 87(1): 46-60.]).

For a model-based therapeutics approach [2Chase JG, Shaw GM, Wong XW, et al. Model-based glycaemic control in critical care - A review of the state of the possible Biomed Signal Process Control [review and tutorial] 2006; 1(1): 3-21., 9Lonergan T, Compte AL, Willacy M, et al. A pilot study of the SPRINT protocol for tight glycemic control in critically Ill patients Diabetes Technol Ther 2006; 8(4): 449-62.-14Chase JG, Shaw GM, Lotz T, et al. Model-based insulin and nutrition administration for tight glycaemic control in critical care Curr Drug Deliv 2007; 4(4): 283-96.], the large number of forward simulations required in the NRLS approach is extremely costly, and is not feasible to implement. Note that an NRLS approach could be applied in the model-based control examples of this paper, and would give similar results to the integral method, but it comes at a considerable computational cost. Therefore, since the model-based therapeutics approach requires minimal computation, this paper focuses entirely on methods that do not require a forward simulation. The two methods considered are the derivative method of Equations (11) and (12) and the integral method of Equations (7)-(10).

The derivative approach is a commonly used concept, for example in gradient descent algorithms, and would therefore most likely be the more easily understood and derived method. It is also perhaps, a more natural way of proceeding, since the original differential equation model is written in terms of derivatives. Therefore, the derivative approach at first sight would appear to be the simplest to implement and potentially a reasonable way of avoiding forward simulations in the parameter identification part of the model-based control algorithm of Fig. (3).

Fig. (3) An algorithm summarizing the method of model-based glucose control which determines the required insulin infusion that brings the blood glucose to a predetermined glucose target Gtarget. Similar approaches can be used in appropriate time frames or intervals for any drug therapy that is similarly modelled with differential equations.

However, as is shown in the results, the integral formulation, which in general is perhaps a less known and accepted way of representing a differential equation model; is in fact fundamental for reliable results. This phenomenon was also investigated in a related approach to parameter identification of a minimal cardiac model on clinical pulmonary embolism animal data [26Starfinger C, Hann CE, Chase JG, et al. Model-based cardiac diagnosis of pulmonary embolism Comput Methods Programs Biomed 2007; 87(1): 46-60.], where even with perfectly smooth, model generated signals, a derivative approach went unstable. The integral approach on the other hand remained stable.

Hence, the main aim of this paper, is to investigate the effect of the two fast parameter identification integral and derivative based methods on the glucose-insulin model; and to better explain the importance of the integral in the formulation. Most importantly, this study is done in the context of model-based therapeutics and glucose control in the Christchurch Hospital Intensive Care Unit.

The glucose control protocol SPRINT [9Lonergan T, Compte AL, Willacy M, et al. A pilot study of the SPRINT protocol for tight glycemic control in critically Ill patients Diabetes Technol Ther 2006; 8(4): 449-62., 10Lonergan T, LeCompte A, Willacy M, et al. A simple insulin-nutrition protocol for tight glycemic control in critical illness: Development and protocol comparison Diabetes Technol Ther 2006; 8(2): 191-206., 14Chase JG, Shaw GM, Lotz T, et al. Model-based insulin and nutrition administration for tight glycaemic control in critical care Curr Drug Deliv 2007; 4(4): 283-96., 32Shaw GM, Chase JG, Wong J, et al. Rethinking glycaemic control in critical illness - from concept to clinical practice change Crit Care Resusc 2006; 8(2): 90-.] is now used extensively in the Christchurch ICU. One of the keys to the success of SPRINT is the significant testing of model-based glucose control algorithms on “virtual” patients prior to implementation. The major physiological variable that is used to represent a “virtual” patient profile is the time varying insulin sensitivity S_I = S_I (t) profile in Equation (1) that can be identified from retrospective data.

The integral-based parameter identification method [23Hann CE, Chase JG, Lin J, et al. Integral-based parameter identification for long-term dynamic verification of a glucose-insulin system model Comput Methods Programs Biomed 2005; 77(3): 259-70.] allowed fast and accurate insulin sensitivity profiles to be constructed for long term patient data. These profiles allowed an accurate physiological representation of a patient’s metabolic dynamics over periods of up to 1-2 weeks [23Hann CE, Chase JG, Lin J, et al. Integral-based parameter identification for long-term dynamic verification of a glucose-insulin system model Comput Methods Programs Biomed 2005; 77(3): 259-70.], and were a fundamental element in the development of SPRINT [9Lonergan T, Compte AL, Willacy M, et al. A pilot study of the SPRINT protocol for tight glycemic control in critically Ill patients Diabetes Technol Ther 2006; 8(4): 449-62., 10Lonergan T, LeCompte A, Willacy M, et al. A simple insulin-nutrition protocol for tight glycemic control in critical illness: Development and protocol comparison Diabetes Technol Ther 2006; 8(2): 191-206., 14Chase JG, Shaw GM, Lotz T, et al. Model-based insulin and nutrition administration for tight glycaemic control in critical care Curr Drug Deliv 2007; 4(4): 283-96., 32Shaw GM, Chase JG, Wong J, et al. Rethinking glycaemic control in critical illness - from concept to clinical practice change Crit Care Resusc 2006; 8(2): 90-.].

The insulin sensitivity profiles provide a means to simulate physiologically realistic time varying glucose response to different insulin and nutrition regimes. This approach thus provides a repeatable cohort for easy comparison of various protocols. It also gives insight into long term clinical performance, and, importantly, lets algorithms and methods be tested safely before clinical implementation.

Fig. (4) shows a comparison of the “virtual clinical trials” versus the clinical data from the SPRINT trial in the Christchurch ICU for the first 16,000 clinical measurements and 24,000 hours of control. The distributions for the “virtual trials” are very close to both the raw SPRINT data and a lognormal fit of the data. The results of the virtual patient trials of other protocols [43Laver S, Preston S, Turner D, McKinstry C, Padkin A. Implementing intensive insulin therapy: development and audit of the Bath insulin protocol Anaesth Intensive Care 2004; 32(3): 311-6.-45Goldberg PA, Siegel MD, Sherwin RS, et al. Implementation of a safe and effective insulin infusion protocol in a medical intensive care unit Diabetes Care 2004; 27(2): 461-7.] (not shown) also match their reports. The tightness of the SPRINT results and good correlation of other protocols serves as a significant validation of the methods and approach.

Fig. (4) A comparison of the virtual trials approach and real clinical ICU results from SPRINT. Also shown are virtual patient simulations of two other well known protocols, Van den Berghe [46Van den Berghe G, Wouters P, Weekers F, et al. Intensive insulin therapy in critically ill patients N Engl J Med 2001; 345(19): 1359-67.] and Krinsley [47Krinsley JS. Decreased mortality of critically ill patients with the use of an intensive glycemic management protocol Crit Care Med 2003; 31: A19.].

To further illustrate the impact of SPRINT, Fig. (5) shows a patient on SPRINT compared to a patient on a previously implemented clinical sliding scale in the Christchurch ICU. The measurements in both Fig. (5a) and (5b) are taken every hour, but the SPRINT patient is significantly better controlled than the patient on the original standard sliding scale. Note that the SPRINT patient also has 2-4 potentially contaminated measurements but still provides better control to a 4-6.1 mmol/L or similar target band than the retrospective data patient who is less acutely ill by APACHE II score.

Fig. (5) (a) A patient (Patient 130) on a typical sliding scale before the use of SPRINT in the Christchurch ICU. APACHE II score = 11. (b) A patient (Patient 5005) on SPRINT in the Christchurch ICU. APACHE II score = 21.

For ease of reference in this section and the following sections, Equations (7)-(10) are referred to as the Integral Method and Equations (11)-(12) are referred to as the Derivative Method. The methods are derived from the same set of differential Equations (1)-(3). Therefore, if no noise is present, it may be reasonable to suggest that they should perform equally well when identifying S_I. In addition, neither requires the forward simulation used in most typical identification approaches. To test this assumption, the following set of parameters is considered:

“Measured” glucose values at t = 0 and t = 60 are generated by numerically solving Equations (1)-(3) for various values of G_E, and initial conditions Q₀ with a fixed initial glucose of G₀ = 5 mmol/L. Two main parameters sets are considered:

In Equation (18), Q₀ is varied in steps of 1 mU/L and in Equation (19), G_E is varied in steps of 1 mmol/L. Figs. (6) and (7) show the results of the identified S_I using the integral and derivative methods for each parameter set of Equations (18) and (19).

Fig. (6) The identified insulin sensitivity SI for the integral method of Equations (7)-(10) and derivative method of Equations (11)-(12) for the parameter set of Equation (18).

Fig. (7) The identified insulin sensitivity SI for the integral method of Equations (7)-(10) and derivative method of Equations (11)-(12) for the parameter set of Equation (19).

Fig. (6) shows that for Q₀ < 5 mU/L the derivative method gives an S_I value that is significantly different from the true value. In fact it becomes non-physiological and negative for Q₀ = 0 and 1. The integral method, in contrast, remains stable. Fig. (7) shows a similar result with the derivative method rapidly diverging after G_E = 4 mmol/L, and the integral method staying virtually constant.

The scenarios of Figs. (6) and (7) can be realized in practice whenever the insulin is cut off, so that Q(t) reaches low levels, followed by an increase of carbohydrate input (see example to follow). In particular, a negative S_I would occur whenever the true S_I is sufficiently low, so that the typical undershooting that occurs with the derivative method goes less than zero, as the example in Fig. (6) demonstrates.

3.2.1. Model-based Glucose Control Example – Minimizing Insulin Infusion

To demonstrate the results of Figs. (6) and (7) in a clinical setting, a patient from the retrospective cohort of [23Hann CE, Chase JG, Lin J, et al. Integral-based parameter identification for long-term dynamic verification of a glucose-insulin system model Comput Methods Programs Biomed 2005; 77(3): 259-70.] is considered. The patient used is Patient 554, who was a female aged 20; type 1 diabetic; medical subgroup – Other Medical; APACHE II score - 26. Seven hours of data is analyzed, and Fig. (8) shows the time-varying insulin sensitivity for this period taken from [23Hann CE, Chase JG, Lin J, et al. Integral-based parameter identification for long-term dynamic verification of a glucose-insulin system model Comput Methods Programs Biomed 2005; 77(3): 259-70.]. Patient 554 also has the parameters:

Fig. (8) Time varying insulin sensitivity for Patient 554 from the retrospective cohort [23Hann CE, Chase JG, Lin J, et al. Integral-based parameter identification for long-term dynamic verification of a glucose-insulin system model Comput Methods Programs Biomed 2005; 77(3): 259-70.].

(20)  GE=4.5 mmol/L and G0=5.4 mmol/L

and all the other parameters are defined in Equation (5).

To begin the model-based control algorithm of Fig. (3), two glucose values are required in the first hour. These values are generated by solving Equations (1)-(4) with the parameters of Equation (5) and (20); an insulin infusion input of u_I = 0.5 U , u_B = 0 for u(t) in Equation (6); and initial conditions for insulin defined, by I₀ = 1 mU/min, Q₀ = 1 mU/min. The target glucose in Step 4 of Fig. (3) is G_target = 5mmol/L. Additional constraints for this example, are that the use of exogenous insulin u_I is minimized and is only in steps of 0.5 U, and that when possible, the carbohydrate input P(t) , is the primary controller with a resolution of 0.01 mmol/L/min. Finally, it is assumed that for hour 6 the feed is increased to 0.06 mmol/L/min.

The identified insulin sensitivity for each of the derivative and integral methods is shown in Fig. (9), along with the true insulin sensitivity of Fig. (8). Notice that even without noise, both parameter identification methods deteriorate at hours 5 and 6, but the integral method is the most accurate. The absolute percentage errors of the methods for hours 1-6 in Fig. (9) are:

Fig. (9) The identified insulin sensitivity values for the derivative and integral methods compared to the true insulin sensitivity for Patient 554.

(21)  errorderivative=1.3,10.5,5.6,0.4,34,55.2%errorintegral=3.5,6.9,3.0,0.8,19.7,13.0%

This deterioration is a result of low insulin levels which progressively removes the effect of S_I on the glucose response, and thus the large errors in S_I have a negligible effect on glucose control, which is shown in hours 1-6 of Fig. (10). This state of no insulin and very little carbohydrate input, of course could not be sustained for any significant period of time, as the patient would face malnutrition/starvation. Thus, the feed is increased at hour 6.

Fig. (10) Model-based glucose control using the algorithm of Fig. (3), with the added constraint of minimizing the exogeneous insulin.

There is no reliable insulin sensitivity value from the prior hour due to the very low insulin levels. Therefore, a conservative infusion of 0.5 mmol/L/min is applied at hour 6 to identify insulin sensitivity so that the algorithm of Fig. (3) can be applied accurately in the following hour. Fig. (9) and Equation (21) show that the integral method identifies the insulin sensitivity quite accurately at hour 6 with an error of 13.0%, where the derivative method has a much larger error of 55.2%.

However, the significant under prediction of insulin sensitivity for the derivative method dramatically affects control. Fig. (10) shows that control in hour 7 for the derivative method is poor, with a 5.5 mU bolus predicted and an undesirable, and potentially dangerous, hypoglycaemia event of 3.61 mmol/L. This result corresponds to an error of 27.8% in the target glucose. On the other hand, the control based on the integral method is good with a final glucose value of 4.83 mmol/L, which corresponds to a 3.4% error. The results of Figs. (9) and (10) further confirm the observations of Figs. (6) and (7).

Note that the scenario of Fig. (10) is not uncommon in critical care and for the extended retrospective data set given in Shaw et al. [48Shaw GM, Chase JG, Lee DS, et al. Peak and range of blood glucose are also associated with ICU mortality Critical Care Med 2005; 32(12): A125.], can occur several times daily. Therefore, the integral method allows more flexibility in the control protocol by not requiring insulin infusion to be on constantly, and is robust to sudden increases in the carbohydrate input.

3.2.2. Model-Based Glucose Control – Constant S_I Approximation

To further test the methods for a longer period and to demonstrate the practical, clinical issues associated with model-based control, another patient from the retrospective data of [23Hann CE, Chase JG, Lin J, et al. Integral-based parameter identification for long-term dynamic verification of a glucose-insulin system model Comput Methods Programs Biomed 2005; 77(3): 259-70.] is used. The patient is Patient 519, who was a male aged 69; type 2 diabetic; medical subgroup - General Surgical; APACHE II score - 29. The integral and derivative methods are compared based on a constant S_I, which is taken to be the mean S_I of patient 519. Similarly, the parameters P(t) and G_E in Equations (1)-(4) are set constant at the mean nutritional input and mean glucose respectively of patient 519. The numerical values of the parameters are thus defined:

(22)  SI=9.28×10−4 L/mU min, Pt=0.049 mmol/Lmin,GE=5.84mmol/L

The rest of the model parameters are defined in Equation (5). Data for the first 3 days of patient 519 is used to test the predictive model-based glucose control of Fig. (3). Note that the protocol of minimizing the insulin, which was implemented in Fig. (9), is not used.

To begin the model-based control algorithm of Fig. (3), two glucose values are required in the first hour. These values are generated by solving Equations (1)-(4) with the parameters of Equation (22); an insulin infusion input of u_I = 1U , u_B = 0 for u(t) in Equation (6); and initial conditions of G₀ = 11.5 mmol/L, I₀ = 0 mU/min, Q₀ = 0 mU/min. The target glucose in Step 4 of Fig. (3) is defined:

(23)  Gtarget=max{G0−1,5}

where for each consecutive hour the initial conditions G₀, I₀ and Q₀ are taken as the previously calculated G₆₀, Q₆₀ and I₆₀, as detailed in Step 2 of Fig. (3). Equation (23) ensures the reductions in glucose are not too large which clinically, may be undesirable for the patient.

Every hour that the algorithm of Fig. (3) is applied, a new infusion u_I,target is defined for the next hour, which in turn defines a new glucose response, and so on as long as required. In this example, the final time is at 3 days or 72 hours, which gives 71 intervention periods since the first period is just a fitting period. Importantly, the size of the infusion cannot be greater than 6 Units [11Chase JG, Shaw GM, Lin J, et al. Adaptive bolus-based targeted glucose regulation of hyperglycaemia in critical care Med Eng Phys 2005; 27(1): 1-11., 13Wong XW, Singh-Levett I, Hollingsworth LJ, et al. A novel, model-based insulin and nutrition delivery controller for glycemic regulation in critically ill patients Diabetes Technol Ther 2006; 8(2): 174-90.] for patient safety. To be physiologically realistic it must be also greater than 0. Therefore, the infusion U_{I, target} in Fig. (3) is constrained:

(24)  0≤UI,target≤6U

The results of the algorithm of Fig. (3) for the integral method of Equations (7)-(10) are shown in Fig. (11a), where 7% uniformly distributed noise is placed on the hourly glucose measurements to mimic the sensor error in the glucometer [11Chase JG, Shaw GM, Lin J, et al. Adaptive bolus-based targeted glucose regulation of hyperglycaemia in critical care Med Eng Phys 2005; 27(1): 1-11., 13Wong XW, Singh-Levett I, Hollingsworth LJ, et al. A novel, model-based insulin and nutrition delivery controller for glycemic regulation in critically ill patients Diabetes Technol Ther 2006; 8(2): 174-90., 23Hann CE, Chase JG, Lin J, et al. Integral-based parameter identification for long-term dynamic verification of a glucose-insulin system model Comput Methods Programs Biomed 2005; 77(3): 259-70.]. All measurements in Fig. (11a) lie in the 4-6.1 mmol/L band showing that very tight glucose control is achieved when S_I is constant. Fig. (11b) shows the results of using the derivative method of Equations (11)-(12) in place of the integral method in Step 3 of Fig. (3). Again all measurements lie in the 4-6.1 mmol/L band, showing there is virtually no difference between the methods.

Fig. (11) (a) Algorithm of Fig. (3), with the integral method of Equations (7)-(10)used in Step 3. (b) Algorithm of Fig. (3), with the derivative method of Equations (11)-(12) used in Step 3.

The result of Figs. (11a,b) shows that for Patient 519, the parameter regimes of Equations (18) and (19) that caused instability for the derivative method in Figs. (6) and (7), were not realized. The mean value of Q(t) during this “virtual trial” of patient 519 was 16.5 mU/min. Examining Fig. (6), it can be seen that for these relatively high Q(t) values the derivative and integral methods behave similarly.

The results of Fig. (11) show that with continual insulin infusions over time the derivative and integral methods perform similarly in glucose control with hourly measurements of glucose. Therefore, since the main differences in control have already been investigated in Fig. (10), the comparison of the derivative and integral methods is discontinued in this section.

The insulin sensitivity profile of Patient 519 as fitted in [23Hann CE, Chase JG, Lin J, et al. Integral-based parameter identification for long-term dynamic verification of a glucose-insulin system model Comput Methods Programs Biomed 2005; 77(3): 259-70.] is highly dynamic, and the first three days are shown in Fig. (12). To demonstrate the practical aspects of model-based glucose control, the algorithm of Fig. (3) is applied to the time varying S_I of Fig. (12) using the integral method. Note that in [38Lin J, Lee DS, Chase JG, Hann CE, Lotz T, Wong XW. Stochastic modelling of insulin sensitivity variability in critical care Biomed Signal Process Control 2006; 1: 229-42., 39Lin J, Lee D, Chase JG, et al. Stochastic modelling of insulin sensitivity and adaptive glycemic control for critical care Comput Methods Programs Biomed 2008; 89(2): 141-52.] and Lin 2007 [49Lin J. Robust modelling and control of the glucose-insulin regulatory system for tight glycemic control of critical care patients PhD Thesis, University of Canterbury, Department of Mechanical Engineering. 2007.], the integral method has been well validated and proven for extensive numbers of virtual patients, therefore no further patients other than Patient 519 are tested in this paper.

Fig. (12) Time varying SI profile over first 3 days for patient 519.

The nutritional input P(t) is again held constant with all other parameters the same as given in Equations (5) and (22).

Fig. (13) gives the result for the integral method, which shows that glucose control is significantly worse than Fig. (11). The mean glucose and standard deviation of 5.58 ± 1.03 mmol/L with 67.57% of glucose values lying in the 4.0 to 6.1 mmol/L band. Very similar results are obtained for the derivative method, so these results are not shown.

Fig. (13) Model-based glucose control with the time varying SI of Fig. (12) and a fixed nutritional input given in Equation (22).

The reason for this decrease in performance is explained by the insulin infusion graph of Fig. (14). There are significant periods in Fig. (14) where the insulin has reached the maximum of 6 Units/hour so effectively no added, but necessary, control is being applied in these periods and insulin effect is saturated [2Chase JG, Shaw GM, Wong XW, et al. Model-based glycaemic control in critical care - A review of the state of the possible Biomed Signal Process Control [review and tutorial] 2006; 1(1): 3-21., 36Chase JG, Shaw GM, Lin J, et al. Impact of insulin-stimulated glucose removal saturation on dynamic modelling and control of hyperglycaemia Intell J Intell Syst Technol Applic (IJISTA) 2004; 1(1/2): 79-94.]. The solution to this problem has been to vary the nutrition, as well as the insulin [9Lonergan T, Compte AL, Willacy M, et al. A pilot study of the SPRINT protocol for tight glycemic control in critically Ill patients Diabetes Technol Ther 2006; 8(4): 449-62., 10Lonergan T, LeCompte A, Willacy M, et al. A simple insulin-nutrition protocol for tight glycemic control in critical illness: Development and protocol comparison Diabetes Technol Ther 2006; 8(2): 191-206.]. A fully developed and validated method for modulating both the nutrition and insulin in a model-based glycemic control system is detailed in [13Wong XW, Singh-Levett I, Hollingsworth LJ, et al. A novel, model-based insulin and nutrition delivery controller for glycemic regulation in critically ill patients Diabetes Technol Ther 2006; 8(2): 174-90., 14Chase JG, Shaw GM, Lotz T, et al. Model-based insulin and nutrition administration for tight glycaemic control in critical care Curr Drug Deliv 2007; 4(4): 283-96.].

Fig. (14) Control infusion input uI,target in Step 4 of Fig. (3), for the model-based glucose control of Fig. (13).

To demonstrate the essential concept the nutrition is dropped to 40% of the original value, whenever the insulin hits the upper limit of 6 units. A new insulin infusion is then calculated in Step 4 of Fig. (3) for this reduced nutrition. This simple rule results in a significant improvement in glucose control as shown in Fig. (15). The mean glucose is 5.32 ± 0.67 mmol/L with 76.14% of values lying between 4 and 6.1 mmol/L.

Fig. (15) Model-based glucose control with the time varying SI of Fig. (12) and a simply varying nutritional input.

To demonstrate a new clinical application of the methods presented and to further investigate the comparison of the integral versus derivative approaches, a CGMS sensor is included in the model-based glucose control algorithm. The CGMS sensor measures glucose every 5 minutes with a measurement error that can be approximated by the formula [34Chase JG, Hann CE, Jackson M, et al. Integral-based filtering of continuous glucose sensor measurements for glycaemic control in critical care Comput Methods Programs Biomed 2006; 82(3): 238-47.]:

Equation (25) gives a mean absolute error of 14%, which is typical for CGMS sensors [41Goldberg PA, Siegel MD, Russell RR, et al. Experience with the continuous glucose monitoring system in a medical intensive care unit Diabetes Technol Ther 2004; 6(3): 339-47., 50Guerci B, Floriot M, Bohme P, et al. Clinical performance of CGMS in type 1 diabetic patients treated by continuous subcutaneous insulin infusion using insulin analogs Diabetes Care 2003; 26(3): 582-9.].

Blood glucose is still assumed to be measured hourly with a glucocard and 7% uniformly distributed noise in addition to the CGMS for comparison. To account for the extra noise in the CGMS and to give the greatest chance for an averaging effect on the errors, insulin sensitivity S_I is fitted over the prior 1½ hours rather than 1 hour, as was presented in Fig. (3). The intervention period is also shortened to ½ hour to take advantage of the extra measurements from CGMS. The 1½ hour periods ensure that 2 glucocard measurements will always be available to fit S_I when stepping along each interval of ½ hour.

The same algorithm of Fig. (3) is applied, except the integral and derivative methods are implemented over the longer 1½ hour period and the infusion U_I,target in Step 4 of Fig. (3) is updated every ½ hour. A further change that is made is that 7% low frequency modelling error is added to the glucose measurements, as well as the normally distributed error in Equation (25). The final expression for noise is thus defined:

Equation (26) reflects the fact that a higher resolution in measurements, trades off with both a higher amount of sensor error and importantly, modelling error.

The modelling error is caused by potentially missed dynamics in the glucose-insulin model of Equations (1)-(3). Simple oscillations are used as an initial proof of concept since low frequency oscillations have been often observed in both glucose and S_I [23Hann CE, Chase JG, Lin J, et al. Integral-based parameter identification for long-term dynamic verification of a glucose-insulin system model Comput Methods Programs Biomed 2005; 77(3): 259-70.]. However, further work must be done on real CGMS data to fully characterize the tradeoff’s in the error.

Fig. (16) shows the resulting glucose control for Patient 519 using the same parameters as used for Fig. (15). A significant improvement can be seen with a mean glucose of 5.03 ± 0.42 mmol/L and 98.55% of glucose values lying between 4 and 6.1 mmol/L.

Fig. (16) Model-based glucose control using the integral method in Step 3 of Fig. (3), with the combination of a CGMS sensor and glucocard.

Fig. (17) shows the first 12 hours of data with both the CGMS data and glucocard data plotted against the true glucose. The “true glucose” is denoted by the solid line and includes the modelling error of Equation (26), but not the sensor error, so that δ = 0 in Equation (24). The widely spread points are the simulated CGMS data which are plotted every 5 minutes using the formula in Equation (24) with δ normally distributed as given in Equation (23). The circles are the simulated glucocard hourly “measurements” which put 7% random uniformly distributed noise on the “true glucose”.

Fig. (17) The first 12 hours (720 minutes) of patient 519, with simulated CGMS data shown as points, glucocard “measurements” shown in circles and a solid line denoting the “true glucose” which includes the modelling error of Equation (26), but not the sensor error.

The derivative method is now used in place of the integral method in Step 3 of Fig. (3). The same data is used, but to potentially assist the derivative method, the data is smoothed several times by a 3 point moving average. However, even with smoothing to remove most of the local noise, a significantly worse result is seen in Fig. (18). The mean glucose is 5.5 ± 1.1 mmol/L with only 64.86% of glucose values lying between 4.0 and 6.1 mmol/L. Thus, the derivative method is unable to take advantage of the extra CGMS data, where the integral method gives significantly better outcomes on glucose control despite the larger noise distribution for these sensors. The derivative method clearly performs better when there is very minimal modelling error and the true glucose is close to a straight line between two points, which was the case in Figs. (11,13,15), but does not occur with significant sensor noise and/or modelling error, both of which typically exist.

Fig. (18) Model-based glucose control using the derivative method instead of the integral method in Step 3 of Fig. (3), and with a CGMS sensor and glucocard.