This paper investigates the impact of fast parameter identification methods, which do not require any forward simulations, on modelbased glucose control, using retrospective data in the Christchurch Hospital Intensive Care Unit. The integralbased identification method has been previously clinically validated and extensively applied in a number of biomedical applications; and is a crucial element in the presented modelbased therapeutics approach. Common nonlinear regression and gradient descent approaches are too computationally intense and not suitable for the glucose control applications presented. The main focus in this paper is on better characterizing and understanding the importance of the integral in the formulation and the effect it has on modelbased drug therapy control. As a comparison, a potentially more natural derivative formulation which has the same computation speed advantages is investigated, and is shown to go unstable with respect to modelling error which is always present clinically. The integral method remains robust.
Open Peer Review Details  

Manuscript submitted on 2832008 
Original Manuscript  The Impact of Parameter Identification Methods on Drug Therapy Control in an Intensive Care Unit 
Therapy guidance using physiological models is a growing trend in bioengineering [1Andreassen S, Benn JJ, Hovorka R, Olesen KG, Carson ER. A probabilistic approach to glucose prediction and insulin dose adjustment: description of metabolic model and pilot evaluation study Comput Methods Programs Biomed 1994; 41(34): 15365.8Vogelzang M, Zijlstra F, Nijsten MW. Design and implementation of GRIP: a computerized glucose control system at a surgical intensive care unit BMC Med Inform Decis Mak 2005; 5(38): 10.]. In general, the idea is to use parameter identification to identify patient specific parameters then use these parameters to predict future dynamics, and in particular, individual patient response to therapy. For example, glucose control in the intensive care unit (ICU), has been dramatically improved by using a glucoseinsulin model to optimize insulin doses and changes of nutrition [2Chase JG, Shaw GM, Wong XW, et al. Modelbased glycaemic control in critical care  A review of the state of the possible Biomed Signal Process Control [review and tutorial] 2006; 1(1): 321., 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): 44962.14Chase JG, Shaw GM, Lotz T, et al. Modelbased insulin and nutrition administration for tight glycaemic control in critical care Curr Drug Deliv 2007; 4(4): 28396.]. A glucose control protocol SPRINT (specialized reduced insulin nutrition table) has changed clinical practice in the Christchurch Intensive Care Unit [14Chase JG, Shaw GM, Lotz T, et al. Modelbased insulin and nutrition administration for tight glycaemic control in critical care Curr Drug Deliv 2007; 4(4): 28396.]. The result is tight control of blood glucose with a 32% mortality reduction. Parameter identification is thus an important part of the overall process, as the identified parameters affect the overall therapy prediction. There are many methods for parameter identification, most of which are some variation of the standard nonlinear regression [15Carson ER, Cobelli C. Modelling methodology for physiology and medicine. San Diego: Academic Press 2001.]. These methods include gradient descent [16Ahearn TS, Staff RT, Redpath TW, Semple SI. The use of the LevenbergMarquardt curvefitting algorithm in pharmacokinetic modelling of DCEMRI data Phys Med Biol 2005; 50(9): N8592., 17Scalia GM, Greenberg NL, McCarthy PM, Thomas JD, Vandervoort PM. Noninvasive assessment of the ventricular relaxation time constant in humans by doppler echocardiography Circulation 1997; 95: 1515.], Bayesian with many starting points [18Aittokallio T, Gyllenberg M, Polo O, Virkki A. Parameter estimation of a respiratory control model from noninvasive carbon dioxide measurements during sleep Math Med Biol 2007; 24(2): 22549., 19Cobelli C, Caumo A, Omenetto M. Minimal model SG overestimation and SI underestimation: improved accuracy by a Bayesian twocompartment model Am J Physiol 1999; 277(3 Pt 1): E4818.] and hybrid approaches [20RodriguezFernandez M, Mendes P , Banga JR. A hybrid approach for efficient and robust parameter estimation in biochemical pathways Biosystems 2006; 83(23): 24865., 21Alonso H, Mendonca T, Rocha P. A hybrid method for parameter estimation and its application to biomedical systems Comput Methods Programs Biomed 2008; 89(2): 11222.].
The problem with these standard nonlinear regression approaches is that they all typically require many forward solutions and starting points to ensure robustness. In a modelbased therapeutics approach [2Chase JG, Shaw GM, Wong XW, et al. Modelbased glycaemic control in critical care  A review of the state of the possible Biomed Signal Process Control [review and tutorial] 2006; 1(1): 321., 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): 44962.14Chase JG, Shaw GM, Lotz T, et al. Modelbased insulin and nutrition administration for tight glycaemic control in critical care Curr Drug Deliv 2007; 4(4): 28396.] parameter identification can occur every 12 hours over periods of up to one or several weeks, for many patients. A Monte Carlo method, taking into account sensor error and error in fixed population parameters to optimize therapy selection, also significantly increases the number of parameter identifications required each time. Similar Monte Carlo approaches to optimizing such protocols in a virtual patient simulated trial [10Lonergan T, LeCompte A, Willacy M, et al. A simple insulinnutrition protocol for tight glycemic control in critical illness: Development and protocol comparison Diabetes Technol Ther 2006; 8(2): 191206., 14Chase JG, Shaw GM, Lotz T, et al. Modelbased insulin and nutrition administration for tight glycaemic control in critical care Curr Drug Deliv 2007; 4(4): 28396., 22Lotz T, Chase JG, Shaw GM, et al. Monte Carlo analysis of a new modelbased method for insulin sensitivity testing Comput Methods Programs Biomed 2008; 89(3): 2155.] are also computationally intense for the same reasons.
An integralbased parameter identification method has been developed [23Hann CE, Chase JG, Lin J, et al. Integralbased parameter identification for longterm dynamic verification of a glucoseinsulin system model Comput Methods Programs Biomed 2005; 77(3): 25970.] and extended to other physiological systems [24Hann CE, Chase JG, Shaw GM. Integralbased identification of patient specific parameters for a minimal cardiac model Comput Methods Programs Biomed 2006; 81(2): 18192.27Hann CE, Chase JG, Andreassen S, Smith B, Shaw GM. Diagnosis using a minimal cardiac model including reflex actions Intensive Care Med 2005; 31: S18.], that avoids the need for any forward simulations. It can thus dramatically reduce the computation required. These integral methods are therefore well suited to modelbased control applications requiring realtime parameter identification. For example, agitation sedation control [28Rudge AD, Chase JG, Shaw GM, Lee DS, Hann CE. Parameter identification and sedative sensitivity analysis of an agitationsedation model Comput Methods Programs Biomed 2006; 83(3): 211., 29Rudge AD, Chase JG, Shaw GM, et al. Impact of control on agitationsedation dynamics Control Eng Pract 2005; 13(9): 113949.], fluid therapy and inotropic drug administration for improved cardiac management [25Shaw GM, Chase JG, Starfinger C, et al. Modelling the cardiovascular system Crit Care Resusc 2007; 9(3): 2639.] and control of neuromuscular blockade in general anaesthesia [21Alonso H, Mendonca T, Rocha P. A hybrid method for parameter estimation and its application to biomedical systems Comput Methods Programs Biomed 2008; 89(2): 11222., 30Mendonca T, Lago P. PID control strategies for the automatics control of neuromuscular blockade Control Eng Pract 1998; 6(10): 122531.]. This paper investigates different computationally fast formulations that don’t require forward simulations and in particular examines the impact of these methods on physiological modelbased glucose control.
These issues are illustrated and tested with respect to noise and modelling error using an exemplar glucoseinsulin model that has been extensively validated over many clinical trials [2Chase JG, Shaw GM, Wong XW, et al. Modelbased glycaemic control in critical care  A review of the state of the possible Biomed Signal Process Control [review and tutorial] 2006; 1(1): 321., 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): 44962.14Chase JG, Shaw GM, Lotz T, et al. Modelbased insulin and nutrition administration for tight glycaemic control in critical care Curr Drug Deliv 2007; 4(4): 28396., 22Lotz T, Chase JG, Shaw GM, et al. Monte Carlo analysis of a new modelbased method for insulin sensitivity testing Comput Methods Programs Biomed 2008; 89(3): 2155., 23Hann CE, Chase JG, Lin J, et al. Integralbased parameter identification for longterm dynamic verification of a glucoseinsulin system model Comput Methods Programs Biomed 2005; 77(3): 25970., 31Doran CV, Hudson NH, Moorhead KT, et al. Derivative weighted active insulin control modelling and clinical trials for ICU patients Med Eng Phys 2004; 26(10): 85566.40Lotz TF, Chase JG, McAuley KA, et al. Transient and steadystate euglycemic clamp validation of a model for glycemic control and insulin sensitivity testing Diabetes Technol Ther 2006; 8(3): 33846.]. The glucoseinsulin model and methods are tested using retrospective clinical data. Several practical issues that arise in clinical implementation are addressed, to highlight issues of performance and stability.
Finally, a new modelbased control method for metabolic control is presented, that combines a noninvasive continuous glucose sensor (CGMS) [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): 33947.] with current standard glucometer sensors [42Arkray GlucocardTM Test Strip 2 Data Sheet Japan Arkray Inc 2001.]. This method is shown to provide a potentially significant improvement in glucose control in simulation that warrants further clinical investigation in the future.
The glucoseinsulin model is defined [11Chase JG, Shaw GM, Lin J, et al. Adaptive bolusbased targeted glucose regulation of hyperglycaemia in critical care Med Eng Phys 2005; 27(1): 111.13Wong XW, SinghLevett I, Hollingsworth LJ, et al. A novel, modelbased insulin and nutrition delivery controller for glycemic regulation in critically ill patients Diabetes Technol Ther 2006; 8(2): 17490., 23Hann CE, Chase JG, Lin J, et al. Integralbased parameter identification for longterm dynamic verification of a glucoseinsulin system model Comput Methods Programs Biomed 2005; 77(3): 25970.]:
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 timevarying insulin sensitivity (L/mU min); k the parameter controlling the effective half life of insulin (min^{1}); and α_{G} the MichaelisMenten 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 bolusbased targeted glucose regulation of hyperglycaemia in critical care Med Eng Phys 2005; 27(1): 111.13Wong XW, SinghLevett I, Hollingsworth LJ, et al. A novel, modelbased insulin and nutrition delivery controller for glycemic regulation in critically ill patients Diabetes Technol Ther 2006; 8(2): 17490., 23Hann CE, Chase JG, Lin J, et al. Integralbased parameter identification for longterm dynamic verification of a glucoseinsulin system model Comput Methods Programs Biomed 2005; 77(3): 25970.].
For the glucoseinsulin Equations (1)(3), a similar integralbased parameter identification method to [23Hann CE, Chase JG, Lin J, et al. Integralbased parameter identification for longterm dynamic verification of a glucoseinsulin system model Comput Methods Programs Biomed 2005; 77(3): 25970.] is applied. The parameters α_{G}, k, n and p_{G} in Equations (1)(3) are held constant at the population values based on prior studies and sensitivity analysis [23Hann CE, Chase JG, Lin J, et al. Integralbased parameter identification for longterm dynamic verification of a glucoseinsulin system model Comput Methods Programs Biomed 2005; 77(3): 25970.]:
Similarly, the parameter G_{E} is held at the mean glucose of each patient. The nutritional carbohydrate input appearance rate, P(t) in Equations (1) and (4) is also held constant, but may change with respect to time for different patients. The exogenous insulin u(t) is defined:
where u_{I} (mU/min) is the constant infusion rate over 1 hour and u_{B} (mU/min) is the amount of bolus given over one minute. The parameter S_{I} is insulin sensitivity and is assumed unknown. Integrating Equation (1) from 0 to t yields:
Choosing n values of time,
where
where G is a continuous approximation to the measured glucose [23Hann CE, Chase JG, Lin J, et al. Integralbased parameter identification for longterm dynamic verification of a glucoseinsulin system model Comput Methods Programs Biomed 2005; 77(3): 25970.] and the integrals are evaluated by the trapezium rule. Equation (10) can be solved by linear least squares to determine S_{I} as a constant over any period. Thus, S_{I} may be identified as piecewise constant.
For glucose control in the Intensive Care Unit (ICU), Equation (1) is utilized over periods of 1 hour [11Chase JG, Shaw GM, Lin J, et al. Adaptive bolusbased targeted glucose regulation of hyperglycaemia in critical care Med Eng Phys 2005; 27(1): 111., 13Wong XW, SinghLevett I, Hollingsworth LJ, et al. A novel, modelbased insulin and nutrition delivery controller for glycemic regulation in critically ill patients Diabetes Technol Ther 2006; 8(2): 17490.] and glucose is measured on the hour. For two glucose measurements G_{0} = G(0) and G_{60} = G(60), the function G(t) in Equation (10) can be approximated by a straight line [23Hann CE, Chase JG, Lin J, et al. Integralbased parameter identification for longterm dynamic verification of a glucoseinsulin system model Comput Methods Programs Biomed 2005; 77(3): 25970.]. For a given infusion u_{I} or bolus u_{B} in Equation (6), nutritional input P(t) and glucose measurements G_{0} and G_{60}, the solution to Equation (10) determines the required insulin sensitivity. However, note that a similar approach could be used if glucose is measured more frequently.
A similar, potentially simpler, approach to the parameter identification of Equations (7)(10) is to use the original differential Equations (1)(3), rather than an integral formulation. For a given set of values,
where t_{0}=0. The analogous matrix system to Equation (10) is defined:
where
This method applies gradients which is similar in concept to typical gradient descent methods. The major difference is that no forward simulations are required so like the integral method [23Hann CE, Chase JG, Lin J, et al. Integralbased parameter identification for longterm dynamic verification of a glucoseinsulin system model Comput Methods Programs Biomed 2005; 77(3): 25970.] it is a computationally fast way of identifying large numbers of S_{I} or other timevarying parameters.
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.]. Modelbased control of glucose typically starts by taking two measurements G_{0} and G_{60 }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. Modelbased glycaemic control in critical care  A review of the state of the possible Biomed Signal Process Control [review and tutorial] 2006; 1(1): 321., 23Hann CE, Chase JG, Lin J, et al. Integralbased parameter identification for longterm dynamic verification of a glucoseinsulin system model Comput Methods Programs Biomed 2005; 77(3): 25970.]. 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
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_{0}=8 and G_{60}=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 G_{target}=5 mmol/L. The true S_{I} in the first and second hours are defined as S_{I,1}= 0.0008 and S_{I,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
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
The correct u_{I} is denoted u_{I,target} and is defined:
where
where u_{I} is treated as a variable on the y axis. Fig. (2) shows the points
Fig. (2) Plotting the points

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 nonlinear regression least squares (NRLS) gradient descent algorithm was tested rigorously in [23Hann CE, Chase JG, Lin J, et al. Integralbased parameter identification for longterm dynamic verification of a glucoseinsulin system model Comput Methods Programs Biomed 2005; 77(3): 25970.]. 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. Integralbased parameter identification for longterm dynamic verification of a glucoseinsulin system model Comput Methods Programs Biomed 2005; 77(3): 25970.]. 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 twocompartment model Am J Physiol 1999; 277(3 Pt 1): E4818.]. However, this dramatically increases the number of forward simulations. For example in [23Hann CE, Chase JG, Lin J, et al. Integralbased parameter identification for longterm dynamic verification of a glucoseinsulin system model Comput Methods Programs Biomed 2005; 77(3): 25970.], the integral method was 1000 times fast than the NRLS algorithm which started from one initial guess. If 10100 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 twocompartment model Am J Physiol 1999; 277(3 Pt 1): E4818.]), the integral method would be 10000100000 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. Integralbased identification of patient specific parameters for a minimal cardiac model Comput Methods Programs Biomed 2006; 81(2): 18192., 26Starfinger C, Hann CE, Chase JG, et al. Modelbased cardiac diagnosis of pulmonary embolism Comput Methods Programs Biomed 2007; 87(1): 4660.]).
For a modelbased therapeutics approach [2Chase JG, Shaw GM, Wong XW, et al. Modelbased glycaemic control in critical care  A review of the state of the possible Biomed Signal Process Control [review and tutorial] 2006; 1(1): 321., 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): 44962.14Chase JG, Shaw GM, Lotz T, et al. Modelbased insulin and nutrition administration for tight glycaemic control in critical care Curr Drug Deliv 2007; 4(4): 28396.], 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 modelbased 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 modelbased 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 modelbased control algorithm of Fig. (3).
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. Modelbased cardiac diagnosis of pulmonary embolism Comput Methods Programs Biomed 2007; 87(1): 4660.], 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 glucoseinsulin model; and to better explain the importance of the integral in the formulation. Most importantly, this study is done in the context of modelbased therapeutics and glucose control in the Christchurch Hospital Intensive Care Unit.
This section reviews the implementation of the integral method for long term modelbased glucose control and compares the method with the similar approach that is based on the derivative. It thus contrasts the difference in using integrals and derivatives for this type of bioengineering inspired parameter identification. The robustness of each formulation is investigated with respect to measurement noise and modelling error, intervention period, and number of measurements used.
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): 44962., 10Lonergan T, LeCompte A, Willacy M, et al. A simple insulinnutrition protocol for tight glycemic control in critical illness: Development and protocol comparison Diabetes Technol Ther 2006; 8(2): 191206., 14Chase JG, Shaw GM, Lotz T, et al. Modelbased insulin and nutrition administration for tight glycaemic control in critical care Curr Drug Deliv 2007; 4(4): 28396., 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 modelbased 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 integralbased parameter identification method [23Hann CE, Chase JG, Lin J, et al. Integralbased parameter identification for longterm dynamic verification of a glucoseinsulin system model Comput Methods Programs Biomed 2005; 77(3): 25970.] 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 12 weeks [23Hann CE, Chase JG, Lin J, et al. Integralbased parameter identification for longterm dynamic verification of a glucoseinsulin system model Comput Methods Programs Biomed 2005; 77(3): 25970.], 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): 44962., 10Lonergan T, LeCompte A, Willacy M, et al. A simple insulinnutrition protocol for tight glycemic control in critical illness: Development and protocol comparison Diabetes Technol Ther 2006; 8(2): 191206., 14Chase JG, Shaw GM, Lotz T, et al. Modelbased insulin and nutrition administration for tight glycaemic control in critical care Curr Drug Deliv 2007; 4(4): 28396., 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): 3116.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): 4617.] (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): 135967.] 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 24 potentially contaminated measurements but still provides better control to a 46.1 mmol/L or similar target band than the retrospective data patient who is less acutely ill by APACHE II score.
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_{0} with a fixed initial glucose of G_{0} = 5 mmol/L. Two main parameters sets are considered:
In Equation (18), Q_{0} 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 S_{I} 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 S_{I} 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_{0} < 5 mU/L the derivative method gives an S_{I} value that is significantly different from the true value. In fact it becomes nonphysiological and negative for Q_{0} = 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.
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. Integralbased parameter identification for longterm dynamic verification of a glucoseinsulin system model Comput Methods Programs Biomed 2005; 77(3): 25970.] 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 timevarying insulin sensitivity for this period taken from [23Hann CE, Chase JG, Lin J, et al. Integralbased parameter identification for longterm dynamic verification of a glucoseinsulin system model Comput Methods Programs Biomed 2005; 77(3): 25970.]. 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. Integralbased parameter identification for longterm dynamic verification of a glucoseinsulin system model Comput Methods Programs Biomed 2005; 77(3): 25970.]. 
and all the other parameters are defined in Equation (5).
To begin the modelbased 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_{0} = 1 mU/min, Q_{0} = 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 16 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. 
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 16 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) Modelbased 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.
To further test the methods for a longer period and to demonstrate the practical, clinical issues associated with modelbased control, another patient from the retrospective data of [23Hann CE, Chase JG, Lin J, et al. Integralbased parameter identification for longterm dynamic verification of a glucoseinsulin system model Comput Methods Programs Biomed 2005; 77(3): 25970.] 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:
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 modelbased 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 modelbased 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_{0 }= 11.5 mmol/L, I_{0} = 0 mU/min, Q_{0} = 0 mU/min. The target glucose in Step 4 of Fig. (3) is defined:
where for each consecutive hour the initial conditions G_{0}, I_{0} and Q_{0} are taken as the previously calculated G_{60}, Q_{60} and I_{60}, 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 bolusbased targeted glucose regulation of hyperglycaemia in critical care Med Eng Phys 2005; 27(1): 111., 13Wong XW, SinghLevett I, Hollingsworth LJ, et al. A novel, modelbased insulin and nutrition delivery controller for glycemic regulation in critically ill patients Diabetes Technol Ther 2006; 8(2): 17490.] 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:
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 bolusbased targeted glucose regulation of hyperglycaemia in critical care Med Eng Phys 2005; 27(1): 111., 13Wong XW, SinghLevett I, Hollingsworth LJ, et al. A novel, modelbased insulin and nutrition delivery controller for glycemic regulation in critically ill patients Diabetes Technol Ther 2006; 8(2): 17490., 23Hann CE, Chase JG, Lin J, et al. Integralbased parameter identification for longterm dynamic verification of a glucoseinsulin system model Comput Methods Programs Biomed 2005; 77(3): 25970.]. All measurements in Fig. (11a) lie in the 46.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 46.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. Integralbased parameter identification for longterm dynamic verification of a glucoseinsulin system model Comput Methods Programs Biomed 2005; 77(3): 25970.] is highly dynamic, and the first three days are shown in Fig. (12). To demonstrate the practical aspects of modelbased 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: 22942., 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): 14152.] and Lin 2007 [49Lin J. Robust modelling and control of the glucoseinsulin 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 S_{I} 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) Modelbased glucose control with the time varying S_{I} 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. Modelbased glycaemic control in critical care  A review of the state of the possible Biomed Signal Process Control [review and tutorial] 2006; 1(1): 321., 36Chase JG, Shaw GM, Lin J, et al. Impact of insulinstimulated glucose removal saturation on dynamic modelling and control of hyperglycaemia Intell J Intell Syst Technol Applic (IJISTA) 2004; 1(1/2): 7994.]. 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): 44962., 10Lonergan T, LeCompte A, Willacy M, et al. A simple insulinnutrition protocol for tight glycemic control in critical illness: Development and protocol comparison Diabetes Technol Ther 2006; 8(2): 191206.]. A fully developed and validated method for modulating both the nutrition and insulin in a modelbased glycemic control system is detailed in [13Wong XW, SinghLevett I, Hollingsworth LJ, et al. A novel, modelbased insulin and nutrition delivery controller for glycemic regulation in critically ill patients Diabetes Technol Ther 2006; 8(2): 17490., 14Chase JG, Shaw GM, Lotz T, et al. Modelbased insulin and nutrition administration for tight glycaemic control in critical care Curr Drug Deliv 2007; 4(4): 28396.].
Fig. (14) Control infusion input u_{I,target} in Step 4 of Fig. (3), for the modelbased 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) Modelbased glucose control with the time varying S_{I} 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 modelbased 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. Integralbased filtering of continuous glucose sensor measurements for glycaemic control in critical care Comput Methods Programs Biomed 2006; 82(3): 23847.]:
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): 33947., 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): 5829.].
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 glucoseinsulin 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. Integralbased parameter identification for longterm dynamic verification of a glucoseinsulin system model Comput Methods Programs Biomed 2005; 77(3): 25970.]. 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) Modelbased 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”.
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) Modelbased glucose control using the derivative method instead of the integral method in Step 3 of Fig. (3), and with a CGMS sensor and glucocard. 
This paper has reviewed the modelbased therapeutics approach to glucose control that has been developed and put into regular use in the Christchurch Hospital, New Zealand ICU, and investigated the impact of two different fast parameter identification methods. The key point with these parameter identification methods is that unlike typical nonlinear regression approaches, they do not require any forward numerical solutions to identify model parameters. They are also not starting point dependent, and thus provide a major advantage in the implementation of modelbased “virtual” clinical trials as well as significant realtime capability. The two methods considered were a previously developed integralbased patient specific parameter identification method and a similar approach based on the derivative. At first sight it might be expected that the integral and derivative approaches would perform similarly given they are derived from the same underlying differential equation model. However, even without noise significant differences were observed for certain parameter sets and glucose control protocols, with the integral method remaining significantly more robust.
A number of tests were performed on clinically derived data from a patient in the Christchurch ICU. Very little differences were observed between the modelbased glucose control using the integral method compared to the derivative method for this patient. This result is due to the fact that the insulin levels remained quite constant and high throughout, so that the insulin dynamics between measurements were minimal. The resulting glucose response was thus very close to a straight line with very little modelling error. However, when adding extra measurements from CGMS along with low frequency modelling error, the derivative method performed very poorly, and had worse results than without the CGMS. The integral method on the other hand remained robust and gave a significant improvement in glucose control.
The overall results are summarized as follows.
The integral method is an important research tool in the modelbased therapeutics approach. For example the addition of simulated CGMS shows that a potentially significant clinical gain could be achieved with this continuous sensor. However, further investigation with real CGMS data is required to validate these results. The derivative method, went unstable and failed to realize this possible clinical gain, further emphasizing the importance of integrals in the formulation.