|Home | About | Journals | Submit | Contact Us | Français|
We investigated the applicability of linear quadratic Gaussian (LQG) methodology to the subcutaneous blood glucose regulation problem. We designed an LQG-based feedback control algorithm using linearization of a previously published metabolic model of type 1 diabetes. A key feature of the controller is a Kalman filter used to estimate metabolic states of the patient based on continuous glucose monitoring. Insulin infusion is computed from linear quadratic regulator feedback gains applied to these estimates, generally seeking to minimize squared deviations from a target glucose concentration and basal insulin rate. We evaluated in silico subject-specific LQG control and compared it to preexisting proportional-integral-derivative control.
In health, blood glucose (BG) is tightly controlled by a hormonal network that includes the gut, liver, pancreas, and brain, ensuring stable fasting BG levels (~80–100 mg/dl) and transient postprandial glucose fluctuations. Diabetes is a combination of disorders characterized by absent or impaired insulin action, resulting in hyperglycemia. Intensive insulin and oral medication treatment to maintain nearly normal levels of glycemia markedly reduces chronic complications in both type 1 diabetes mellitus (T1DM1) and type 2 diabetes mellitus,2 but may increase the risk of hypoglycemia or even potentially life-threatening severe hypoglycemia, which could result from imperfect insulin replacement reducing warning symptoms and hormonal defenses.3 Consequently, hypoglycemia has been identified as the primary barrier to optimal diabetes management.4 Thus, the primary purpose of diabetes treatment is optimal control of postprandial hyperglycemia while avoiding hypoglycemia, which naturally formulates an engineering optimization problem.
Glucose control has been studied for more than three decades and widely different solutions have been proposed. The earliest work was based on intravenous (IV) glucose measurements and both positive (glucose) and negative (insulin) control actuations. Studies by Pfeiffer and Clemens created systems such as a glucose-controlled insulin infusion system5 or the more well-known Biostator6 that have been used in hospital settings. Both of these regulators were based on a proportional derivative strategy, where the injected insulin is proportional to the difference between a fixed plasma glucose target and the measured plasma glucose, as well as to the rate of change of plasma glucose. At that time, different types of controllers were also designed based on the prediction of glucose, therefore counteracting the inherent inertia of exogenous insulin compared to the endogenous hormones. The major designs can be found elsewhere.7–11 More work followed, spanning a broader range of control theoretic approaches. All systems were based on IV sensing and IV action, and most of them relied on modeling of human physiology. Techniques such as pole placement,12 adaptive control,13 time domain,14 worst-case frequency domain,15 (H∞) and optimization of linear quadratic (LQ) costs16–20 were adapted to the problem of glucose control. More recently, there is significant interest in applying model predictive techniques to the control of T1DM.21,22 For a review, see Bequette.23
With the advent of new technologies in glucose sensing and insulin infusion, it is now possible to observe and act upon glucose levels using real-time measurements (the sampling period of most continuous glucose monitors (CGM) is 5 minutes or less). Therefore, increasing scientific and industrial effort is focused on the development of regulation systems (e.g., artificial pancreas) to control insulin delivery in people with diabetes. While these new technologies open the way to both open- and closed-loop BG control, they also suffer from certain drawbacks.
Advances in surgically implantable intravenous sensors and intraperitoneal (IP) insulin pumps have triggered great interest in the control community.34–36 However, while it is believed that implantable sensors are closer to intravenous sensing (via blood draws) and are therefore less vulnerable to delays and errors, studies have shown that these sensors suffer from delays nearly equivalent to subcutaneous sensors.37 Implantable pumps are also believed to be more efficient than SC pumps, mimicking the natural route of insulin delivery (peritoneal injections) more closely.37 However, all implantable devices require surgery for insertion and have a limited lifetime.36
In summary, recent efforts in regulating glucose homeostasis have explored three major routes. Results on the IV–SC route have been published by Hovorka et al.38 and El-Khatib et al.,39 both focusing on SC insulin injection but accessing glucose concentration via IV measurements and both utilizing model-predictive control methodologies. Hovorka and co-workers38 focused on strictly negative actuation (insulin only). El-Khatib and colleagues39 developed a double actuation scheme (insulin + glucagon). Renard36 (University of Montpellier, France) is developing a glucose control scheme based on the use of an implanted sensor and pump (IP–IP route). Finally, Steil and co-workers40 are developing a fully SC–SC-based glucose regulator based on proportional-integral-derivative (PID) methodology.
In this article we proposed a methodology for the design of closed-loop feedback controllers based on linear quadratic Gaussian (LQG) control, in which optimal insulin injection are computed based on CGM. The LQG control design methodology, developed in the 1960s, has been used in many application domains.41 In LQG control, an actuation signal is computed to minimize squared-error deviations from a nominal operating point, which in the control of diabetes corresponds to tight glycemic regulation around a reference glucose concentration (e.g., 100 mg/dl). An LQG controller comprises two main components—a state observer and a set of feedback gains—both of which are derived from a linear dynamic model of the system being controlled. In modeling the system, the choice of the state vector is made via two antagonistic criteria: the higher the dimension of the model, the more precisely the model can describe observed dynamics and estimate the modeled quantities; however, a high dimension also renders the estimation procedure sensitive to noise in the observed signal (BG), lowering the precision of the state estimate, which can preclude the use of subject-specific regulation due to unobservable states and nonidentifiable parameters.
The Augmented Meal Model (AMM). The baseline model for glucose–insulin kinetics employed in this article is the oral glucose “meal” model of Dalla Man et al.,42 which by construction represents glucose and insulin fluxes during a meal. The structure of the meal model42 includes a nonlinear gastrointestinal submodel (three states), two glucose compartments, and five insulin-related states and was validated with triple tracer data collected from more than 200 subjects without T1DM. In Dalla Man et al.,43 the model is modified to reflect the lack of pancreatic insulin production in T1DM. Finally, as described previously,43 we augmented the model to take into account the transport of insulin from subcutaneous injection to blood circulation, and further to the interstitium. We refer to the resulting set of differential equations as the augmented meal model and use the AMM for in silico testing of controllers synthesized from the LQG design methodology.
The Reduced Meal Model (RMM). For control design purposes, we used a reduced version of the AMM, in which we captured the effect of oral glucose via a single equation. Specifically, we modeled the glucose rate of appearance Ra (mg/kg/min) as a first-order lag of the meal disturbance D (mg/kg/min):
where τmeal is the time constant associated with glucose absorption from the gut. The glucose rate of appearance term appears in the type 1 differential equation for plasma glucose Gp (mg/kg) as
where Gt (mg/kg) refers to tissue glucose, Uii (mg/kg/min) refers to insulin-independent glucose utilization, and Id (pmol/liter) refers to the delayed insulin signal associated with endogenous glucose production, where we have assumed there is no glucose renal excretion for the target range of Gp under closed-loop control. The resulting set of differential equations has 11 state variables.
Subject-Specific Parameters for AMM and RMM. Because people are widely different from one another, it is important to tailor both AMM and RMM to the physiology of a particular person. Some model parameters are readily available (e.g., body weight), but most are highly model specific and require data collection before the controller can be activated. For a complete person-specific model we need to estimate most of the model parameters, which will allow the estimator not only to be bias free at steady state, but also capture the complete dynamics of the system. In practice, parameter estimation is cumbersome and must be based on the analysis of glucose, insulin, and BG data collected during standard clinical glucose tolerance tests. The simulation results of this article present the ideal case where, for individual subjects, we assumed complete knowledge of all AMM parameters.
The LQG controller developed consists of two main parts: (i) a state observer (Kalman filter) based on the RMM with subject-specific parameters, linearized around the desired BG operating point, and (ii) a set of linear quadratic regulator (LQR) feedback gains computed from the linearized RMM to minimize a least-squares criterion. Both aspects of the controller are outlined.
Physiological State Estimation. The continuous-time dynamics of the linearized RMM can be expressed in succinct form, as follows:
Estimates of the state vector X can be computed dynamically from the observer equation:
where L represents the innovation gain that causes the state estimate to deviate from the open-loop prediction based on the difference between what is actually measured and what was predicted to be measured, . The choice of the innovation gain matrix L is critical to the method. In LQG control, we computed L as the optimal Kalman filter gain, assuming (for control synthesis only) that both w and v are white noise processes with covariance matrices chosen to reflect the magnitude of the disturbances and sensor noise. Computationally, this amounts to the solution of an algebraic Riccati equation.
LQ Regulator. The second component of the LQG controller is the set of LQR feedback gains, which transform our estimate of the physiological state of the system into insulin dosing recommendations, as follows
Note that since u(t) represents the deviation from the reference insulin injection, the actual command to the insulin pump is , where [x]+ refers to the nonnegative part of x [i.e., [x]+ = max(0,x)], accounting for the constraint that negative boluses cannot be injected. The choice of Ka is critical. Following LQG methodology, we computed Ka through the solution of another Riccati equation so that minimizes the deterministic objective function:
where (i) the state is constrained according to the dynamical equations and (ii) Q is a positive semidefinite matrix of weights that penalize state deviations away from the reference operating point of the controller. In general, matrix Q allows the control designer to specify the states that are of interest (glucose vs insulin, plasma concentration vs other compartment) and the aggressiveness of the regulator (how fast it will try to reach equilibrium, possibly undershooting the target value). For the experimental results of this article, we have chosen Q to be a diagonal matrix, with a common weight q applied to all glucose states, and a weight of one applied to all insulin states. The best value of q depends on the reliability of the sensor and on subject characteristics.
Discretization. Our discussion so far has focused on LQG control as a continuous-time process; however, in practice, the LQG controller must be implemented as a sampled data format because CGM measurements are only available at discrete intervals and because existing insulin pumps only allow discrete changes to basal patterns or discrete bolus injections. Envisioning scenarios in which sampling intervals change over time, our approach was to design a continuous time LGQ controller, as outlined above, which we then discretized for final implementation. (We did not synthesize a discrete time LQG controller based on a discrete representation of the RMM.) To discretize the state observer, we zero order held the samples and used a piecewise constant measurement signal as input to the continuous time Kalman filter. Next, in computing the bolus associated with each CGM sample, we used a predictive technique in which the bolus amount is computed as an estimate of the amount of insulin that the LQG controller would inject over the subsequent sampling interval.
We evaluated, in silico, subject-specific LQG control and compared it to a previously reported PID controller.40 The purpose of this study was to evaluate the ideal limits of performance associated with LQG control, assuming complete knowledge of all AMM parameters for an in silico population of type 1 diabetic individuals. The evaluation of LQG and comparison to PID was performed in a simulation environment based on the AMM and is similar to the GIM simulator of Dalla Man and colleagues.43 Our simulation environment included a random CGM sensor noise model developed at the University of Virginia. Our simulation experiments followed the in vivo protocol described by Steil et al.,40 with both CGM sampling and closed-loop control boluses every minute.
Figure 1 illustrates a typical scenario in which both the PID controller and the LQG controller are designed around a target glucose concentration of 100 mg/dl. The in silico subject is described by a set of parameters for the AMM, including a basal plasma glucose mass of 250.1 (mg/kg) and a basal plasma insulin mass of 3.539 (pmol/kg). The insulin clearance of this subject is 1.240 (min-1) and, while not a parameter of the AMM, the insulin sensitivity is approximately 0.0003 (ml/min/μU). Gains of the PID controller, along with computation of the proportional, integral, and derivative signals, are implemented as closely to the methodology of Steil et al.40 as possible using the AMM to assess subject-specific metabolic parameters, as needed. The LQG controller for this subject was synthesized with q = 1 (so that deviations in glucose and insulin states are weighted equally) and r = 1, assuming a process noise covariance of 1 and a sensor noise covariance of 0.001. For both controllers, the subject is assumed to be at target BG, in steady state, at hour zero of the protocol. Then, as shown in the third subplot, the subject experiences a glucose load at breakfast (hour 2, 55 grams carbohydrate), lunch (hour 7, 88 grams), dinner (hour 12, 69 grams), and a snack (hour 16, 45 grams). While the simulated meal disturbance is deterministic, neither the LQG nor the PID controllers have prior knowledge of meals (i.e., no “meal announcement”). More specifically, for the LQG controller, the meal disturbance signal is not an input to the Kalman filter. The first subplot in Figure 1 shows the BG concentration resulting from each controller, with PID shown in blue and LQG shown in red. Both controllers are subject to the same sensor noise signal, shown in the second subplot (Figure 1). Note that whereas both controllers experience similar postprandial excursions, the LQG controller, by virtue of its ability to estimate insulin states, generally manages to avoid hypoglycemic dips after meals. (An exception to this occurs in the recovery period after the snack, where the subject experiences a nadir of 90 mg/dl.) The fourth subplot of Figure 1 shows commanded insulin as a pump rate for both the LQG controller (in red) and the PID controller (in blue). Note that the LQG controller tends to stop insulin injections sooner than the PID controller in the recovery period after each meal. Finally, note that both controllers start injecting insulin in hour 22 of the protocol, when the subject is already below the target glucose concentration of 100 (mg/dl); this is a consequence of degradation of the sensor signal at hour 22.
To understand the relative performance of LQG versus PID, we designed a simulation experiment involving 10 independent trials (as described earlier) for each of 100 simulated subjects. For each subject, a PID controller was designed according to the specifications of Steil et al.40 using the AMM to estimate subject-specific metabolic parameters for determining PID coefficients. After running all 1000 trials of the PID controller, we computed the average closed-loop glucose concentration (across all subjects, for the duration of the protocol). Next, for each of the 100 subjects, again assuming full knowledge of each subject's AMM parameters, we designed an LQG controller. We used a common target glucose concentration for all subjects, identified through an iterative process, so that the LQG approach would achieve the same average glucose concentration achieved by the PID approach. Results from this study are presented for four traditional indices of glucose control:
Results from the 100 subjects (1000 total trials) are shown in Table 1. Note that the subjects spent less time in hyperglycemia on average (p < 0.001) under PID control compared to LQG. However, the LQG controller achieved much lower hypoglycemic excursions, with a smaller average PERCVL (p < 0.001), smaller average LBGI (p < 0.001), and higher Min_BG on average (p < 0.001). Average glycemia for the PID controller was 128.2 (mg/dl), whereas average glycemia for the LQG controller was 128.7, not significantly different (p = 0.65). [Note that while the average Min_BG for the LQG controller was 81 (mg/dl), a small number of in silico subjects under LQG control experienced hypoglycemic excursions below 70 (mg/dl), resulting in an average PERCVL of 0.3%.]
The current management of diabetes typically uses SMBG to adjust the dosing of insulin delivered via injections or an insulin pump. Glucose is measured at infrequent (5/day) and irregular times, and insulin is injected subcutaneously according to these measures and the estimated amount of ingested carbohydrates. Depending on the treatment strategy, insulin is administered either via an insulin pump (as basal rate and premeal boluses) or via injections typically containing fast-acting and long-acting insulin. In both cases, the relationship between the amount of insulin injected and measured plasma glucose is determined by both the health care practitioner and the patient, based on past experience and an initial rule of thumb (e.g., 1800 rule and 450 rule). Insulin boluses are traditionally calculated in two phases. First, the amount of insulin needed by a person to compensate for the carbohydrate content of an incoming meal is computed. This is done by estimating the amount of carbohydrates to be ingested and multiplying by each person's insulin/carbohydrate ratio. Second, based on the difference between measured BG concentration and the target level, the amount of insulin needed to reach the target is computed by multiplying the (BG - target) difference by the individual's insulin correction factor.
With the advent of CGM technology, these well-established empirical calculations are now being enhanced and reformulated for automatic insulin delivery. The most recent control efforts are focusing on the subcutaneous (SC–SC) route, which relies on available technologies. In addition, advances in implantable sensors and insulin pumps have triggered great interest in the glucose control community.34–36 Both subcutaneous and implantable monitoring and insulin delivery devices have the advantage of frequent BG measurements and nearly continuous insulin adjustment. As the timescale of observation changes from three to four BG readings per day to a BG reading every few minutes, new mathematical methods are employed, generally based on a description of system behavior via a set of differential equations. The most straightforward feedback algorithm is PID control, in which SC insulin injections are computed as a proportional, integral, and differential response to the difference between SC glucose measurements and the target glucose concentration.
This article presented a different control approach based on LQG methodology applied to the T1DM meal model42,43,45 of glucose–insulin kinetics. The rationale for this model-based approach is that it may provide a means to avoid hypoglycemic excursions after meals through the estimation of important metabolic states, particularly insulin states. From our simulation experiments, we observed that, at least with complete knowledge of all AMM parameters, LQG control achieves tight glycemic regulation with minimum hypoglycemic events. In this ideal setting, LQG compares favorably to PID: equal average BG comparable maximum, and significantly lower risks for hypoglycemia. Postprandial excursions were comparable in both methods, hinting at limitations of purely reactive control algorithms.
Further study is required to evaluate LQG control in a nonideal case where AMM parameters for individual subjects must be estimated from clinically available data. In ongoing research, we are developing methods for estimating AMM parameters from clinical data using a combination of population average values for some parameters, along with adjusted parameters values that reflect the subject's steady-state glucose and insulin characteristics, insulin clearance, and correction factors.