Robust closed-loop control of plasma glycemia: A discrete-delay model approach

The paper investigates the problem of tracking a desired plasma glucose evolution by means of intra-venous insulin administration. A model-based approach is followed, according to a recent model of the glucose/insulin regulatory system which consists of a discrete-delay nonlinear differential equation model. A disturbance is added to the insulin kinetics in order to model uncertainties concerning both the insulin delivery rate and the mechanism actuating the insulin pump. A feedback control law which yields input-to-state stability of the closed loop error system with respect to such a disturbance is provided, which depends on the glucose and insulin measurements at the present and at a delayed time. In silico simulations validate the theoretical results.


1.
Introduction. The design of insulin infusion devices able to control plasma glucose concentration is of great importance when attempting to gain control of decompensated hyperglycemia in selected clinical situations (e.g. perioperative control of glycemia in a decompensated, acutely ill diabetic patient undergoing emergency surgery). From an applicative point of view, different therapeutic schemes can be considered, according to the accuracy of the glucose-insulin model adopted and to the technology available in actuating the designed control law. Glucose control strategies are mainly actuated by subcutaneous or intravenous injections or infusions. Other drug delivery methods are still under investigation, even if recently the Food and Drug Administration (FDA) has approved a device for the delivery of a powder form of insulin by inhalation, as an alternative to subcutaneous injections, [26]. Control of glycemia by means of subcutaneous insulin injections, with the dose adjusted on the basis of capillary plasma glucose concentration measurements, is by far more widespread than control by means of intravenous insulin, since the dose is administered by the patients themselves (see [4] and references therein). However, only open loop or semiclosed loop control strategies can be used, mainly due to the problem of modeling accurately the absorption from the subcutaneous depot in the plasma circulation (see e.g. [28] for a critical review of subcutaneous absorption models, [9] for a model of intra/inter subject variability of the absorption of subcutaneous insulin preparations and the recent paper [22] dealing with ODE models of rapid/long-acting insulin analogues). On the other hand, the use of intravenous insulin administration, delivered by automatic, variable speed pumps under the direct supervision of a physician, provides a wider range of possible strategies and ensures a rapid delivery with negligible delays (see [37] and references therein for a survey of the intravenous route to plasma glucose control).
A closed loop control strategy may be implemented according to a model-less or model-based approach. The former approach does not use a mathematical model of the glucose-insulin system, and provides a control rule for the insulin infusion rate basically according to the experimental data: the earlier empirical algorithms are due to Albisser et al. [2,3]; recent papers on this topic are mainly devoted to the application of PID controllers aiming to mimic the pancreatic glucose response (see e.g. [6,44,16,27]). On the other hand a model-based approach involves a more complete knowledge of the physiology of the system under investigation. The advantages of a model-based approach are evident since, by using a glucose/insulin model, the control problem may be treated mathematically and optimal strategies may be determined. Clearly, the more accurate the model is, the more effective will the control law be. Model-based glucose control has been mainly developed for the Ackerman's linear model [1] (e.g. adaptive control [30], optimal control [46,13], ∞ control [20]); more recently, different approaches have been proposed, based on nonlinear models such as the Minimal Model [5,47], or more exhaustive compartmental models [10,43,18] (e.g. Model Predictive Control [38], nonlinear Model Predictive Control [17], Neural Predictive Control [48], ∞ control [39], non-standard ∞ control [7,40], feedback linearization [32,31]). It has to be stressed that most of these approaches are based on the approximation of the original nonlinear model, provided by linearization, discretization and model reduction (balanced truncation). An excellent review of the available models presently adopted for blood glucose regulation as well as the closed loop control methodologies and technical devices (blood glucose sensors and insulin pumps) may be found in [8] and references therein.
In this paper, a model-based closed loop control scheme is proposed. Differently from the previously mentioned model-based approaches, which use nonlinear Ordinary Differential Equation (ODE) models, the one presented here uses a nonlinear discrete Delay Differential Equation (DDE) model of the glucose/insulin system [33,36]. Since [11], several DDE models have been published, mainly devoted to better represent pancreatic Insulin Delivery Rate (IDR): see e.g. the survey [25] or the very recent [50], where a family of DDE model is used to explain low-frequency glucose/insulin oscillations and their influences in insulin-dependent glucose uptake (last of a chain of papers on such a topic [45,23,24,21,49]). It has to be stressed that when attempting to design a closed loop glucose control, the works published so far have concentrated on Type 1 diabetic patients (who have essentially no endogenous insulin production), avoiding in this way the need to take IDR into account. In the present work we can take into account spontaneous pancreatic IDR, thereby treating healthy, Type 2 diabetic and Type 1 diabetic patients, because the glucose/insulin model we use to represent the natural dynamics of the system [36] does include a spontaneous glucose-dependent insulin delivery rate. This model has been shown to exhibit a number of desirable characteristics: it conforms to established physiological concepts (e.g. pancreatic insulin secretion rate is limited), exhibits satisfactory properties of the solutions [33] (e.g. positivity and boundedness of solutions, local attractivity of a single positive equilibrium), is statistically robust, in that its parameters are statistically identifiable with satisfactory precision by means of standard perturbation experiments, such as the Intra-Venous Glucose Tolerance Test (IVGTT) [36].
The proposed control aims to track a lower desired glucose reference level by means of intravenous insulin infusion, according to a given smooth glucose trajectory. A disturbance affecting the insulin kinetics is also assumed. The followed approach is based on the theory of the input-output feedback linearization with delay cancellation [14,29], and on the theory of the input-to-state stabilization with respect to disturbances adding to the control law [41]. The proposed feedback control law depends on glucose and insulin measurements at the present and at a delayed time, and yields input-to-state stability of the closed loop error system with respect to the disturbance.
The paper is organized as follows. Immediately below basic definitions and notations are reported to clarify the mathematical tools required by the proposed control law. Section 2 deals with the adopted DDE glucose-insulin model, detailing its properties. The control law is designed in Section 3. Simulations are finally reported in Section 4, employing sets of parameters derived from real data. Conclusions follow.
A preliminary version of the paper has been accepted for presentation in [34].
A function : ℝ + → ℝ + is: positive definite if it is continuous, zero at zero and ( ) > 0 for all > 0; of class G if it is continuous, zero at zero, and nondecreasing; of class K if it is of class G and strictly increasing; of class K ∞ if it is of class K and it is unbounded; of class L if it monotonically decreases to zero as its argument tends to +∞.
A function : is of class K for each ≥ 0 and ( , ⋅) is of class L for each ≥ 0. As usual, ISS stands for both input-to-state stable and input-to-state stability.
2. The DDE glucose-insulin model. The glucose-insulin model considered here belongs to the family of DDE models described in [33], which has been proven to provide persistent positive bounded solutions for any admissible initial condition.

Denote ( ), [mM], ( ), [pM]
, plasma glycemia and insulinemia, respectively. The model considered consists of a single discrete-delay differential equation system: where: , is the rate of glucose uptake by tissues (insulin-dependent) per pM of plasma insulin concentration; ℎ , [min −1 (mmol/kgBW)], is the net balance between hepatic glucose output and insulin-independent zero-order glucose tissue uptake (mainly by the brain); -, [L/kgBW], is the apparent distribution volume for glucose; -, [min −1 ], is the apparent first-order disappearance rate constant for insulin; -, [min −1 (pmol/kgBW)], is the maximal rate of second-phase insulin release; -, [L/kgBW], is the apparent distribution volume for insulin; -is the apparent delay with which the pancreas varies secondary insulin release in response to varying plasma glucose concentrations. The nonlinear function (⋅) models the Insulin Delivery Rate as: where: -, [♯], is the progressivity with which the pancreas reacts to circulating glucose concentrations. If were zero, the pancreas would not react to circulating glucose at all; if were 1, the pancreas would respond according to a Michaelis-Menten dynamics, where * is the glucose concentration of halfmaximal insulin secretion; when is greater than 1 (as is usually the case), the pancreas responds according to a sigmoidal function; - * [mM] is the glycemia at which the insulin release is the half of its maximal rate; at a glycemia equal to * corresponds an insulin secretion equal to /2.

Remark 1.
The original model analyzed in [36] considered also a first-order insulinindependent glucose uptake term and a second discrete delay affecting the insulindependent glucose uptake in the glucose kinetics, that is: However, as it is detailed in that paper [36], neither nor appeared to improve the model fit to observations. Remark 2. Note that model (1)-(2) may be thought to belong to the family of DDE models of [50], if the glucose/insulin concentrations are within a reasonable range of values, comprising plasma glycemia far from very low values and plasma insulinemia far from very high values. This way, the first assumption allows to approximate the brain/nerve cells insulin-independent glucose uptake of [50] by a constant function, while the second assumption allows to approximate the hepatic glucose release and the insulin degradation of [50] by a constant function and a linear function of insulin, respectively.
It has to be stressed that the DDE model (1) represents equally well healthy subjects and insulin-resistant or severe diabetic patients, changing the parameter values as appropriate. Moreover, it does belong to the class of "minimal models", in the sense that according to a "minimal" set of independent parameters, it allows to very well resemble the physiology of the glucose/insulin kinetics, and is identifiable from data according to standard perturbation experiments (IVGTT), [36].
Stability properties of (1) have been investigated in [33], and are briefly reported below: -there exists a unique equilibrium point ( , ), corresponding to the basal glycemia and insulinemia; -the equilibrium point is globally absolutely stable (i.e. global asymptotic stability independent of the delay ) provided a sufficient condition is satisfied; however, such a condition appears to be rarely satisfied according to a reasonable physiological setting of the model parameters, supported by plausible statistical results -the equilibrium point is locally asymptotically stable provided the delay is lower than a bound computed according to the other model parameters. Differently from the global absolute stability, such a condition is satisfied according to a very wide range of model parameters (in fact, the whole admissible parameter space).
3. The control law. The proposed control law aims to reduce a high basal plasma glucose concentration to a lower level, according to a reference glucose trajectory, by means of intra-venous insulin administration. Basal glycemia/insulinemia will be referred in the following as and , respectively. To do so, a control input ( ) is added as well as a disturbance ( ) (due to often unavoidable actuator errors) to the insulin kinetics, so that system (1) may be rewritten as: The disturbance ( ) is assumed Lebesgue measurable and locally essentially bounded. Note that by applying the theory of input-output feedback linearization with delay cancellation, with respect to the output ( ) = ( ) and the input ( ), the system (4) has full (equal to 2) relative degree of type III (see e.g. [14], Def.2.4, [29], Th.2). Indeed, by ready computation, it is: with: ) .
Let us consider The aim is that the output ( ) = 1 ( ) tracks a desired smooth reference glycemia, named ref ( ). By setting the error , the following equation holds for the error dynamics: By designing the control law ( ) as: where ( ) is a new input which will be chosen in the following, the error dynamics becomes:˙ Remark 3. From a theoretical point of view, the control law (10) can always be computed because, following the same reasoning of [33], Th.1, the glucose evolution is strictly positive for any admissible initial condition.
Then, set: with the gain ∈ ℝ 1×2 such that the matrix: With the choice of the input ( ) given by (10), (12), the error dynamics becomes: According to the case of a continuous and bounded reference ref ( ), as it is reasonable and desirable, a new disturbance ( ) is defined, with the same characteristics of ( ) (Lebesgue measurable and locally essentially bounded), as follows: Then, the error dynamics is described by the following equation: If the disturbance is zero, then the error decays exponentially to zero. If the disturbance is not zero, system (16) is a bilinear system and, by a well known result in [42], it is integral input-to-state stable with respect to the disturbance ( ), that is, there exist a K ∞ function , a KL function and a K function such that the following inequality holds for all ≥ 0: A further significant improvement can be achieved by the result in [41] on the input-to-state stabilization of nonlinear systems. According to the main theorem in [41], consider for system (4) the following feedback control law: is the previously designed control law (10), (12), and 2 ( ) is given by: ] System (21) is not time-invariant because of the presence of the time function ref ( ). Nevertheless here the results in [41], which concern with time-invariant systems, can be applied as well, since it can be easily proven by using Theorem 4.19, p.176 in [19], dealing with time-varying systems, that the control law (18)- (20) yields input-to-state stability of system (21). The time-invariant function ( ) = can be used as an ISS-Liapunov function for system (21). In this case, the functions 1 , 2 , in the above mentioned Theorem 4.19 are given by: with (⋅) and (⋅) the minimum and maximum eigenvalues of a symmetric, positive-definite matrix. Then, according to the given feedback control law (18)- (20), it results that the following ISS inequality holds for the error ( ): where is a KL function and is a K function. The function is very important since it describes some attenuation or amplification of the disturbance effect. In this case, by the above Theorem 4.19: Note that and depend on the choice of the control designer. With the choice of as in the following section, we have seen that, if = , with positive real, the identity matrix, then the larger is, the more little the term is.
4. Simulation results. Simulations have been performed in order to test in silico the proposed methodology. It has to be stressed that the closed-loop control law has been designed without facing the problem of saturations (no upper/lower bounds are assumed for the exogenous insulin delivery); on the contrary it can happen, from a mathematical point of view, that the control law be not feasible due to an unrealistic insulin delivery rate (e.g. because it is too high, or because it takes negative values). Nevertheless, saturation constraints are in fact considered when attempting to validate the proposed insulin therapy by closing the loop in an in silico simulation. For instance, all simulations performed (including the ones reported in this section) have been carried out by switching off the insulin infusion at any time when the suggested insulin delivery rate is negative (a simple empirical solution which is widely used, especially for model-based control algorithms aiming at satisfying a mathematical criterion, see e.g. [13]). At the same time, however, the control parameters are set in order to reduce such a drawback, and indeed we manage to avoid negative insulin infusions as it appears in the following simulations (e.g. by choosing a desired glucose reference trajectory which decreases from a hyperglycemic level down to a reasonable normoglycemia in a smooth fashion within a reasonable time span of about 3 hours). Parameter values are those of an obese, insulin-resistant subject, identified by means of Generalized Linear Square fitting of an IVGTT perturbation experiment [35]. For system (1), some parameters are directly measured, such as and ; others are known and kept constant, such as and * ; others are estimated, such as , , , , ; others are computed according to the algebraic steady-state conditions, such as and ℎ . Case 1. The following parameters have been used (unit measurements are the ones defined in Section 2, therefore, they are here omitted): In the present simulation case, all parameters used were those actually estimated from the IVGTT test conducted on an obese patient (Body Mass Index ≃ 50) [35]: they show high-normal glycemia ( = 5.611) a substantial degree of insulin resistance ( ≪ 10 −4 ) and a sub-normal Insulin Delivery Rate. This picture (moderate hyperglycemia, obesity, insulin resistance, low IDR) is consistent with the picture of a pre-diabetic patient, whose long-standing obesity has induced such a state of insulin resistance for such a long time that pancreatic glucose toxicity is apparent and insulin delivery (which should be above normal to compensate for increased insulin resistance) is progressively failing. This subject would be expected to develop frank Type 2 Diabetes Mellitus within a relatively short time, unless therapeutic maneuvers (first of all weight loss) are not vigorously employed. As far as the control law, we choose such that the matrix has eigenvalues −0.05, The reference signal is chosen such to obtain the plasma glycemia decreasing exponentially from the value of 5.611 to the new value 4.7. We choose the disturbance as ( ) = 0.1 sin(0.05 ). The pictures in Fig.1 are obtained when the control law is designed as in (18)- (20). With this choice of and , the function in (24) is given by ( ) = 0.2283 √ , thus ensuring a good attenuation of the disturbances effects. Note that real and desired plasma glycemias are almost indistinguishable in the picture.  Case 2. In this case we are considering the same patient as in case 1, but with an acute increase in insulin resistance, such as may be due to infection or surgery: in this case the infectious process and the accompanying hormonal changes (such as increased plasma concentrations of cortisol and catecholamines) contrast the action of insulin at the periphery and the falls. This produces a consequent hyperglycemia, to which the pancreas responds by increasing (to the extent possible) insulin secretion and insulinemia. Therefore, the following parameters are changed, keeping unchanged the others of case 1: = 10 −5 , = 7.856, = 204.11.
We choose such that the matrix has eigenvalues −0.07, −0.08, with: The reference signal is chosen such to obtain the plasma glycemia decreasing exponentially from the value of 7.856. to the new value 4.7. We choose the disturbance as ( ) = 0.1 sin(0.05 ). The plots of Fig.2 are obtained when the control law is designed as in (18)- (20). With this choice of and , the function in (24) is given by ( ) = 0.1351 √ , thus ensuring a good attenuation of the disturbances effects.
Also in this case real and desired plasma glycemias are almost indistinguishable in the picture.  Case 3. In this case we again consider the same patient as in Case 1, but allow for a certain length of time (one or two years, say) to have passed without any effective therapy having been followed. In this case, the natural progression of disease has determined the failure of pancreatic insulin secretion and, in the face of unchanged insulin resistance, a dropping insulin concentration. This in turn determines the emergence of severe hyperglycemia and the establishment of a state of frank Type 2 Diabetes Mellitus. Therefore, the following parameters are changed, keeping unchanged the others of case 1: The reference signal is chosen such to obtain the plasma glycemia decreasing exponentially from the value of 10.37. to the new value 4.7. We choose the disturbance as ( ) = 0.1 sin(0.05 ). The plots of Fig.3 are obtained when the control law is designed as in (18)- (20). With this choice of and , the function in (24)   In order to further validate the proposed control algorithm, a final in silico experiment has been carried out, taking into account glucose and insulin noisy measurements with coefficients of variation equal to 1.5% (glycemia) and 7% (insulinemia), as usually happens in clinical experiments [36]. The virtual patient is the same as for Case 3. It has to be stressed that, if we perturb the measurements with noise, the theory developed in the paper does not apply any more, since it is based on a deterministic model, taking into account actuator disturbances and not sensor disturbances; nevertheless, the proposed simulations confirm the effectiveness of the control scheme in this case also (see Fig.4).

Remark 4.
As a final remark, we would like to stress that it would be very interesting to see what happens when applying the proposed control law to many different models of the glucose-insulin regulatory system. The idea of using one model to simulate and another one to design the control law could in fact be tested in silico. In this case a complete parametric identification step should be conducted on the model used to design the control law [36], according to appropriate data sets generated by the model working as the patient model (e.g. [10], [18], [43]).

5.
Conclusions. The control problem of tracking a desired plasma glucose evolution by means of insulin administration has been investigated. A feedback control  law which yields input-to-state stability of the closed loop error system with respect to a disturbance adding to the control law is provided. The feedback control law depends on the glucose and insulin measurements at the present and at a delayed time. The novelty of the paper relies in the use of a DDE model for the glucoseinsulin system, instead of standard ODE. Moreover, no approximations have been adopted to simplify or reduce the original nonlinear model to design the modelbased control law. Performed simulations validate the theoretical results. Future developments will concern the use of a glucose/insulin model closer to reality: for instance, a more reasonable Hepatic Glucose Production (HGP) should depend on the plasma insulin concentration in a decreasing fashion, especially in frameworks different from the IVGTT [50]. Also the construction of a feedback control law by means of the measurement of plasma glucose only, which is the more common situation, is work in progress. To this end, a nonlinear observer will have to be built for the system (4): some important results have already been achieved in [12,15].