Glucose control with incomplete information

Glucose control with incomplete information Alessandro Borri Simona Panunzi Andrea De Gaetano CNR-IASI Biomathematics Laboratory Italian National Research Council (CNR) 00168 Rome, Italy Email: alessandro.borri@biomatematica.it simona.panunzi@biomatematica.it andrea.degaetano@biomatematica.it Pasquale Palumbo Costanzo Manes DISIM Department CNR-IASI Italian National Research Council (CNR) University of L’Aquila 00185 Rome, Italy 67100 L’Aquila, Italy Email: pasquale.palumbo@iasi.cnr.it Email: costanzo.manes@univaq.it Abstract—This paper addresses the problem of glucose control in the presence of sampled measurements. This topic is important in the study and development of the Artificial Pancreas (AP), which is a general expression to describe a set of techniques for the control of the glucose behaviour by means of exogenous insulin administration in diabetic individuals. Differently from most of the approaches available in the literature, we not only assume the lack of insulin measurements, but also the availability of the glucose measurements just at sampling times. An observer is designed for the model-based reconstruction of glucose and insulin trajectories from the glucose samples. On top of that, a feedback algorithm (based on the estimated state) is designed to continuously deliver exogenous intra-venous insulin to the patient. Simulations have been performed in-silico on models of virtual patients identified from real data and in the presence of quantization errors. The preliminary results highlight the potential of the proposed approach. I. I NTRODUCTION Artificial Pancreas (AP) is a general expression to describe a set of techniques for the control of the glucose behaviour by means of exogenous insulin administration in diabetic individuals (see e.g. [1], [2]), usually via subcutaneous or intravenous infusions. Diabetes Mellitus (DM) is a widespread disease (affecting millions of people worldwide) characterized by the lack of insulin production or resistance to insulin action. As a consequence, a number of techniques have been made available in the last few years for the automatic control of glucose levels, which need to cope with the several nonidealities of the underlying dynamics, among which we recall the nonlinear and delayed insulin response [3], [4] and the difficulty in obtaining precise and frequent measurements of the current levels of insulin [5]. A major role in this field is given by the model-based techniques, when the regulator is synthesized by explicitly exploiting the structure of the model equations. To this aim, small-scale minimal models are usually preferred since they allow to provide the analytical solution to the control problem under investigation [6]. Such methods require an internal model of the system to reconstruct the behavior of the unknown variables and to simulate (in-silico) the effect of different control strategies. In this framework, our attention mainly focuses on the DDE minimal model introduced in [7], [8] and exploited e.g. in [9], [10], [11], [12] in different clinical settings. The use of DDE models when seeking model-based glucose control laws applicable not only to Type-1 (lacking insulin release), but also to Type-2 diabetic patients (the great majority of people suffering from DM, in which the exogenous insulin administration complements the endogenous insulin production), is due to the need for modeling the pancreatic Insulin Delivery Rate (IDR), which is known to show irregular variations [4]. In this context, the results in [9] allow to track plasma glycemia down to a safe euglycemic level for a Type-2 diabetic patient, while [10] performs a validation of the DDE model closing the observer-based control loop on a population of virtual patients modeled by a different maximal model, accepted by the Food and Drug Administration (FDA) as an alternative to animal trials for the preclinical testing of control strategies in AP. In this paper, we build on top of the literature reviewed above by considering the problem of glucose control in the presence of sampled measurements. Differently from most of the approaches available in the literature, we not only assume the non-availability of the insulin measurements, but also the availability of the glucose measurements just at sampling times (this happens also for the so-called Continuous Glucose Monitoring, CGM [13]). An observer, based on the construction illustrated in [14], [15], is designed for the modelbased reconstruction of glucose and insulin trajectories from the glucose samples. In particular, with respect to the approach followed in [9], we do not assume the availability of glycemia at all times. On the other hand, the construction in [14], [15] cannot be applied to DDE models, so we restrict our attention to the cases in which the apparent delay in the insulin response is negligible (less than 1 minute). As a consequence, we consider a particular case of the glucose-insulin model developed in [7], [8]. On top of that, a feedback algorithm (based on the estimated state) is designed to continuously deliver exogenous intravenous insulin to the patient. We assume that the delivery rate is changed at the sampling times. In the simulation setting, we consider further non-idealities, by introducing quantization (modeling the possible lack of accuracy of the instrument as well as the analog-to-digital and digital-to-analog conversion processes) both in the measurement and in the control phases, where the goal is to track some desired good trajectories of glucose and insulin. The paper is structured as follows: In Section II, we describe the ODE model of the glucose-insulin system. In Section III, we describe an observer reconstructing the continuous glucose-insulin behaviour from sparse glucose measurements. In Section IV, we construct a control law for the exogenous insulin rate input in order to track desired glucose trajectories. Section V shows some simulations, based on data taken from real patients, in an experimental setting. Section VI offers concluding remarks. II. A N ODE MODEL OF THE GLUCOSE - INSULIN SYSTEM WITH SAMPLED MEASUREMENTS The work in [7], [8] introduces a DDE model of the glucoseinsulin system, characterized by a nonlinear delayed secondary insulin release in response to varying plasma glucose concentration, quantified in a constant delay τg > 0. The interested reader is referred to [7] for identification issues and statistical robustness of this model and to [8] for further details on its mathematical properties and qualitative behavior. In general, it is well known that the state space of systems with delays in the state equations has infinite dimension, which leads to difficulties in the system analysis, controller/observer synthesis and physical implementation. Hence, although the aforementioned model is more accurate, it introduces complications which are not worthwhile in those cases where the apparent delay in the pancreatic response is negligible. The behaviour of such individuals can still be well represented by setting the delay τg = 0, which leads to the following ODE model: ( dG(t) dt dI(t) dt • h(·) is a nonlinear map modeling the endogenous pancreatic Insulin Delivery Rate (IDR) as γ h(G) = where γ (dimensionless) is the progressivity with which the pancreas reacts to circulating glucose concentrations and G∗ [mM] is the glycemia at which the insulin release reaches half of its maximal rate; • u(t) is the exogenous intra-venous insulin delivery rate at time t, i.e. the control input [pM/min]. As already said, the previous model is less general than the one with delay (whose qualitative behavior is described in [8]), since fitting the ODE model on the real data for some patients showing a non-negligible τg leads to progressively larger estimation errors for respectively larger values of the delay τg , which is not considered in (1). On the other hand, an ODE model allows to apply the results in [15] and obtain the state reconstruction when the measurements are available just at some instants. So, given the sampling sequence {ti }, we define the piecewise-constant output function as y(t) = G(ti ) t ∈ [ti , ti+1 ), i = 0, 1, ... which can be equivalently restated as a delayed output y(t) = G(t − δ(t)) t≥0 (2) with the time varying delay δ(t) defined within any two consecutive sampling instants as δ(t) = t − ti , t ∈ [ti , ti+1 ), i = 0, 1, ... (3) where we take t0 = 0. The upper bound on the sampling interval is given by ∆ := max (ti+1 − ti ). i III. O BSERVER FOR THE GLUCOSE - INSULIN SYSTEM WITH SAMPLED MEASUREMENTS T = −Kxgi G(t)I(t) + Vgh , G
TiGmax = −Kxi I(t) + VI h G(t) + u(t), (G/G∗ ) γ, 1 + (G/G∗ ) t≥0 (1) with initial conditions G(0) = G0 , I(0) = I0 , where: • G(t) is the Plasma Glycemia at time t [mM ]; • I(t) is the Plasma Insulinemia at time t [pM ]; • Kxgi is the rate of (insulin-dependent) glucose uptake by tissues per unit of plasma insulin concentration [min−1 pM −1 ]; • Tgh is the net balance between hepatic glucose output and insulin-independent zero-order glucose tissue uptake [min−1 (mmol/KgBW )]; • VG is the apparent distribution volume for glucose [L/kgBW ]; • Kxi is the apparent first-order disappearance rate constant for insulin [min−1 ]; • TiGmax is the maximal rate of second-phase insulin release [min−1 (pmol/kgBW )]; • VI is the apparent distribution volume for insulin [L/kgBW ]; By collecting the state variables within the state vector x(t) = [x1 (t), x2 (t)]T = [G(t), I(t)]T , and denoting by ẋ(t) := dx(t) dt the time derivative of x(t), state equations (1)– (3) can be rewritten in the compact form with   ẋ(t) y(t)   δ(t) = f (x(t)) + Bu(t), t≥0 = Cx(t − δ(t)), = t − ti , t ∈ [ti , ti+1 ), i = 0, 1, ... (4) #  " T −Kxgi x1 x2 + Vgh f1 (x) G , f (x) = = f2 (x) −Kxi x2 + TiGmax h(x1 ) VI  δ(t) ∈ [0, ∆], B = [0 1]T , C = [1 0]. Generally speaking, the problem of asymptotic state observation consists in finding a causal system (asymptotic observer), driven by the pair (u(t), y(t)), that produces a vector variable x̂(t) (observed state) asymptotically converging to the state x(t) (i.e., kx(t) − x̂(t)k → 0). Furthermore, an observer is said to be an exponential observer if there exist µ > 0 and α > 0 such that kx(t) − x̂(t)k ≤ µ e−αt kx(0) − x̂(0)k, (5) for any x(0) and x̂(0) in IRn . In order to design such an observer, we preliminarily define the drift-observability map z = φ(x) as       Cx z1 x1 = φ(x) := z= = z2 f1 (x) Cf (x)   x1 (6) = T −Kxgi x1 x2 + Vgh G and its Jacobian  ∂φ(x) 1 = Q(x) := −Kxgi x2 ∂x  0 , −Kxgi x1 (7) which is non-singular for x1 6= 0. The observer in [14], [15] takes the following form: ˙ x̂(t) = f (x̂(t)) + Bu(t) ∆. Otherwise, a chain of sampled observer is required (full details are given in [15]). IV. G LUCOSE CONTROL ALGORITHM To develop the glucose control algorithm, we follow an input-output linearization approach [16]. The dynamics of the observability map in (6) is the following:   ∂φ(x) ż ẋ = Q(x)(f (x) + Bu) ż = 1 = ż2 ∂x    1 0 f1 (x) = (10) −Kxgi x2 −Kxgi x1 f2 (x) + u    z2 1 0 = , h(x1 ) + u −Kxgi x2 −Kxgi x1 −Kxi x2 + TiGmax VI from which one gets  ż = z2 ,       1 T + K x ż2 = −Kxgi x2 − Kxgi x1 x2 + Vgh Kxi x2 xgi 1 G    T  h(x1 ) − Kxgi x1 u. − iGmax VI (11) By imposing ż2 := v, we get the linearizing feedback as  + e−ηδ(t) Q−1 (x̂(t))K y(t) − C x̂(t − δ(t)) , (8)   −(λ1 + λ2 ) where K = assigns the eigenvalues λ1 < 0, λ1 λ2 λ2 < 0 of (A − KC) to obtain convergence to zero of the estimation error x(t) − x̂(t), with   0 1 A= , 0 0 and η > 0 being a design parameter, giving more weight to recent measurements with respect to the older ones. ˆ T , the By recalling x̂(t) = [x̂1 (t), x̂2 (t)]T = [Ĝ(t), I(t)] explicit form of the observer (8), for t ∈ [ti , ti+1 ), for i = 0, 1, ... is the following:  dĜ(t) ˆ + Tgh  = −Kxgi Ĝ(t)I(t)  dt VG    +e−ηδ(t) (λ1 + λ2 )(G(ti ) − Ĝ(ti )),
ˆ dI(t) ˆ + TiGmax h Ĝ(t) + u(t)  = −Kxi I(t)  dt V I   ˆ K (λ +λ )I(t)−λ  1 λ2 +e−ηδ(t) xgi 1 2 (G(t ) − Ĝ(t )) Kxgi Ĝ(t) i i (9) Under hypotheses that are known to hold for the glucoseinsulin system in (1) (see the technical assumptions described in [9] for the more general DDE case), it is possible to prove the following result, establishing the exponential convergence to zero of the observation error. Theorem 1: [15] Consider the system (1)–(3), with δ(t) ∈ [0, ∆]. Then, for any assigned η > 0, there exists K and a ¯ such that, if ∆ < ∆, ¯ then the system in (9) is positive ∆ a global exponential observer for the system in (1)–(3) such that η is the decay rate of the estimation error (i.e., eq. (5) is verified for some µ > 0 and α = η). Note that the previous result enables the use of the observer (9), provided that the sampling interval is upper-bounded by T u = Kxi x2 − v + Kxgi x2 (−Kxgi x1 x2 + Vgh ) TiGmax G h(x1 ) − VI Kxgi x1 (12) which is computable for positive values of the glycemia x1 , consistently with the fact that in those cases the Jacobian matrix Q(x) in (7) is non-singular. We now need to choose v in order to achieve a desirable behavior. Consider a reference signal for the glucose given by yref (t) = Gref (t) = Gd + (Gb − Gd )e−λt , with λ > 0, which aims at driving the basal glycemia Gb of a subject to a desired value Gd . By setting     yref z1,ref zref = := , ẏref z2,ref one easily obtains       ż1,ref (t) z2,ref (t) −λ(Gb − Gd )e−λt żref (t) = = . = ż2,ref (t) ż2,ref (t) λ2 (Gb − Gd )e−λt We define the error e := z − zref , whose dynamics is     ż − ż1,ref z − z2,ref ė = 1 = 2 = Ae + B(v − ż2,ref ). ż2 − ż2,ref v − ż2,ref Since the pair (A, B) is reachable, it suffices to impose v = He + ż2,ref (13) to obtain the convergence to zero T of the linearized error  −λ3 λ4 assigns the eigenvalues dynamics, where H = (λ3 + λ4 ) λ3 < 0, λ4 < 0 of matrix (A + BH). Note that the control law given in Eq. (12)–(13) is a continuous state-feedback control law, depending on the continuous states x1 , x2 , which are not available in our case, since we only know their estimates x̂1 , x̂2 from (9). In the linear case, the separation principle would guarantee the convergence of the output y(t) = G(t) to its reference value yref (t) = Gref (t). Unfortunately, this is not guaranteed in the non-linear case, in general. However, inspired by the procedure and results followed in [9] (which proves local convergence results), we restate the control law in (12)–(13) in terms of a feedback from the reconstructed state, which leads to the following continuous control law  TiGmax u = max 0, Kxi x̂2 − h(x̂1 )+ (14) VI ) T ) H(ẑ − zref ) + ż2,ref + Kxgi x̂2 (−Kxgi x̂1 x̂2 + Vgh G − Kxgi x̂1       x̂1 x̂ Ĝ where ẑ := , and x̂ = 1 = ˆ is the observer f1 (x̂) x̂2 I output in (9). Note that in (14) we explicitly prevent the possibility of a negative exogenous insulin rate. Parameter Gb Ib Kxgi Tgh VG Kxi TiGmax VI γ G∗ Patient 1 8.89 28.29 7.60 · 10−5 0.0025 0.13 0.10 1.45 0.25 2.28 9 Patient 2 8.67 23.70 9.86 · 10−5 0.0026 0.13 0.06 0.74 0.25 2.55 9 Patient 3 8.56 6.90 5.50 · 10−5 0.0026 0.10 0.25 0.89 0.25 1.55 9 TABLE I N UMERICAL VALUES OF THE MODEL PARAMETERS ( IN THE RESPECTIVE UNITS OF MEASUREMENT ) FOR THE 3 PATIENTS CONSIDERED . Fig. 1. Time course of glucose concentration (top panel) and insulin concentration (bottom panel) for the 3 virtual patients considered. The dashdotted lines represent the basal values of glycemia and insulinemia, while the solid lines are obtained by closing the loop with an Artificial Pancreas. V. S IMULATION R ESULTS In this section, the performance of the algorithm developed in the previous section is evaluated in an experimental setting. We consider the data coming from 3 healthy subjects, whose glucose and insulin data are a subset of the data collected in [7]. Such patients underwent an Intra-Venous Glucose Tolerance Test (IVGTT), which is a perturbation experiment consisting in administering intra-venously a glucose bolus after an overnight fasting period and then sampling plasma glucose and serum insulin concentration during the following 3 hours, with varying sampling time. Glucose and insulin measurements from this experiment are used to identify the parameters of the DDE model in [7], [8], by imposing τg = 0 in the identification procedure. In fact, as already mentioned at the beginning of Section II, we focus on the subjects for which the delay in the glucose action on pancreatic IDR is negligible, since the sample-based observer illustrated in Section III follows the approach in [15] and, unfortunately, there are no theoretical results on such a method applied to systems with state delays. After the preliminary identification, since the considered subjects are not diabetic (but some of them show a picture of pre-diabetic condition), we follow the conceptual procedure in [9] and consider a perturbation of the parameters that simulates a natural progression of the disease towards diabetes. In particular, we reduced the insulin resistance (up to values Kxgi < 104 ) as well as the pancreatic glucose sensitivity TiGmax , and then recomputed consistently (via the algebraic steady-state conditions obtained from the model in Eq. (1)) some of the other parameters, in particular the basal values of glycemia Gb and insulinemia Ib , which constitute the equilibria of (1) in absence of exogenous insulin administration (u = 0). The values of the parameters for the three individuals are summarized in Table I. In the spirit of personalized medicine, the parameters of each model are exploited for the design of the artificial pancreas (the observer-based control law) tailored to the particular subject. On top of the assumptions considered so far, we build a more realistic simulation setup by imposing a quantization error both in the sampling and in the actuation, accounting for the processes of analog-to-digital and digital-to-analog conversion in digital devices. We consider a quantization step of 0.1 mM for the glycemia measurements and 20pM/min for the exogenous Insulin Delivery Rate (IDR). As a consequence, the observer-based controller is initialized with quantization errors. Regarding the sampling time of the glucose measurements, we impose a constant sampling time ti+1 − ti = ∆, for all i. So we can write more simply ti = i · ∆, with ∆ = 5 [min]. This choice is consistent with many devices currently on the market [17]. We assume that each control sample is held for the same interval. Simulations have been carried out by designing the AP for each of the three subjects by properly exploiting the correct model parameters. The observer gain K in (8) and the control gain H in (14) are set to obtain the same closed-loop eigenvalues for all patients: λ1 = −0.5, λ2 = −1, λ3 = −0.4, λ4 = −0.8. The parameter η in (9) is set equal to 5, the target glycemia is equal to Gd = 5 mM, with decay rate λ = 1/30. The results are shown in Figures 1–2 in terms of glucoseinsulin behavior and IDR input. It is readily seen that, starting from the individual basal glycemia levels, the target value Gd is approached in the three subjects within the considered time horizon of 3 hours. Fig. 2. Applied control input in terms of Exogenous Insulin Delivery Rate (IDR) for the 3 virtual patients considered. VI. C ONCLUSION In this paper, we addressed the problem of glucose control with incomplete information. Starting from a minimal ODE model, well representing the dynamic behavior of the glucoseinsulin system in the individuals showing a negligible delay in the pancreatic response, we first designed a sample-based observer reconstructing the full glucose-insulin behavior from sparse glucose measurements. Then, we closed the loop by designing a control law (based on the reconstructed state) with the goal of tracking a desired behavior of glycemia. Preliminary simulations have been performed in silico on virtual patients whose parameters have been identified from real data, in a realistic setting with additional non-idealities given by quantization errors. The results show that the approach can constitute a promising tool in the context of the development of the artificial pancreas. Future research effort will be devoted to possible extensions of the observer-based sampled data glucose control to formally take into account the more general cases of state delays and quantized inputs and outputs. R EFERENCES [1] S.A. Weinzimer, G.M. Steil, K.L. Swan, J. Dziura, N. Kurtz, W.V. Tamborlane, Fully automated closed-loop insulin delivery versus semiautomated hybrid control in pediatric patients with type 1 diabetes using an artificial pancreas, Diabetes care, 31(5), 934-939, 2008. [2] A.M. Albisser, B.S. Leibel, T.G. Ewart, Z. Davidovac, C.K. Botz, W. Zingg, H. Schipper, R. Gander, Clinical control of diabetes by the artificial pancreas, Diabetes, 23(5), 397-404, 1974. [3] A. Makroglou, J. Li, and Y. Kuang, Mathematical models and software tools for the glucose-insulin regulatory system and diabetes: An overview, Appl. Numer. Math., 56, 559-573, 2006. [4] P. Palumbo, S. Ditlevsen, A. Bertuzzi, and A. De Gaetano, Mathematical modeling of the glucose-insulin system: A review, Math. Biosci., 244, 69-81, 2013. [5] P. Palumbo, P. Pepe, S. Panunzi, and A. De Gaetano, Recent results on glucose-insulin predictions by means of a state observer for time-delay systems, “Prediction Methods for Blood Glucose Concentration: Design, Use and Evaluation”, H. Kirchsteiger, J. Jørgensen, E. Renard and L. Del Re Editors, Springer, 2015. [6] F. Chee, T. Fernando, Closed loop control of blood glucose, Berlin, Heidelberg, Springer-Verlag, 2007. [7] S. Panunzi, P. Palumbo, and A. De Gaetano, A discrete single delay model for the intra-venous glucose tolerance test, Theor. Biol. Med. Model., 4, 2007. [8] P. Palumbo, S. Panunzi, and A. De Gaetano, Qualitative behavior of a family of delay-differential models of the glucose-insulin system, Discrete Contin. Dyn. Syst. Ser. B., 7(2), 399-424, 2007. [9] P. Palumbo, P. Pepe, S. Panunzi, and A. De Gaetano, Time-Delay ModelBased Control of the Glucose-Insulin System, by Means of a State Observer, Eur. J. Control, 6, 591-606, 2012. [10] P. Palumbo, G. Pizzichelli, S. Panunzi, P. Pepe, and A. De Gaetano, Model-based control of plasma glycemia: Tests on populations of virtual patients, Math. Biosci., 257, 2-10, 2014. [11] P. Palumbo, P. Pepe, S. Panunzi, A. De Gaetano, Closed-loop glucose control: application to the Euglycemic Hyperinsulinemic Clamp, in Proc. 52nd Conf. Decision and Control, Dec. 2013, 4461–4466. [12] P. Palumbo, G. Pizzichelli, P. Pepe, S. Panunzi, A. De Gaetano, Closedloop control scheme for the Euglycemic Hyperinsulinemic Clamp: validation on virtual patients, in Proc. 19th IFAC World Congress of Automatic Control, 2088–2093, Cape Town, South Africa, 2014. [13] W.V. Tamborlane, R.W. Beck, B.W. Bode, B. Buckingham, H.P. Chase, R. Clemons, et al., Continuous glucose monitoring and intensive treatment of type 1 diabetes, The New England journal of medicine, 359(14), 14641476, 2008. [14] F. Cacace, A. Germani, and C. Manes, An observer for a class of nonlinear systems with time varying observation delay, Syst. Control Lett., 59, 305–312, 2010. [15] F. Cacace, A. Germani, and C. Manes, A chain observer for nonlinear systems with multiple time-varying measurement delays, SIAM J. Control and Optim., 52.3, 1862-1885, 2014. [16] A. Isidori, Nonlinear control systems, Springer Science & Business Media, 1995. [17] D.B. Keenan, J.J. Mastrototaro, G. Voskanyan, G.M. Steil, Delays in Minimally Invasive Continuous Glucose Monitoring Devices: A Review of Current Technology, Journal of diabetes science and technology, 3.5, 1207–1214, 2009.