- Email: [email protected]
Diabetic Gastroparesis Modeling and Observer Design
Proceedings, 2nd IFAC Conference on Modelling, Identification and Controlon of Nonlinear Systems Proceedings, 2nd IFAC Conference Modelling, Identification and Controlon of Nonlinear Systems Proceedings, 2nd IFAC Conference Guadalajara, Mexico, June 20-22, 2018 Modelling, Identification and Control of Nonlinear Systems Available online at www.sciencedirect.com Guadalajara, Mexico, June 20-22, 2018 Modelling, Identification and Control of Nonlinear Systems Proceedings, 2nd IFAC Conference on Guadalajara, Mexico, June 20-22, 2018 Guadalajara, Mexico, June 20-22, 2018 Modelling, Identification and Control of Nonlinear Systems Guadalajara, Mexico, June 20-22, 2018
ScienceDirect
IFAC PapersOnLine 51-13 (2018) 97–102 Diabetic Gastroparesis Modeling and Observer Diabetic Gastroparesis Modeling and Observer Diabetic Gastroparesis Modeling and Observer Diabetic Gastroparesis Modeling and Observer Design Design Design Diabetic Gastroparesis Modeling and Observer Design ∗ ∗∗ Claudia Califano ∗ Emeric Scharbarg ∗∗ Design Claudia Claudia Califano Califano∗∗∗∗ Emeric Emeric Scharbarg Scharbarg ∗∗ ∗∗∗
Nicolas Magdelaine Claude H. Moog ∗ and ∗∗ ∗∗∗ Nicolas Magdelaine and H. Claudia Califano∗∗∗ Scharbarg ∗∗∗ Emeric ∗∗∗ Nicolas Magdelaine and Claude Claude H. Moog Moog ∗∗∗ ∗ ∗∗ ∗∗∗ Nicolas Magdelaine and Claude H. Moog Claudia Califano Emeric Scharbarg ∗ Claudia Califano with DIAG, Universit` aa di Roma, via Ariosto 25, 00184 ∗ Claudia Califano is ∗∗∗ and Claude ∗∗∗ is with DIAG, Universit` via Nicolas H. Moog ∗ Claudia Califano is Magdelaine with DIAG, Universit` a di di Roma, Roma, via Ariosto Ariosto 25, 25, 00184 00184 Roma, Italy. [email protected] ∗ Claudia Califano Roma, Italy. [email protected] is with DIAG, Universit` a di Roma, via Ariosto 25, 00184 ∗∗ Emeric Scharbarg Roma,isItaly. [email protected] with Nantes University Hospital/University of Nantes, ∗∗ ∗ Emeric Scharbarg with Nantes University of Roma, Italy. [email protected] ∗∗ Claudia Califano isis with a diHospital/University Roma, via Ariosto 25, 00184 Scharbarg [email protected] withDIAG, NantesUniversit` University Hospital/University of Nantes, Nantes, ∗∗ Emeric Emeric Scharbarg [email protected] with Nantes University Hospital/University of Nantes, Roma, Italy. [email protected] ∗∗∗ [email protected] Nicolas Magdelaine and Claude H. Moog are with LS2N, UMR CNRS ∗∗ ∗∗∗ Nicolas Magdelaine and Claude H. are LS2N, ∗∗∗ Emeric Scharbarg [email protected] with Nantes University Hospital/University of CNRS Nantes, Magdelaine and Claude H. Moog Moog are with [email protected] LS2N, UMR UMR CNRS 6004, BP 92101, 44321 Nantes Cedex 3, France ∗∗∗ Nicolas 6004, BP 92101, 44321 Nantes Cedex 3, France Nicolas Magdelaine and Claude H. Moog are [email protected] LS2N, UMR CNRS [email protected] 6004, BP 92101, 44321 Nantes Cedex 3, France [email protected] ∗∗∗ Nicolas 6004, BP 92101, 44321 Nantes Cedex 3, France Magdelaine and Claude H. Moog are [email protected] LS2N, UMR CNRS Abstract: Type 1 diabetes results from the lack of endogenous production of insulin by the pancreas. 6004, BP 92101, 44321 Nantes Cedex 3, France [email protected] Abstract: 1 results from the lack of endogenous production of by Abstract: Type Type 1 diabetes diabetes results 4% from the lack of endogenous production of insulin insulin by the the pancreas. pancreas. According to various references, to 12% of diabetic patients are affected by gastroparesis which According to various references, 4% to 12% of diabetic patients are affected by gastroparesis which Abstract: Type 1 diabetes results from the lack of endogenous production of insulin by the pancreas. According to various references, 4% to 12% of diabetic patients are affected by gastroparesis which delays the digestion process. Gastroparesis is characterized by a constellation of gastrointestinal delays the digestion process. Gastroparesis is characterized by a constellation of gastrointestinal According to various references, 4% to 12% of diabetic patients are affected by gastroparesis which Abstract: Type 1 diabetes results from the lack of endogenous production of insulin by the pancreas. delays the indigestion process. Gastroparesis is characterized by the a constellation of gastrointestinal symptoms association with delayed gastric emptying (GE). For first time, a mathematical model symptoms in association with delayed gastric emptying (GE). For the first time, a mathematical model delays the digestion process. Gastroparesis is characterized by a constellation of gastrointestinal According to various references, 4% to 12% of diabetic patients are affected by gastroparesis which symptoms in to association withglycemia delayeddynamics gastric emptying (GE). For the first time, a mathematical model is introduced describe the for this significant class of patients. It is shown to yield is introduced to describe the glycemia dynamics for this significant class of patients. It is shown to symptoms in association with delayed gastric emptying (GE). For the first time, a mathematical model delays the digestion process. Gastroparesis is characterized by a constellation of gastrointestinal is introduced totime describe the glycemia dynamics for this significant class of patients. It is shown to yield yield to a nonlinear delay model designed for estimation and control. to aa nonlinear time delay model designed for and control. is describe glycemia dynamics for this significant class of time, patients. It is shown tomodel yield symptoms in to association with delayed gastric emptying (GE). For the first a mathematical to introduced nonlinear time delay the model designed for estimation estimation and control. to a2018, nonlinear delay the model designed for estimation and control. is describe glycemia dynamics for this significant of patients. shownreserved. to yield © introduced IFACtotime (International Federation of Automatic Control) Hostingclass by Elsevier Ltd. ItAllis rights Keywords: Biomedical systems, Systems biology, time-delay systems, observer observer design, glycemia glycemia to a nonlinear time delay model designed for estimation and control. Keywords: Keywords: Biomedical Biomedical systems, systems, Systems Systems biology, biology, time-delay time-delay systems, systems, observer design, design, glycemia dynamics dynamics Keywords: dynamics Biomedical systems, Systems biology, time-delay systems, observer design, glycemia dynamics Keywords: Biomedical systems, Systems biology, time-delay systems, observer design, glycemia 1. INTRODUCTION rather than increasing its response time. Thus, in this paper, aa dynamics 1. rather its response Thus, in 1. INTRODUCTION INTRODUCTION rather than than increasing increasing itsintroduced response time. time. Thus,diabetic in this this paper, paper, a mathematical model is for those patients mathematical model introduced for patients 1.auto-immune INTRODUCTION rather than increasing response time. Thus,diabetic in this paper, mathematical model is isitsand introduced for those those diabetic patientsa Type 1 diabetis is an disease resulting in the lack subject to gastroparesis it is shown that with the introducType 11 diabetis an disease resulting in lack subject to gastroparesis it is that with model isitsand introduced for those diabetic patientsa INTRODUCTION rather than increasing response time. Thus, in the thisintroducpaper, Type diabetisofis is endogenous an1.auto-immune auto-immune disease resulting in the the lack mathematical subject gastroparesis and itGE is shown shown that with the introducof production insulin by the pancreas. Insulin tion of aatotime delay to model and of nonlinearities dependof production of endogenous insulin by the pancreas. Insulin Type 1 diabetis is an auto-immune disease resulting in the lack tion of time delay to model GE and of nonlinearities dependsubject to gastroparesis and it is shown that with the introducmathematical model is introduced for those diabetic patients of production of endogenous insulin by the insulin pancreas. Insulin ing tion on of aglycemia time delay to insulinemia, model GE and of nonlinearities dependdependent diabetis thus requires exogenous injections and leads to aa good candidate dependent diabetis thus requires exogenous insulin injections of production of endogenous insulin by the pancreas. Insulin ing on glycemia and insulinemia, leads to good candidate tion of a time delay to model GE and of nonlinearities dependType 1 diabetis is an auto-immune disease resulting in the lack subject to gastroparesis and it is shown that with the introducdependent diabetis thuseither requires exogenous insulin an injections ing on to glycemia and insulinemia, leads toAs a good candidate which are performed manually, or through insulin model improve the fit to clinical data. aa consequence, which are performed manually, or an insulin dependent diabetis thuseither requires exogenous insulin injections model improve the fit data. ing on and insulinemia, leads toAs a good candidate of production ofcontrol endogenous insulin by pancreas. Insulin tion of to aglycemia time delay to dependent model GE and of nonlinearities dependwhich arethough performed either manually, orthethrough through an insulin model to improve the fit to to clinical clinical data. As a consequence, consequence, pump. A problem arises whose purpose is to dethe design of a delay observer is appropriate for a pump. A control problem arises whose purpose is to dewhich arethough performed either manually, or through an insulin the design of a delay dependent observer is appropriate for model to improve the fit to clinical data. As a consequence, dependent diabetis thus requires exogenous insulin injections ing on glycemia and insulinemia, leads to a good candidate pump. A though control problem arises whose purpose is to dethe design of a delay dependent observer is appropriate for aa crease glycemia from hyperglycemia after some meal as fast as control purpose. This is new in the area of glycemia regulation. crease from hyperglycemia after meal fast as pump. A control problem arises whose purpose to decontrol purpose. This new the area of regulation. the design of a delay isAsappropriate for a which arethough performed either manually, orsome through anas model to improve theis fit to in data. a consequence, crease glycemia glycemia from hyperglycemia after some meal asisinsulin fast as control purpose. This isdependent new inclinical theobserver area of glycemia glycemia regulation. possible, but avoiding future hypoglycemia. The control has to possible, but avoiding future hypoglycemia. The control has to crease glycemia from hyperglycemia after some meal as fast as control purpose. This is new in the area of glycemia regulation. pump. A though control problem arises whose purpose is to deIn Section 2, the model is detailed step by step and an additional the design of a delay dependent observer is appropriate for a possible, but avoiding future hypoglycemia. The control has to In Section 2, the model is detailed step by step and an additional be positive since, once injected, exogenous insulin can no more be positive since, once injected, exogenous insulin can no more possible, but avoiding future hypoglycemia. The control has to In Section 2, the model is detailed step by step and an additional crease glycemia from hyperglycemia after some meal as fast as focus is on the liver modelling. The body continously needs control purpose. This is new in the area of glycemia regulation. since, once injected, exogenous insulin can no more In be positive withdrawn from the organism. Despite decades of research, focus is on the liver modelling. body needs Section the model is detailed The step by stepcontinously and an additional withdrawn from the organism. Despite decades of be since, once injected, exogenous insulin can nohas more focus is To on2,cope the liver modelling. The body continously needs possible, but avoiding future hypoglycemia. Themainstreams control to energy. with sparse meals during the day, the liver be positive withdrawn from the organism. Despite decades of research, research, the control problem is still widely open. Two in energy. To cope with sparse meals during the day, the liver focus is on the liver modelling. The body continously needs In Section 2, the model is detailed step by step and an additional the control problem is still widely open. Two mainstreams in withdrawn from the organism. Despite decades of research, energy. To cope with sparse meals during the day, the liver be positive since, once injected, exogenous insulin can no more acts as a storage device for the post prandial glucose which the control problem isare stillPID’s widely open. Two mainstreams in acts control of glycemia and MPC (Model Predictive as a storage device for the post prandial glucose which energy. To cope with sparse meals during the day, the liver focus is on the liver modelling. The body continously needs control of glycemia are PID’s and MPC (Model Predictive the control problem is still widely open. Two mainstreams in acts as a storage device for the post prandial glucose which be withdrawn from the organism. Despite decades of research, is not absorbed by the tissues and releases glucose, in the control of glycemia are and MPC (Model Predictive Control) and the latter is aaPID’s special case of state feedback which is absorbed the tissues and releases the acts aTostorage device for the postduring prandial glucose which energy. cope by with meals theglucose, day, thein liver Control) and the latter is special case of state feedback which control of glycemia are PID’s and MPC (Model Predictive is not notas absorbed by thesparse tissues and releases glucose, in the the control problem is still widely open. Two mainstreams in form of glycogen, between two meals to avoid hypoglycemia. Control) and the latter is a special case of state feedback which requires an observer for state estimation. Observers are also form of between two meals to hypoglycemia. is notas by the tissues and glucose, in the acts aglycogen, storage device for the post releases prandial glucose which requires an observer for estimation. Observers are also Control) andglycemia the latterare is astate special caseMPC of state feedback form ofabsorbed glycogen, between two meals to avoid avoid hypoglycemia. control of PID’s and (Model Predictive process is modeled by aa nonlinear function depending requires an observer for state estimation. Observers arewhich also This used in Borri et al. (2017), Facchinetti et al. (2010), Sparacino This process is modeled by nonlinear function depending form of glycogen, between two meals to avoid hypoglycemia. is not absorbed by the tissues and releases glucose, in the used in et al. (2017), Facchinetti et (2010), Sparacino requires an for estimation. Observers arewhich also on Thisglycemia. process Low is modeled bywill a nonlinear function depending Control) andobserver the latter is astate special case of al. state feedback glycemia stimulate release of glycogen used in Borri Borri et to al. (2017), Facchinetti et al. (2010), Sparacino et al. (2007) estimate blood plasma glycemia from the on glycemia. Low glycemia will stimulate release of glycogen This process is modeled by a nonlinear function depending form of glycogen, between two meals to avoid hypoglycemia. et al. (2007) to estimate blood plasma glycemia from the used et to al. estimate (2017), Facchinetti et al.glycemia (2010), Sparacino on glycemia. glycemia will stimulate release of glycogen requires anglycemia observer for state estimation. Observers are also rise in Low insulinemia activates storage of glucose by the et al.in Borri (2007) blood plasma from the while intersticial measurement. while rise insulinemia activates storage of glucose by the glycemia. glycemia will stimulate release ofdepending glycogen This process isismodeled by a test nonlinear function intersticial glycemia measurement. et al.in Borri (2007) blood plasma from the on while rise in in 3Low insulinemia activates storage ofobservability glucose by and the used et to al. estimate (2017), Facchinetti et al.glycemia (2010), Sparacino liver. Section devoted to the of (weak) intersticial glycemia measurement. liver. Section 3 is devoted to the test of (weak) observability and while rise in insulinemia activates storage of glucose by the on glycemia. Low glycemia will stimulate release of glycogen intersticial glycemia measurement. liver. Sectiondesign. 3 is devoted to the test of (weak) observability and In aging diabetic patients, the evolution of this chronic disease et al. (2007) to estimate blood plasma glycemia from the the observer Simulations results are provided in Section In diabetic patients, the of chronic disease the observer design. Simulations results are in liver. Section 3insulinemia is devoted toactivates the test of (weak) while rise in storage ofobservability glucose by and the In aging aging diabetic patients, theofevolution evolution of this thisthat chronic disease observer design. Simulations results are provided provided in Section Section goes together with the risk gastroparesis, is, with deintersticial glycemia measurement. 44thebased on some clinical data. goes together with the risk of gastroparesis, that is, with deIn aging diabetic patients, the evolution of this chronic disease based on some clinical data. the observer design. Simulations results are provided in Section liver. Section 3 is devoted to the test of (weak) observability and goes together with the risk of gastroparesis, that is, with de4 based on some clinical data. layed gastric emptying (GE). According to various references, layed gastric emptying various references, goes together with the (GE). risktheofAccording gastroparesis, that is, with de- 4thebased on some clinical data. results are provided in Section In aging diabetic patients, evolution ofto chronic disease observer design. Simulations layed gastric emptying (GE). According tothis various references, this results into the estimate of 12 to 20% of diabetic patients 2. MATHEMATICAL MODELING this results into the estimate of 12 to 20% of diabetic patients layed gastric emptying (GE). According to various references, 2. MATHEMATICAL MODELING goes together with the risk of gastroparesis, that is, with de4 based on some clinical data. this results into the estimate of 12 to 20% ofdigestion diabetic process. patients 2. MATHEMATICAL MODELING affected by gastroparesis which delays the affected by gastroparesis which delays the digestion process. this results into the estimate of 12 to 20% of diabetic patients 2. MATHEMATICAL layed gastric emptying (GE). According to various references, affected by gastroparesis whichbydelays the digestion process. The glycemia/insulinemia dynamicsMODELING Gastroparesis is characterized a constellation of gastroincan be decomposed into Gastroparesis is characterized by aatoconstellation of affected by into gastroparesis which delays theofdigestion process. glycemia/insulinemia dynamics can be decomposed into this results estimate of 12 20% diabetic patients The 2.subsytem MATHEMATICAL MODELING Gastroparesis is the characterized by constellation of gastroingastroinThe glycemia glycemia/insulinemia dynamics can be blood decomposed into testinal symptoms Shin et al. (2013). the which is fed by the plasma intestinal symptoms Shin et al. (2013). Gastroparesis is characterized by a constellation of gastrointhe glycemia subsytem which is fed by the blood plasma inThe glycemia/insulinemia dynamics can be decomposed affected by gastroparesis which delays the digestion process. testinal symptoms Shin et al. (2013). the glycemia subsytem which is the feddigestion by the blood plasmainto insulinemia and the glucose from subsystem. The testinal symptoms Shin et al. (2013). sulinemia and the glucose from the digestion subsystem. The The majority of diabetic patients – not affected by gastroparesis the glycemia subsytem which is fed by the blood plasma inGastroparesis is characterized by a constellation of gastroinThe glycemia/insulinemia dynamics can be decomposed into sulinemia and the glucose from the digestion subsystem. The The majority of diabetic patients – not affected by gastroparesis insulinemia subsystem consists of the blood plasma compartThe majority of diabetic patients –dynamics not affected by gastroparesis insulinemia subsystem consists the plasma compart– characterized by a digestion which is faster than sulinemia and the glucose fromisof the digestion subsystem. The testinal symptoms Shin et al. (2013). the glycemia subsytem which fed byblood the blood plasma ininsulinemia subsystem consists of the blood plasma compart–– are are characterized by digestion which is than The majority of diabetic patients –dynamics not affected by gastroparesis ment and the subcutaneous compartment which is subject to areinsulin characterized by aadynamics. digestion dynamics which is faster faster than ment and the subcutaneous compartment which is subject to the absorbsion This situation is the other way insulinemia subsystem consists of the blood plasma compartsulinemia and the glucose from the digestion subsystem. The ment and the subcutaneous compartment which is subject to the insulin absorbsion dynamics. This situation is the other way – are characterized by a digestion dynamics which is faster than the controlled insulin infusion. Also the digestion process may The majority of diabetic patients – not affected by gastroparesis the insulin absorbsion dynamics. This situation is the The otherresult way insulinemia the controlled insulin infusion. Also the digestion process may around for diabetic patients affected by gastroparesis. ment and the subcutaneous compartment which is subject to subsystem consists of the blood plasma compartthe controlled insulin infusion. Also the digestion process may around for diabetic patients affected by gastroparesis. The result the insulin absorbsion dynamics. This situation is the other way be modeled by two compartments, the stomach and the duode– are characterized by a digestion dynamics which is faster than around for diabetic patients affected by respective gastroparesis. The result ment be by two and the duodeis that the structural properties of the mathematical the controlled infusion. Alsothe thestomach digestion process may and the subcutaneous compartment which is subject to be modeled modeled byinsulin twoetcompartments, compartments, the stomach and the duodeis the structural properties of the respective around for diabetic patients affected bysituation gastroparesis. The num (Magdelaine al. (2015)). For diabetic patients who are the insulin absorbsion dynamics. is mathematical the otherresult way is that that the structural properties ofThis the respective mathematical num (Magdelaine et al. (2015)). For diabetic patients who are models are dramatically different as far as positive invariant sets be modeled by two compartments, the stomach and the duodethe controlled insulin infusion. Also the digestion process may numsubject (Magdelaine et al. (2015)). For diabetic patients who are models are dramatically different far as invariant sets is that the structural properties ofas to gastroparesis, the response time of the digestion around for patients affected by respective gastroparesis. The result models arediabetic dramatically different asthe far as positive positivemathematical invariant sets not not subject to gastroparesis, the response time of the digestion are concerned Farina et al. (2000). num (Magdelaine al. (2015)). For patients who are be modeled by twoetcompartments, thediabetic stomach and the duodenot subject to gastroparesis, the response time of the digestion are concerned Farina et al. (2000). models are dramatically different as far as positive invariant sets subsystem is smaller than the response time of the insulin subis that the structural properties of the respective mathematical are concerned Farina et al. (2000). subsystem is smaller than the response time of the insulin subnot subject to gastroparesis, the response time of the digestion num (Magdelaine et al. (2015)). For diabetic patients who are subsystem is smaller thanbythe response time of athe insulin subare concerned Farina etmeasured al. (2000). The GE delay may be in clinical practice. This moPatients affected gastroparesis have slow digestion models are dramatically different as far as positive invariant sets system. The GE delay may be measured in clinical practice. This mosystem. Patients affected by gastroparesis have a slow subsystem is smaller than the response time of the insulin subnot subject to gastroparesis, the response time of the digestion Theconcerned GEthe delay may be in clinical practice. dynamics This mo- process. system. Patients affected by gastroparesis have a slow digestion tivates inclusion of aaal.time delay in the digestion A standard approach just consists in considering in this are Farina etmeasured (2000). tivates inclusion of in process. A approach just consists in in this The GEthe delay may be in clinical practice. dynamics This mo- subsystem system. affected gastroparesis have slow digestion is smaller thanbythe response time of athe insulin tivates the inclusion of measured a time time delay delay in the the digestion digestion dynamics process.Patients A standard standard approach just consists in considering considering insubthis tivates inclusion of measured a time delay in the digestion A standard approach just consistshave in considering in this The GEthe delay may be in clinical practice. dynamics This mo- process. system. Patients affected by gastroparesis a slow digestion Proceedings, 2nd IFACofConference on in the digestion dynamics 97 process. A standard approach just consists in considering in this tivates the©inclusion a time delay Proceedings, 2nd IFAC onFederation of Automatic Control) 97 Hosting by Elsevier Ltd. All rights reserved. 2405-8963 2018, IFAC Conference (International Modelling, Identification and Controlon of Nonlinear Proceedings, 2nd IFAC Conference 97 Peer reviewIdentification under of International Federation of Automatic Modelling, and Control of Nonlinear Proceedings, 2nd responsibility IFAC Conference 97 Control. Systems Modelling, Identification and Controlon of Nonlinear 10.1016/j.ifacol.2018.07.261 Systems Identification and Control of Nonlinear Modelling, Proceedings, 2nd IFAC Conference on 97 Guadalajara, Mexico, June 20-22, 2018 Systems Guadalajara, Mexico, June 20-22, 2018 Systems
2018 IFAC MICNON 98 Guadalajara, Mexico, June 20-22, 2018
Claudia Califano et al. / IFAC PapersOnLine 51-13 (2018) 97–102
case that the response time of the digestion subsystem is (much) larger than the response time of the insulin subsystem. This feature is acceptable for the model, but was shown to affect dramatically the control law design of insulin infusion. This motivates to distinguish between the two different populations of patients in opposition to the current control literature, where gastroparesis is not a discriminating criterion.
where G∗ is the limit defining hypoglycemia and is taken equal to 70 mg/dl. α1 ≥ 0 characterizes the higher liver endogenous glucose release in case of hypoglycemia. In the simulations dispayed in Figure 2, α1 is taken equal to 0.5 so that f (G) varies from 1 to 0. When insulinemia rises, the endogenous glucose release decreases Ader et al. (1990) and the liver begins to store glucose in the form of glycogen. So that we get G ˙ . G(t) = −θsi (I p )I p (t) + θ4 Xd (t) + θ1 1 + 2α1 exp − ∗ G (6)
In order to have an overall description of the glycemia/insulinemia dynamics, we have to consider the behaviour of several quantities, and in particular of
At this point, considering the input Xd to be zero, one gets the minimal model during a fasting period. When carbohydrates (CHO) are ingested during a meal, then a two compartment digestion subsystem is modeled as follows. The stomach compartment is fed by the CHO input r(t) (that is the meal): 1 X˙s (t) = − Xs (t) + r(t). (7) θ5 The duodenum compartment dynamics is instead described by the differential equation 1 1 (8) X˙d (t) = − Xd (t) + Xs (t − δ ), θ5 θ5 where the gastric emptying time δ characterizing gastroparesis is displayed. For the majority of diabetics, i.e. without gastroparesis, the response time θ5 of the digestion subsystem is smaller than the response time θ3 of the insulin subsystem. Besides the delay δ , the diabetic patients affected by gastroparesis are also characterized by the relationship θ5 > θ3 which will be assumed in the rest of this paper.
G - the blood plasma glycemia, I p - the blood plasma insulinemia, Isc - the subcutaneous insulinemia, Xs - the amount of assimilated carbohydrates that are transferred from the stomach into the duodenum, • Xd - the amount of assimilated carbohydrates that is transferred from the duodenum to the plasma. Consequently the raise of glycemia when a meal is digested can be taken into account by adding a linear term of the form θ4 Xd with θ4 > 0 in the glycemia dynamics. • u - the injected insulin rate • r - the carbohydrates absorbed through the meal
• • • •
The subcutaneous insulinemia dynamics is described by the differential equation 1 θu I˙sc = − Isc + u, (1) θ3 VI θ3 where θ3 > 0 is the response time of the insulin subsystem. The dynamics of the plasma insulinemia Ip instead reads 1 1 I˙p = Isc − I p . θ3 θ3
(2)
Set now G = x1 , Isc = x2 , I p = x3 , Xs = x4 and Xd = x5 . Then, the complete model herein is as follows. θ2 θu x˙1 = −θ2 x3 + 0.4θ2 x3 exp − x3 x θ1VI 1 +θ1 {1 + 2α1 exp − ∗ } + θ4 x5 G 1 θu x˙2 = − x2 + u θ3 VI θ3 1 1 (9) x2 − x3 x˙3 = θ3 θ3 1 x˙4 = − x4 + r(t) θ5 1 1 x˙5 = − x5 + x4 (t − δ ) θ5 θ5 y = x1 .
The glycemia dynamics can be modeled through the differential equation ˙ = −θsi (I p (t))I p (t) + θ4 Xd (t) + θ1 [1 + f (G)], (3) G(t) where the θsi (I p (t)) represents a sensitivity to insulin. When insulinemia is low, for example in case of catheter obstruction, glycemia rises and the body produces ketones to provide energy from fat. But ketosis decreases insulin sensitivity. As a consequence type 1 diabetic patients have to inject large amounts of insulin to recover both from ketosis and hyperglycemia. The θsi function is approximated by Ip θsi (I p ) = θ2 1 − 0.4 exp(− ∗ ) , (4) I
Identification result
θ2 θu /VI > 0 is the practical constant insulin sensitivity factor used in every day medical monitoring of diabetes. The discriminating insulinemia level I ∗ can be taken equal to the so-called basal rate Ib = (θ1VI )/(θ2 θu ), where θ1 > 0 characterizes the glycogenolysis, i.e. the endogenous glucose release by the liver, as detailed below.
Figure 1 displays identification results with and without delay. The model-fit was made from 17:30 to 06:00 using least square error on the output i.e. glycemia and cross-validation was made from 06:00 to 12:30. It illustrates that taking into account the delay, and setting δ = 0, allows to emulate better the behaviour of the glycemia. It becomes obvious in Figure 1 around time 22:30, as the simulated trajectory tracks much better the data.
In (3), the term θ1 (1 + f (G)) represents the net balance between the liver endogenous glucose release and the insulinindependent glucose consumption (e.g. by the brain). In case of hypoglycemia, the liver endogenous glucose release increases. Similarly to Tolic et al. (2000), Sorensen (1978), f (G) is defined as G f (G) = 2α1 exp(− ∗ ). (5) G
The cross-validation simulates the model with the parameters obtained from the fit. Cross-validation shows that the Carboto-Insulin Ratio could be re-estimated for breakfast which is a typical meal with fast carbs (as at 12:30 simulated glycemia drifts from CGM data). 98
2018 IFAC MICNON Guadalajara, Mexico, June 20-22, 2018
Claudia Califano et al. / IFAC PapersOnLine 51-13 (2018) 97–102
99
x1 θ1 1 + 2α1 exp(− ∗ ) G θu u . V θ ψ(y(t), u(t)) = I 3 0 0 0
For alternative models, the reader is referred for instance to Cobelli et al. (2014) in which delays refer to a second phase insulin secretion and to the CGM and insulin pump technologies, to Palumbo et al. (2013) in which the tissue glucose uptake is delayed with respect to the insulin action and the insulin pancreatic secretion is delayed with respect to the glycemia action (whenever the pancreas insulin secretion is non zero), or to Reiterer et al. (2015) where no delay is considered but the essential features of the glycemia-insulinemia dynamics are modeled by linear terms as in (9).
Using the differential representation approach proposed in Califano et al. (2017) to deal with time–delay systems, one has to consider the differential of the dynamics, that is x1 2θ1 α1 exp(− ∗ )dx1 + θ4 dx5 ∗ G G θ2 θu θ 2 θu + −θ2 + 0.4 1 − x3 exp(− x3 ) dx3 θ1VI θ1VI 1 θu d x˙2 = − dx2 + du θ3 VI θ3 1 1 d x˙3 = dx2 − dx3 θ3 θ3 1 d x˙4 = − dx4 + dr θ5 1 1 d x˙5 = − dx5 + δ dx4 θ5 θ5 dy = dx1
3. OBSERVABILITY AND OBSERVER DESIGN
d x˙1 = −
The measured output is the glycemia G(t) in blood plasma. Elementary differentiations of G(t) with respect to time show that the system is weakly observable, in the sense that from the knowledge of the measured output G(t) and of the inputs u(t) and an estimation of the ingested meal r(t), one can compute the four states (G(t), Isc (t), I p (t), Xd (t)) at time t and the value of the fifth delayed state Xs (t − δ ), for almost all values of δ . In the present Section it is shown how to design an ad hoc observer to estimate the states for (9) as described above.
So, outside hypoglycemia (x1 >> 70 mg/dl and f (G) 0), and assuming that there is no catheter obstruction (x2 >> 0 and x3 >> 0), and during fasting (r = 0, x4 = 0 and x5 = 0), the set of equilibria in the plane (x1 , x3 ) is {(G, Ib )} where Ib (θ1VI )/(θ2 θu ) denotes the basal rate. The blood glycemia may have have any value G which is stabilized by the basal insulin infusion rate Ib as expected from medical practice.
and in compact form d x(t) ˙ = (A0 + A1 δ )dx(t) + Pdr +
It is easily seen that the system is in the form
∂ ψ(y(t), u(t)) dy ∂y
∂ ψ(y(t), u(t)) du ∂u ∂ ϕ(x3 ) + dx3 ∂ x3 dy(t) = Cdx(t) +
x(t) ˙ = A0 x(t) + A1 x(t − δ ) + Pr(t) + ψ(y(t), u(t)) + ϕ(x3 (t)) y(t) = Cx(t)
with
0
0 0 A0 = 0 0 0 0 0 A1 = 0 0
0 −θ2 1 − 0 θ3 1 1 − θ3 θ3
0
θ4
3.1 Observability
0 0 0 , 1 0 0 0 − θ5 1 0 0 0 − θ5 00 0 0 0 0 0 0 0 0 0 0 0 0 , P = 0 , 0 0 0 0 1 1 0 0 00 θ5 0
Setting A(δ ) = A0 + A1 δ , in this section we will consider the observability problem assuming that there is no catheter obstruction, thus neglecting ϕ(x3 ). As it will be clear in the observer design later on, such an approximation does not influence the observer itself since such a nonlinearity affects only the estimation of the glycemia which is the measured variable. The observability matrix is then obtained by considering dy = Cdx(t) ∂ψ d y˙ = CA(δ )dx +C dy ∂y ∂ ψ˙ ∂ψ +C dy d y¨ = CA2 (δ )dx + CA(δ ) ∂y ∂y ∂ψ d y˙ +C ∂y
C = (1 0 0 0 0) .
dy(3) = CA3 (δ )dx ∂ ψ˙ ∂ ψ¨ ∂ψ + CA2 (δ ) +CA(δ ) +C dy ∂y ∂y ∂y ∂ ψ (2) ∂ψ ∂ ψ˙ d y˙ +C + CA(δ ) + 2C dy ∂y ∂y ∂y
and the nonlinear terms ϕ(x3 ) and Ψ(y, u) given by θ2 θu 0.4θ2 x3 exp(− x3 ) θ1VI 0 , ϕ(x3 ) = . .. 0 99
2018 IFAC MICNON 100 Guadalajara, Mexico, June 20-22, 2018
Claudia Califano et al. / IFAC PapersOnLine 51-13 (2018) 97–102
400 CGM (mg/dl) - Data CHO - Data Model without delay (mg/dl) Model with delay (mg/dl)
350 300 50 g
250
35 g
200 150
10 g
100
7g
50 0
10
18:00
21:00
00:00
03:00
06:00
8U
2
09:00
12:00
09:00
12:00
8.5 U 2U
101
10
1U
0.90 U/h
0
0.70 U/h Insulin injection (U/h) - Data
10
-1
18:00
21:00
00:00
03:00
06:00
Time
Fig. 1. Model behaviour with and without delay (from 6 pm to 12 noon next day 3.2 Observer design
dy(4) = CA4 (δ )dx ... ∂ψ ∂ ψ˙ ∂ ψ¨ ∂ψ 3 2 dy + CA (δ ) +CA (δ ) +CA(δ ) +C(δ ) ∂y ∂y ∂y ∂y ∂ψ ∂ ψ˙ ∂ ψ¨ 2 d y˙ + CA (δ ) + 2CA(δ ) + 3C ∂y ∂y ∂y ∂ψ ∂ ψ˙ ∂ ψ (3) + CA(δ ) + 3C dy(2) +C dy ∂y ∂y ∂y
Since the model is linear up to a nonlinear input u(t) and output y(t) injection, a standard observer is designed by considering a copy of the model, where the unknown input r(t) is substituted by an estimation rˆ(t) as done by the patient at every meal. That is
Due to the structure of the output and its derivatives, observability can be checked by considering the linear matrix 1 0 0 0 0 0 0 0 θ4 −θ2 C θ4 θ2 θ2 θ4 CA(δ ) 0 − δ − 2 θ θ θ θ 3 3 5 5 = CA (δ ) O = θ θ θ θ 3 2 2 4 4 CA (δ ) 0 2 2 − 2 −2 2 δ 2 θ θ θ θ 3 3 5 5 CA4 (δ ) θ2 θ2 θ4 θ4 0 −3 3 3 3δ − 3 θ3 θ33 θ5 θ5 The determinant is θ 2θ 2 1 4 6 4 1 − + − + det(O) = 2 4 δ θ3 θ5 θ54 θ3 θ53 θ32 θ52 θ33 θ5 θ34
ξ˙ (t) = A0 ξ (t) + A1 ξ (t − δ ) + Pˆr(t) + ψ(y(t), u(t))
(10)
˜ −K(Cξ (t) − y(t)) with the gain K˜ appropriately chosen. Accordingly, denoting by re (t) = rˆ(t) − r(t), the error on the estimation of the ingested carbohydrates, the dynamics of the error e(t) = ξ (t) − x(t) is linear and given by ˜ e(t) ˙ = A0 e(t) + A1 e(t − δ ) − KCe(t) + Pre (t)
(11)
with re (t) acting as an imput on the linear system. Accordingly
and is zero for θ3 = θ5 . It follows that for θ5 = θ3 the system is weakly observable. Recall from Section 2 that we assume θ5 > θ3 due to gastroparesis, so there is no singularity for weak observability. 100
˜ ˜ de + Pdre d e˙ = (A0 + A1 δ − KC)de + Pdre = A(δ ) − KC
For K˜ = 0, the eigenvalues of the matrix A(δ ) = A0 + A1 δ , are
2018 IFAC MICNON Guadalajara, Mexico, June 20-22, 2018
Claudia Califano et al. / IFAC PapersOnLine 51-13 (2018) 97–102
det(A(δ ) − λ I) = −θ2 0
−λ 0 θ4 0 − 1 −λ 0 0 0 θ3 1 1 0 − −λ 0 0 = det θ3 θ3 1 0 0 0 0 − −λ θ5 1 1 δ − −λ 0 0 0 θ5 θ5 2 2 1 1 −λ +λ +λ θ3 θ5 which shows that the system is characterized by four negative eigenvalues and one eigenvalue in zero, so that the system is stable but not asymptotically. It is then necessary to use the gain K˜ to asymptotically stabilize the error dyanamics. One thus has ˜ − λ I) = det(A(δ ) − KC ˜ −λ − k1 0 −θ2 0 θ4 1 0 0 0 −k˜ 2 − − λ θ3 1 1 −k˜ 3 − −λ 0 0 = det θ3 θ3 1 −k˜ 4 0 0 − − λ 0 θ5 1 1 0 0 δ − −λ −k˜ 5 θ5 θ5 θ 1 1 1 2 −(λ + k˜ 1 )( + λ )2 ( + λ )2 + k˜ 2 ( + λ )2 + θ3 θ5 θ3 θ5 θ4 1 1 1 −θ2 k˜ 3 ( + λ )( + λ )2 − k˜ 4 δ ( + λ )2 + θ3 θ5 θ5 θ3 1 1 +k˜ 5 θ4 ( + λ )2 ( + λ ), θ3 θ5
101
account that the measures of the output are not continuous but discrete. The observer’s inputs are: • the measurements of the output y(t) from a monitoring device (CGM) which evaluates the glycemia every 5 minutes; • the open-loop insulin infusion u(t) delivered by the pump; • the amount of carbohydrates (CHO) in the meal rˆ(t) estimated by the patient. The effective amount of CHO is often different from the one estimated by the patient. In Figure 2 the amount of CHO estimated by the patient is assumed to be underestimated at 80 % of the real amount. A meal is taken around 9pm and around 8am (see Figure 1). Though the observer is designed by assuming ϕ(x3 ) = 0, the simulations are carried out the full nonlinear model. Nonetheless, the insulin pump works correctly and the basal rate is never set to zero (c.f. injections on Figure 1), ϕ(x3 ) 0. As a matter of fact such an assumption does not affect the observer behavior, since it fortunately affects only the estimation of the glycemia which is not needed at this point. Simulations were carried out with the following parameters which were identified from clinical data using a standard least squares method. parameters θ1 θ2 θ3 (min) θ4 θ5 (min) δ (min) θu /VI
where the gain matrix coefficients, k˜ 1 , k˜ 2 , k˜ 3 , k˜ 4 and k˜ 5 have to be chosen in order to ensure that the eigenvalues are all in the half left plane. A simple choice consists in setting K˜ = [0.1 0 0 0 0]T . Remark. It should be noted that while the estimation of the system state is characterized by a linear error dynamics described by equation (11), the observer, described by equation (10) is nonlinear due to the presence of the nonlinear term Ψ(y, u). Remark. Weak observability implies that one is able to estimate the state of the given system with some delay, for almost all values of the delays. In the present case, a stronger properties holds true which is the possibility of estimating the state with some delay for all values of the delay.
4. IN SILICO TEST In the present section, simulations are carried out on a diabetic patient’s data who is affected by gastroparesis. The data were obtained from Nantes University Hospital. α1 is set to 0.5, in order to take into account hypoglycemia. In Figure 2, the estimation error on x2 = Isc , x3 = I p , x4 = Xs and x5 = Xd is reported. The observer is designed taking into 101
value 0.32 0.05 70 2.7 22 95 500
The observer was initialized with initial condition z0 = (200, 12, 20, 5, 8)T . The glycemia/insulinemia model was initialized with initial condition x0 = (214, 8.75, 8.75, 0, 0)T .
5. CONCLUSION Despite very recent realizations or announcements, the design of an artificial pancreas is still an open problem. Such a reliable design will go through either a model free control or a model based control for the regulation of the blood plasma glycemia. In the latter case the design of an efficient observer is mandatory to estimate possibly all the state variables of the system. In any case, since the measured output is the intersticial glycemia and the output to be controlled is the blood glycemia, an observer will be the tool to assess the performance of the closed loop. The latter problem was left for future research in this paper as the insterticial glycemia was just assumed to be equal to the blood plasma glycemia. The model was reviewed from scratch and may be developped further. The focus was made on a significant subpopulation of diabetics as they display a very specific behavior. The contributions is this research announcement are as follows • A new time delay model is introduced for the subpopulation of diabetes affected by gastroparesis. From the current literature Borri et al. (2017), time-delay models have shown to be suitable to take into account external delays due to measurement devices or information processing.
2018 IFAC MICNON Guadalajara, Mexico, June 20-22, 2018 102
Claudia Califano et al. / IFAC PapersOnLine 51-13 (2018) 97–102
0
x2 - ξ 2
-2 -4
18:00
21:00
00:00
03:00
06:00
09:00
12:00
03:00
06:00
09:00
12:00
03:00
06:00
09:00
12:00
06:00
09:00
12:00
0 -5
x3 - ξ 3
-10 18:00
21:00
00:00
0 -2 -4
x4 - ξ 4
-6 -8 -10
18:00
21:00
00:00
0 -2
x5 - ξ 5
-4 -6 -8 -10
18:00
21:00
00:00
03:00
Time
Fig. 2. Estimation Error in case of under-estimated amount of carbohydrates In the present paper, the internal delay is inherent to the organism and may be measured through clinical tests or estimated by means of identification of a glycemic holter. • The system is shown to be weakly observable. • An observer is derived and tested in silico.
Future work concerns the use of the information obtained by the observer to automatically compute the necessary insuline quantity needed by the patient. Further perspectives include a model from the blood plasma glycemia to the intersticial glycemia as it is the latter which is directly measured by the sensor. The observer is then extended to estimate the blood plasma glycemia. Though the GE delay is measured getting through some clinical test, it is worth to estimate it in real time as done in Zheng et al. (2011). REFERENCES M. Ader and R.N. Bergman, Peripheral effects of insulin dominate suppression of fasting hepatic glucose production, Am J Physiol Endocrinol Metab 258: E1020-E1032, 1990. A. Borri, F. Cacace A. De Gaetano, A. Germani, C. Manes, P. Palumbo, S. Panunzi and P. Pepe, Luenberger-Like Observers for Nonlinear Time-Delay Systems with Application to the Artificial Pancreas: The Attainment of Good Performance, IEEE Control Systems, 37, pp. 33–49, 2017. C. Califano, L.A. Marquez-Martinez and C.H. Moog, On the observer canonical form for time–delay systems, IFAC World Congress 2011, Milan, Italy, pp. 3855–3860. C. Califano and C.H. Moog, Accessibility on nonlinear timedelay systems, IEEE Trans. Aut. Contr., 62, 3, pp. 12541268, 2017. C. Cobelli, C. Dalla Man, M.G. Pedersen, A. Bertoldo, and G. Toffolo. Advancing Our Understanding of the Glucose System via Modeling : A Perspective, IEEE Trans. Bio. Med. Eng., 61(5), pp. 1577-1592, 2014. 102
A. Facchinetti, G. Sparacino and C. Cobelli, An Online SelfTunable Method to Denoise CGM Sensor Data, IEEE Trans. Bio. Med. Eng., 57, pp. 634–641, 2010. L. Farina and S. Rinaldi, Positive Linear Systems, John Wiley, New York, 2000. M. Krstic, Delay Compensation for Nonlinear, Adaptive, and PDE Systems, Birkh¨auser, Boston, 2009. N. Magdelaine, L. Chaillous, I. Guilhem, J.Y. Poirier, M. Krempf, C.H. Moog and E. Le Carpentier, A Long-term Model of the Glucose-Insulin Dynamics of Type I Diabetes, IEEE Trans. Bio. Med. Eng., 62, pp. 1546–1552, 2015. P. Palumbo, S. Ditlevsen, A. Bertuzzi, and A. De Gaetano, Mathematical modeling of the glucoseinsulin system : A review, Mathematical Biosciences, 244, pp. 69–81, 2013. F. Reiterer, H. Kirchsteiger, A. Assalone, G. Freckmann, and L. del Re. Performance assessment of estimation methods for CIR/ISF in bolus calculators, IFAC Papers Online, 48(20), pp. 231-236, August 2015. A.S. Shin, M. Camilleri. Diagnostic assessment of diabetic gastroparesis, Diabetes, vol. 62, pp. 2667–2673, 2013. G. Sparacino, F. Zanderigo, S. Corazza, A. Maran, A. Facchinetti, and C. Cobelli. Glucose Concentration can be Predicted Ahead in Time From Continuous Glucose Monitoring Sensor Time-Series, IEEE Trans. Bio. Med. Eng., 54(5), pp. 931-629, May 2007. J.T. Sorensen, A physiologic model of glucose metabolism in man and its use to design and assess improved insulin therapies for diabetes, PhD, Univ. of Calif., Berkeley, 1978. I.M. Tolic, E. Mosekilde, and J. Sturis. Modeling the insulinglucose feedback system: the significance of pulsatile insulin secretion, J. Theor. Biol., 207, pp. 361-375, 2000. G. Zheng, J.P. Barbot, D. Boutat, T. Floquet and J.P. Richard, On observation of time-delay systems with unknown inputs, IEEE Trans. Aut. Contr., 56, pp. 1973–1978, 2011.
ScienceDirect
IFAC PapersOnLine 51-13 (2018) 97–102 Diabetic Gastroparesis Modeling and Observer Diabetic Gastroparesis Modeling and Observer Diabetic Gastroparesis Modeling and Observer Diabetic Gastroparesis Modeling and Observer Design Design Design Diabetic Gastroparesis Modeling and Observer Design ∗ ∗∗ Claudia Califano ∗ Emeric Scharbarg ∗∗ Design Claudia Claudia Califano Califano∗∗∗∗ Emeric Emeric Scharbarg Scharbarg ∗∗ ∗∗∗
Nicolas Magdelaine Claude H. Moog ∗ and ∗∗ ∗∗∗ Nicolas Magdelaine and H. Claudia Califano∗∗∗ Scharbarg ∗∗∗ Emeric ∗∗∗ Nicolas Magdelaine and Claude Claude H. Moog Moog ∗∗∗ ∗ ∗∗ ∗∗∗ Nicolas Magdelaine and Claude H. Moog Claudia Califano Emeric Scharbarg ∗ Claudia Califano with DIAG, Universit` aa di Roma, via Ariosto 25, 00184 ∗ Claudia Califano is ∗∗∗ and Claude ∗∗∗ is with DIAG, Universit` via Nicolas H. Moog ∗ Claudia Califano is Magdelaine with DIAG, Universit` a di di Roma, Roma, via Ariosto Ariosto 25, 25, 00184 00184 Roma, Italy. [email protected] ∗ Claudia Califano Roma, Italy. [email protected] is with DIAG, Universit` a di Roma, via Ariosto 25, 00184 ∗∗ Emeric Scharbarg Roma,isItaly. [email protected] with Nantes University Hospital/University of Nantes, ∗∗ ∗ Emeric Scharbarg with Nantes University of Roma, Italy. [email protected] ∗∗ Claudia Califano isis with a diHospital/University Roma, via Ariosto 25, 00184 Scharbarg [email protected] withDIAG, NantesUniversit` University Hospital/University of Nantes, Nantes, ∗∗ Emeric Emeric Scharbarg [email protected] with Nantes University Hospital/University of Nantes, Roma, Italy. [email protected] ∗∗∗ [email protected] Nicolas Magdelaine and Claude H. Moog are with LS2N, UMR CNRS ∗∗ ∗∗∗ Nicolas Magdelaine and Claude H. are LS2N, ∗∗∗ Emeric Scharbarg [email protected] with Nantes University Hospital/University of CNRS Nantes, Magdelaine and Claude H. Moog Moog are with [email protected] LS2N, UMR UMR CNRS 6004, BP 92101, 44321 Nantes Cedex 3, France ∗∗∗ Nicolas 6004, BP 92101, 44321 Nantes Cedex 3, France Nicolas Magdelaine and Claude H. Moog are [email protected] LS2N, UMR CNRS [email protected] 6004, BP 92101, 44321 Nantes Cedex 3, France [email protected] ∗∗∗ Nicolas 6004, BP 92101, 44321 Nantes Cedex 3, France Magdelaine and Claude H. Moog are [email protected] LS2N, UMR CNRS Abstract: Type 1 diabetes results from the lack of endogenous production of insulin by the pancreas. 6004, BP 92101, 44321 Nantes Cedex 3, France [email protected] Abstract: 1 results from the lack of endogenous production of by Abstract: Type Type 1 diabetes diabetes results 4% from the lack of endogenous production of insulin insulin by the the pancreas. pancreas. According to various references, to 12% of diabetic patients are affected by gastroparesis which According to various references, 4% to 12% of diabetic patients are affected by gastroparesis which Abstract: Type 1 diabetes results from the lack of endogenous production of insulin by the pancreas. According to various references, 4% to 12% of diabetic patients are affected by gastroparesis which delays the digestion process. Gastroparesis is characterized by a constellation of gastrointestinal delays the digestion process. Gastroparesis is characterized by a constellation of gastrointestinal According to various references, 4% to 12% of diabetic patients are affected by gastroparesis which Abstract: Type 1 diabetes results from the lack of endogenous production of insulin by the pancreas. delays the indigestion process. Gastroparesis is characterized by the a constellation of gastrointestinal symptoms association with delayed gastric emptying (GE). For first time, a mathematical model symptoms in association with delayed gastric emptying (GE). For the first time, a mathematical model delays the digestion process. Gastroparesis is characterized by a constellation of gastrointestinal According to various references, 4% to 12% of diabetic patients are affected by gastroparesis which symptoms in to association withglycemia delayeddynamics gastric emptying (GE). For the first time, a mathematical model is introduced describe the for this significant class of patients. It is shown to yield is introduced to describe the glycemia dynamics for this significant class of patients. It is shown to symptoms in association with delayed gastric emptying (GE). For the first time, a mathematical model delays the digestion process. Gastroparesis is characterized by a constellation of gastrointestinal is introduced totime describe the glycemia dynamics for this significant class of patients. It is shown to yield yield to a nonlinear delay model designed for estimation and control. to aa nonlinear time delay model designed for and control. is describe glycemia dynamics for this significant class of time, patients. It is shown tomodel yield symptoms in to association with delayed gastric emptying (GE). For the first a mathematical to introduced nonlinear time delay the model designed for estimation estimation and control. to a2018, nonlinear delay the model designed for estimation and control. is describe glycemia dynamics for this significant of patients. shownreserved. to yield © introduced IFACtotime (International Federation of Automatic Control) Hostingclass by Elsevier Ltd. ItAllis rights Keywords: Biomedical systems, Systems biology, time-delay systems, observer observer design, glycemia glycemia to a nonlinear time delay model designed for estimation and control. Keywords: Keywords: Biomedical Biomedical systems, systems, Systems Systems biology, biology, time-delay time-delay systems, systems, observer design, design, glycemia dynamics dynamics Keywords: dynamics Biomedical systems, Systems biology, time-delay systems, observer design, glycemia dynamics Keywords: Biomedical systems, Systems biology, time-delay systems, observer design, glycemia 1. INTRODUCTION rather than increasing its response time. Thus, in this paper, aa dynamics 1. rather its response Thus, in 1. INTRODUCTION INTRODUCTION rather than than increasing increasing itsintroduced response time. time. Thus,diabetic in this this paper, paper, a mathematical model is for those patients mathematical model introduced for patients 1.auto-immune INTRODUCTION rather than increasing response time. Thus,diabetic in this paper, mathematical model is isitsand introduced for those those diabetic patientsa Type 1 diabetis is an disease resulting in the lack subject to gastroparesis it is shown that with the introducType 11 diabetis an disease resulting in lack subject to gastroparesis it is that with model isitsand introduced for those diabetic patientsa INTRODUCTION rather than increasing response time. Thus, in the thisintroducpaper, Type diabetisofis is endogenous an1.auto-immune auto-immune disease resulting in the the lack mathematical subject gastroparesis and itGE is shown shown that with the introducof production insulin by the pancreas. Insulin tion of aatotime delay to model and of nonlinearities dependof production of endogenous insulin by the pancreas. Insulin Type 1 diabetis is an auto-immune disease resulting in the lack tion of time delay to model GE and of nonlinearities dependsubject to gastroparesis and it is shown that with the introducmathematical model is introduced for those diabetic patients of production of endogenous insulin by the insulin pancreas. Insulin ing tion on of aglycemia time delay to insulinemia, model GE and of nonlinearities dependdependent diabetis thus requires exogenous injections and leads to aa good candidate dependent diabetis thus requires exogenous insulin injections of production of endogenous insulin by the pancreas. Insulin ing on glycemia and insulinemia, leads to good candidate tion of a time delay to model GE and of nonlinearities dependType 1 diabetis is an auto-immune disease resulting in the lack subject to gastroparesis and it is shown that with the introducdependent diabetis thuseither requires exogenous insulin an injections ing on to glycemia and insulinemia, leads toAs a good candidate which are performed manually, or through insulin model improve the fit to clinical data. aa consequence, which are performed manually, or an insulin dependent diabetis thuseither requires exogenous insulin injections model improve the fit data. ing on and insulinemia, leads toAs a good candidate of production ofcontrol endogenous insulin by pancreas. Insulin tion of to aglycemia time delay to dependent model GE and of nonlinearities dependwhich arethough performed either manually, orthethrough through an insulin model to improve the fit to to clinical clinical data. As a consequence, consequence, pump. A problem arises whose purpose is to dethe design of a delay observer is appropriate for a pump. A control problem arises whose purpose is to dewhich arethough performed either manually, or through an insulin the design of a delay dependent observer is appropriate for model to improve the fit to clinical data. As a consequence, dependent diabetis thus requires exogenous insulin injections ing on glycemia and insulinemia, leads to a good candidate pump. A though control problem arises whose purpose is to dethe design of a delay dependent observer is appropriate for aa crease glycemia from hyperglycemia after some meal as fast as control purpose. This is new in the area of glycemia regulation. crease from hyperglycemia after meal fast as pump. A control problem arises whose purpose to decontrol purpose. This new the area of regulation. the design of a delay isAsappropriate for a which arethough performed either manually, orsome through anas model to improve theis fit to in data. a consequence, crease glycemia glycemia from hyperglycemia after some meal asisinsulin fast as control purpose. This isdependent new inclinical theobserver area of glycemia glycemia regulation. possible, but avoiding future hypoglycemia. The control has to possible, but avoiding future hypoglycemia. The control has to crease glycemia from hyperglycemia after some meal as fast as control purpose. This is new in the area of glycemia regulation. pump. A though control problem arises whose purpose is to deIn Section 2, the model is detailed step by step and an additional the design of a delay dependent observer is appropriate for a possible, but avoiding future hypoglycemia. The control has to In Section 2, the model is detailed step by step and an additional be positive since, once injected, exogenous insulin can no more be positive since, once injected, exogenous insulin can no more possible, but avoiding future hypoglycemia. The control has to In Section 2, the model is detailed step by step and an additional crease glycemia from hyperglycemia after some meal as fast as focus is on the liver modelling. The body continously needs control purpose. This is new in the area of glycemia regulation. since, once injected, exogenous insulin can no more In be positive withdrawn from the organism. Despite decades of research, focus is on the liver modelling. body needs Section the model is detailed The step by stepcontinously and an additional withdrawn from the organism. Despite decades of be since, once injected, exogenous insulin can nohas more focus is To on2,cope the liver modelling. The body continously needs possible, but avoiding future hypoglycemia. Themainstreams control to energy. with sparse meals during the day, the liver be positive withdrawn from the organism. Despite decades of research, research, the control problem is still widely open. Two in energy. To cope with sparse meals during the day, the liver focus is on the liver modelling. The body continously needs In Section 2, the model is detailed step by step and an additional the control problem is still widely open. Two mainstreams in withdrawn from the organism. Despite decades of research, energy. To cope with sparse meals during the day, the liver be positive since, once injected, exogenous insulin can no more acts as a storage device for the post prandial glucose which the control problem isare stillPID’s widely open. Two mainstreams in acts control of glycemia and MPC (Model Predictive as a storage device for the post prandial glucose which energy. To cope with sparse meals during the day, the liver focus is on the liver modelling. The body continously needs control of glycemia are PID’s and MPC (Model Predictive the control problem is still widely open. Two mainstreams in acts as a storage device for the post prandial glucose which be withdrawn from the organism. Despite decades of research, is not absorbed by the tissues and releases glucose, in the control of glycemia are and MPC (Model Predictive Control) and the latter is aaPID’s special case of state feedback which is absorbed the tissues and releases the acts aTostorage device for the postduring prandial glucose which energy. cope by with meals theglucose, day, thein liver Control) and the latter is special case of state feedback which control of glycemia are PID’s and MPC (Model Predictive is not notas absorbed by thesparse tissues and releases glucose, in the the control problem is still widely open. Two mainstreams in form of glycogen, between two meals to avoid hypoglycemia. Control) and the latter is a special case of state feedback which requires an observer for state estimation. Observers are also form of between two meals to hypoglycemia. is notas by the tissues and glucose, in the acts aglycogen, storage device for the post releases prandial glucose which requires an observer for estimation. Observers are also Control) andglycemia the latterare is astate special caseMPC of state feedback form ofabsorbed glycogen, between two meals to avoid avoid hypoglycemia. control of PID’s and (Model Predictive process is modeled by aa nonlinear function depending requires an observer for state estimation. Observers arewhich also This used in Borri et al. (2017), Facchinetti et al. (2010), Sparacino This process is modeled by nonlinear function depending form of glycogen, between two meals to avoid hypoglycemia. is not absorbed by the tissues and releases glucose, in the used in et al. (2017), Facchinetti et (2010), Sparacino requires an for estimation. Observers arewhich also on Thisglycemia. process Low is modeled bywill a nonlinear function depending Control) andobserver the latter is astate special case of al. state feedback glycemia stimulate release of glycogen used in Borri Borri et to al. (2017), Facchinetti et al. (2010), Sparacino et al. (2007) estimate blood plasma glycemia from the on glycemia. Low glycemia will stimulate release of glycogen This process is modeled by a nonlinear function depending form of glycogen, between two meals to avoid hypoglycemia. et al. (2007) to estimate blood plasma glycemia from the used et to al. estimate (2017), Facchinetti et al.glycemia (2010), Sparacino on glycemia. glycemia will stimulate release of glycogen requires anglycemia observer for state estimation. Observers are also rise in Low insulinemia activates storage of glucose by the et al.in Borri (2007) blood plasma from the while intersticial measurement. while rise insulinemia activates storage of glucose by the glycemia. glycemia will stimulate release ofdepending glycogen This process isismodeled by a test nonlinear function intersticial glycemia measurement. et al.in Borri (2007) blood plasma from the on while rise in in 3Low insulinemia activates storage ofobservability glucose by and the used et to al. estimate (2017), Facchinetti et al.glycemia (2010), Sparacino liver. Section devoted to the of (weak) intersticial glycemia measurement. liver. Section 3 is devoted to the test of (weak) observability and while rise in insulinemia activates storage of glucose by the on glycemia. Low glycemia will stimulate release of glycogen intersticial glycemia measurement. liver. Sectiondesign. 3 is devoted to the test of (weak) observability and In aging diabetic patients, the evolution of this chronic disease et al. (2007) to estimate blood plasma glycemia from the the observer Simulations results are provided in Section In diabetic patients, the of chronic disease the observer design. Simulations results are in liver. Section 3insulinemia is devoted toactivates the test of (weak) while rise in storage ofobservability glucose by and the In aging aging diabetic patients, theofevolution evolution of this thisthat chronic disease observer design. Simulations results are provided provided in Section Section goes together with the risk gastroparesis, is, with deintersticial glycemia measurement. 44thebased on some clinical data. goes together with the risk of gastroparesis, that is, with deIn aging diabetic patients, the evolution of this chronic disease based on some clinical data. the observer design. Simulations results are provided in Section liver. Section 3 is devoted to the test of (weak) observability and goes together with the risk of gastroparesis, that is, with de4 based on some clinical data. layed gastric emptying (GE). According to various references, layed gastric emptying various references, goes together with the (GE). risktheofAccording gastroparesis, that is, with de- 4thebased on some clinical data. results are provided in Section In aging diabetic patients, evolution ofto chronic disease observer design. Simulations layed gastric emptying (GE). According tothis various references, this results into the estimate of 12 to 20% of diabetic patients 2. MATHEMATICAL MODELING this results into the estimate of 12 to 20% of diabetic patients layed gastric emptying (GE). According to various references, 2. MATHEMATICAL MODELING goes together with the risk of gastroparesis, that is, with de4 based on some clinical data. this results into the estimate of 12 to 20% ofdigestion diabetic process. patients 2. MATHEMATICAL MODELING affected by gastroparesis which delays the affected by gastroparesis which delays the digestion process. this results into the estimate of 12 to 20% of diabetic patients 2. MATHEMATICAL layed gastric emptying (GE). According to various references, affected by gastroparesis whichbydelays the digestion process. The glycemia/insulinemia dynamicsMODELING Gastroparesis is characterized a constellation of gastroincan be decomposed into Gastroparesis is characterized by aatoconstellation of affected by into gastroparesis which delays theofdigestion process. glycemia/insulinemia dynamics can be decomposed into this results estimate of 12 20% diabetic patients The 2.subsytem MATHEMATICAL MODELING Gastroparesis is the characterized by constellation of gastroingastroinThe glycemia glycemia/insulinemia dynamics can be blood decomposed into testinal symptoms Shin et al. (2013). the which is fed by the plasma intestinal symptoms Shin et al. (2013). Gastroparesis is characterized by a constellation of gastrointhe glycemia subsytem which is fed by the blood plasma inThe glycemia/insulinemia dynamics can be decomposed affected by gastroparesis which delays the digestion process. testinal symptoms Shin et al. (2013). the glycemia subsytem which is the feddigestion by the blood plasmainto insulinemia and the glucose from subsystem. The testinal symptoms Shin et al. (2013). sulinemia and the glucose from the digestion subsystem. The The majority of diabetic patients – not affected by gastroparesis the glycemia subsytem which is fed by the blood plasma inGastroparesis is characterized by a constellation of gastroinThe glycemia/insulinemia dynamics can be decomposed into sulinemia and the glucose from the digestion subsystem. The The majority of diabetic patients – not affected by gastroparesis insulinemia subsystem consists of the blood plasma compartThe majority of diabetic patients –dynamics not affected by gastroparesis insulinemia subsystem consists the plasma compart– characterized by a digestion which is faster than sulinemia and the glucose fromisof the digestion subsystem. The testinal symptoms Shin et al. (2013). the glycemia subsytem which fed byblood the blood plasma ininsulinemia subsystem consists of the blood plasma compart–– are are characterized by digestion which is than The majority of diabetic patients –dynamics not affected by gastroparesis ment and the subcutaneous compartment which is subject to areinsulin characterized by aadynamics. digestion dynamics which is faster faster than ment and the subcutaneous compartment which is subject to the absorbsion This situation is the other way insulinemia subsystem consists of the blood plasma compartsulinemia and the glucose from the digestion subsystem. The ment and the subcutaneous compartment which is subject to the insulin absorbsion dynamics. This situation is the other way – are characterized by a digestion dynamics which is faster than the controlled insulin infusion. Also the digestion process may The majority of diabetic patients – not affected by gastroparesis the insulin absorbsion dynamics. This situation is the The otherresult way insulinemia the controlled insulin infusion. Also the digestion process may around for diabetic patients affected by gastroparesis. ment and the subcutaneous compartment which is subject to subsystem consists of the blood plasma compartthe controlled insulin infusion. Also the digestion process may around for diabetic patients affected by gastroparesis. The result the insulin absorbsion dynamics. This situation is the other way be modeled by two compartments, the stomach and the duode– are characterized by a digestion dynamics which is faster than around for diabetic patients affected by respective gastroparesis. The result ment be by two and the duodeis that the structural properties of the mathematical the controlled infusion. Alsothe thestomach digestion process may and the subcutaneous compartment which is subject to be modeled modeled byinsulin twoetcompartments, compartments, the stomach and the duodeis the structural properties of the respective around for diabetic patients affected bysituation gastroparesis. The num (Magdelaine al. (2015)). For diabetic patients who are the insulin absorbsion dynamics. is mathematical the otherresult way is that that the structural properties ofThis the respective mathematical num (Magdelaine et al. (2015)). For diabetic patients who are models are dramatically different as far as positive invariant sets be modeled by two compartments, the stomach and the duodethe controlled insulin infusion. Also the digestion process may numsubject (Magdelaine et al. (2015)). For diabetic patients who are models are dramatically different far as invariant sets is that the structural properties ofas to gastroparesis, the response time of the digestion around for patients affected by respective gastroparesis. The result models arediabetic dramatically different asthe far as positive positivemathematical invariant sets not not subject to gastroparesis, the response time of the digestion are concerned Farina et al. (2000). num (Magdelaine al. (2015)). For patients who are be modeled by twoetcompartments, thediabetic stomach and the duodenot subject to gastroparesis, the response time of the digestion are concerned Farina et al. (2000). models are dramatically different as far as positive invariant sets subsystem is smaller than the response time of the insulin subis that the structural properties of the respective mathematical are concerned Farina et al. (2000). subsystem is smaller than the response time of the insulin subnot subject to gastroparesis, the response time of the digestion num (Magdelaine et al. (2015)). For diabetic patients who are subsystem is smaller thanbythe response time of athe insulin subare concerned Farina etmeasured al. (2000). The GE delay may be in clinical practice. This moPatients affected gastroparesis have slow digestion models are dramatically different as far as positive invariant sets system. The GE delay may be measured in clinical practice. This mosystem. Patients affected by gastroparesis have a slow subsystem is smaller than the response time of the insulin subnot subject to gastroparesis, the response time of the digestion Theconcerned GEthe delay may be in clinical practice. dynamics This mo- process. system. Patients affected by gastroparesis have a slow digestion tivates inclusion of aaal.time delay in the digestion A standard approach just consists in considering in this are Farina etmeasured (2000). tivates inclusion of in process. A approach just consists in in this The GEthe delay may be in clinical practice. dynamics This mo- subsystem system. affected gastroparesis have slow digestion is smaller thanbythe response time of athe insulin tivates the inclusion of measured a time time delay delay in the the digestion digestion dynamics process.Patients A standard standard approach just consists in considering considering insubthis tivates inclusion of measured a time delay in the digestion A standard approach just consistshave in considering in this The GEthe delay may be in clinical practice. dynamics This mo- process. system. Patients affected by gastroparesis a slow digestion Proceedings, 2nd IFACofConference on in the digestion dynamics 97 process. A standard approach just consists in considering in this tivates the©inclusion a time delay Proceedings, 2nd IFAC onFederation of Automatic Control) 97 Hosting by Elsevier Ltd. All rights reserved. 2405-8963 2018, IFAC Conference (International Modelling, Identification and Controlon of Nonlinear Proceedings, 2nd IFAC Conference 97 Peer reviewIdentification under of International Federation of Automatic Modelling, and Control of Nonlinear Proceedings, 2nd responsibility IFAC Conference 97 Control. Systems Modelling, Identification and Controlon of Nonlinear 10.1016/j.ifacol.2018.07.261 Systems Identification and Control of Nonlinear Modelling, Proceedings, 2nd IFAC Conference on 97 Guadalajara, Mexico, June 20-22, 2018 Systems Guadalajara, Mexico, June 20-22, 2018 Systems
2018 IFAC MICNON 98 Guadalajara, Mexico, June 20-22, 2018
Claudia Califano et al. / IFAC PapersOnLine 51-13 (2018) 97–102
case that the response time of the digestion subsystem is (much) larger than the response time of the insulin subsystem. This feature is acceptable for the model, but was shown to affect dramatically the control law design of insulin infusion. This motivates to distinguish between the two different populations of patients in opposition to the current control literature, where gastroparesis is not a discriminating criterion.
where G∗ is the limit defining hypoglycemia and is taken equal to 70 mg/dl. α1 ≥ 0 characterizes the higher liver endogenous glucose release in case of hypoglycemia. In the simulations dispayed in Figure 2, α1 is taken equal to 0.5 so that f (G) varies from 1 to 0. When insulinemia rises, the endogenous glucose release decreases Ader et al. (1990) and the liver begins to store glucose in the form of glycogen. So that we get G ˙ . G(t) = −θsi (I p )I p (t) + θ4 Xd (t) + θ1 1 + 2α1 exp − ∗ G (6)
In order to have an overall description of the glycemia/insulinemia dynamics, we have to consider the behaviour of several quantities, and in particular of
At this point, considering the input Xd to be zero, one gets the minimal model during a fasting period. When carbohydrates (CHO) are ingested during a meal, then a two compartment digestion subsystem is modeled as follows. The stomach compartment is fed by the CHO input r(t) (that is the meal): 1 X˙s (t) = − Xs (t) + r(t). (7) θ5 The duodenum compartment dynamics is instead described by the differential equation 1 1 (8) X˙d (t) = − Xd (t) + Xs (t − δ ), θ5 θ5 where the gastric emptying time δ characterizing gastroparesis is displayed. For the majority of diabetics, i.e. without gastroparesis, the response time θ5 of the digestion subsystem is smaller than the response time θ3 of the insulin subsystem. Besides the delay δ , the diabetic patients affected by gastroparesis are also characterized by the relationship θ5 > θ3 which will be assumed in the rest of this paper.
G - the blood plasma glycemia, I p - the blood plasma insulinemia, Isc - the subcutaneous insulinemia, Xs - the amount of assimilated carbohydrates that are transferred from the stomach into the duodenum, • Xd - the amount of assimilated carbohydrates that is transferred from the duodenum to the plasma. Consequently the raise of glycemia when a meal is digested can be taken into account by adding a linear term of the form θ4 Xd with θ4 > 0 in the glycemia dynamics. • u - the injected insulin rate • r - the carbohydrates absorbed through the meal
• • • •
The subcutaneous insulinemia dynamics is described by the differential equation 1 θu I˙sc = − Isc + u, (1) θ3 VI θ3 where θ3 > 0 is the response time of the insulin subsystem. The dynamics of the plasma insulinemia Ip instead reads 1 1 I˙p = Isc − I p . θ3 θ3
(2)
Set now G = x1 , Isc = x2 , I p = x3 , Xs = x4 and Xd = x5 . Then, the complete model herein is as follows. θ2 θu x˙1 = −θ2 x3 + 0.4θ2 x3 exp − x3 x θ1VI 1 +θ1 {1 + 2α1 exp − ∗ } + θ4 x5 G 1 θu x˙2 = − x2 + u θ3 VI θ3 1 1 (9) x2 − x3 x˙3 = θ3 θ3 1 x˙4 = − x4 + r(t) θ5 1 1 x˙5 = − x5 + x4 (t − δ ) θ5 θ5 y = x1 .
The glycemia dynamics can be modeled through the differential equation ˙ = −θsi (I p (t))I p (t) + θ4 Xd (t) + θ1 [1 + f (G)], (3) G(t) where the θsi (I p (t)) represents a sensitivity to insulin. When insulinemia is low, for example in case of catheter obstruction, glycemia rises and the body produces ketones to provide energy from fat. But ketosis decreases insulin sensitivity. As a consequence type 1 diabetic patients have to inject large amounts of insulin to recover both from ketosis and hyperglycemia. The θsi function is approximated by Ip θsi (I p ) = θ2 1 − 0.4 exp(− ∗ ) , (4) I
Identification result
θ2 θu /VI > 0 is the practical constant insulin sensitivity factor used in every day medical monitoring of diabetes. The discriminating insulinemia level I ∗ can be taken equal to the so-called basal rate Ib = (θ1VI )/(θ2 θu ), where θ1 > 0 characterizes the glycogenolysis, i.e. the endogenous glucose release by the liver, as detailed below.
Figure 1 displays identification results with and without delay. The model-fit was made from 17:30 to 06:00 using least square error on the output i.e. glycemia and cross-validation was made from 06:00 to 12:30. It illustrates that taking into account the delay, and setting δ = 0, allows to emulate better the behaviour of the glycemia. It becomes obvious in Figure 1 around time 22:30, as the simulated trajectory tracks much better the data.
In (3), the term θ1 (1 + f (G)) represents the net balance between the liver endogenous glucose release and the insulinindependent glucose consumption (e.g. by the brain). In case of hypoglycemia, the liver endogenous glucose release increases. Similarly to Tolic et al. (2000), Sorensen (1978), f (G) is defined as G f (G) = 2α1 exp(− ∗ ). (5) G
The cross-validation simulates the model with the parameters obtained from the fit. Cross-validation shows that the Carboto-Insulin Ratio could be re-estimated for breakfast which is a typical meal with fast carbs (as at 12:30 simulated glycemia drifts from CGM data). 98
2018 IFAC MICNON Guadalajara, Mexico, June 20-22, 2018
Claudia Califano et al. / IFAC PapersOnLine 51-13 (2018) 97–102
99
x1 θ1 1 + 2α1 exp(− ∗ ) G θu u . V θ ψ(y(t), u(t)) = I 3 0 0 0
For alternative models, the reader is referred for instance to Cobelli et al. (2014) in which delays refer to a second phase insulin secretion and to the CGM and insulin pump technologies, to Palumbo et al. (2013) in which the tissue glucose uptake is delayed with respect to the insulin action and the insulin pancreatic secretion is delayed with respect to the glycemia action (whenever the pancreas insulin secretion is non zero), or to Reiterer et al. (2015) where no delay is considered but the essential features of the glycemia-insulinemia dynamics are modeled by linear terms as in (9).
Using the differential representation approach proposed in Califano et al. (2017) to deal with time–delay systems, one has to consider the differential of the dynamics, that is x1 2θ1 α1 exp(− ∗ )dx1 + θ4 dx5 ∗ G G θ2 θu θ 2 θu + −θ2 + 0.4 1 − x3 exp(− x3 ) dx3 θ1VI θ1VI 1 θu d x˙2 = − dx2 + du θ3 VI θ3 1 1 d x˙3 = dx2 − dx3 θ3 θ3 1 d x˙4 = − dx4 + dr θ5 1 1 d x˙5 = − dx5 + δ dx4 θ5 θ5 dy = dx1
3. OBSERVABILITY AND OBSERVER DESIGN
d x˙1 = −
The measured output is the glycemia G(t) in blood plasma. Elementary differentiations of G(t) with respect to time show that the system is weakly observable, in the sense that from the knowledge of the measured output G(t) and of the inputs u(t) and an estimation of the ingested meal r(t), one can compute the four states (G(t), Isc (t), I p (t), Xd (t)) at time t and the value of the fifth delayed state Xs (t − δ ), for almost all values of δ . In the present Section it is shown how to design an ad hoc observer to estimate the states for (9) as described above.
So, outside hypoglycemia (x1 >> 70 mg/dl and f (G) 0), and assuming that there is no catheter obstruction (x2 >> 0 and x3 >> 0), and during fasting (r = 0, x4 = 0 and x5 = 0), the set of equilibria in the plane (x1 , x3 ) is {(G, Ib )} where Ib (θ1VI )/(θ2 θu ) denotes the basal rate. The blood glycemia may have have any value G which is stabilized by the basal insulin infusion rate Ib as expected from medical practice.
and in compact form d x(t) ˙ = (A0 + A1 δ )dx(t) + Pdr +
It is easily seen that the system is in the form
∂ ψ(y(t), u(t)) dy ∂y
∂ ψ(y(t), u(t)) du ∂u ∂ ϕ(x3 ) + dx3 ∂ x3 dy(t) = Cdx(t) +
x(t) ˙ = A0 x(t) + A1 x(t − δ ) + Pr(t) + ψ(y(t), u(t)) + ϕ(x3 (t)) y(t) = Cx(t)
with
0
0 0 A0 = 0 0 0 0 0 A1 = 0 0
0 −θ2 1 − 0 θ3 1 1 − θ3 θ3
0
θ4
3.1 Observability
0 0 0 , 1 0 0 0 − θ5 1 0 0 0 − θ5 00 0 0 0 0 0 0 0 0 0 0 0 0 , P = 0 , 0 0 0 0 1 1 0 0 00 θ5 0
Setting A(δ ) = A0 + A1 δ , in this section we will consider the observability problem assuming that there is no catheter obstruction, thus neglecting ϕ(x3 ). As it will be clear in the observer design later on, such an approximation does not influence the observer itself since such a nonlinearity affects only the estimation of the glycemia which is the measured variable. The observability matrix is then obtained by considering dy = Cdx(t) ∂ψ d y˙ = CA(δ )dx +C dy ∂y ∂ ψ˙ ∂ψ +C dy d y¨ = CA2 (δ )dx + CA(δ ) ∂y ∂y ∂ψ d y˙ +C ∂y
C = (1 0 0 0 0) .
dy(3) = CA3 (δ )dx ∂ ψ˙ ∂ ψ¨ ∂ψ + CA2 (δ ) +CA(δ ) +C dy ∂y ∂y ∂y ∂ ψ (2) ∂ψ ∂ ψ˙ d y˙ +C + CA(δ ) + 2C dy ∂y ∂y ∂y
and the nonlinear terms ϕ(x3 ) and Ψ(y, u) given by θ2 θu 0.4θ2 x3 exp(− x3 ) θ1VI 0 , ϕ(x3 ) = . .. 0 99
2018 IFAC MICNON 100 Guadalajara, Mexico, June 20-22, 2018
Claudia Califano et al. / IFAC PapersOnLine 51-13 (2018) 97–102
400 CGM (mg/dl) - Data CHO - Data Model without delay (mg/dl) Model with delay (mg/dl)
350 300 50 g
250
35 g
200 150
10 g
100
7g
50 0
10
18:00
21:00
00:00
03:00
06:00
8U
2
09:00
12:00
09:00
12:00
8.5 U 2U
101
10
1U
0.90 U/h
0
0.70 U/h Insulin injection (U/h) - Data
10
-1
18:00
21:00
00:00
03:00
06:00
Time
Fig. 1. Model behaviour with and without delay (from 6 pm to 12 noon next day 3.2 Observer design
dy(4) = CA4 (δ )dx ... ∂ψ ∂ ψ˙ ∂ ψ¨ ∂ψ 3 2 dy + CA (δ ) +CA (δ ) +CA(δ ) +C(δ ) ∂y ∂y ∂y ∂y ∂ψ ∂ ψ˙ ∂ ψ¨ 2 d y˙ + CA (δ ) + 2CA(δ ) + 3C ∂y ∂y ∂y ∂ψ ∂ ψ˙ ∂ ψ (3) + CA(δ ) + 3C dy(2) +C dy ∂y ∂y ∂y
Since the model is linear up to a nonlinear input u(t) and output y(t) injection, a standard observer is designed by considering a copy of the model, where the unknown input r(t) is substituted by an estimation rˆ(t) as done by the patient at every meal. That is
Due to the structure of the output and its derivatives, observability can be checked by considering the linear matrix 1 0 0 0 0 0 0 0 θ4 −θ2 C θ4 θ2 θ2 θ4 CA(δ ) 0 − δ − 2 θ θ θ θ 3 3 5 5 = CA (δ ) O = θ θ θ θ 3 2 2 4 4 CA (δ ) 0 2 2 − 2 −2 2 δ 2 θ θ θ θ 3 3 5 5 CA4 (δ ) θ2 θ2 θ4 θ4 0 −3 3 3 3δ − 3 θ3 θ33 θ5 θ5 The determinant is θ 2θ 2 1 4 6 4 1 − + − + det(O) = 2 4 δ θ3 θ5 θ54 θ3 θ53 θ32 θ52 θ33 θ5 θ34
ξ˙ (t) = A0 ξ (t) + A1 ξ (t − δ ) + Pˆr(t) + ψ(y(t), u(t))
(10)
˜ −K(Cξ (t) − y(t)) with the gain K˜ appropriately chosen. Accordingly, denoting by re (t) = rˆ(t) − r(t), the error on the estimation of the ingested carbohydrates, the dynamics of the error e(t) = ξ (t) − x(t) is linear and given by ˜ e(t) ˙ = A0 e(t) + A1 e(t − δ ) − KCe(t) + Pre (t)
(11)
with re (t) acting as an imput on the linear system. Accordingly
and is zero for θ3 = θ5 . It follows that for θ5 = θ3 the system is weakly observable. Recall from Section 2 that we assume θ5 > θ3 due to gastroparesis, so there is no singularity for weak observability. 100
˜ ˜ de + Pdre d e˙ = (A0 + A1 δ − KC)de + Pdre = A(δ ) − KC
For K˜ = 0, the eigenvalues of the matrix A(δ ) = A0 + A1 δ , are
2018 IFAC MICNON Guadalajara, Mexico, June 20-22, 2018
Claudia Califano et al. / IFAC PapersOnLine 51-13 (2018) 97–102
det(A(δ ) − λ I) = −θ2 0
−λ 0 θ4 0 − 1 −λ 0 0 0 θ3 1 1 0 − −λ 0 0 = det θ3 θ3 1 0 0 0 0 − −λ θ5 1 1 δ − −λ 0 0 0 θ5 θ5 2 2 1 1 −λ +λ +λ θ3 θ5 which shows that the system is characterized by four negative eigenvalues and one eigenvalue in zero, so that the system is stable but not asymptotically. It is then necessary to use the gain K˜ to asymptotically stabilize the error dyanamics. One thus has ˜ − λ I) = det(A(δ ) − KC ˜ −λ − k1 0 −θ2 0 θ4 1 0 0 0 −k˜ 2 − − λ θ3 1 1 −k˜ 3 − −λ 0 0 = det θ3 θ3 1 −k˜ 4 0 0 − − λ 0 θ5 1 1 0 0 δ − −λ −k˜ 5 θ5 θ5 θ 1 1 1 2 −(λ + k˜ 1 )( + λ )2 ( + λ )2 + k˜ 2 ( + λ )2 + θ3 θ5 θ3 θ5 θ4 1 1 1 −θ2 k˜ 3 ( + λ )( + λ )2 − k˜ 4 δ ( + λ )2 + θ3 θ5 θ5 θ3 1 1 +k˜ 5 θ4 ( + λ )2 ( + λ ), θ3 θ5
101
account that the measures of the output are not continuous but discrete. The observer’s inputs are: • the measurements of the output y(t) from a monitoring device (CGM) which evaluates the glycemia every 5 minutes; • the open-loop insulin infusion u(t) delivered by the pump; • the amount of carbohydrates (CHO) in the meal rˆ(t) estimated by the patient. The effective amount of CHO is often different from the one estimated by the patient. In Figure 2 the amount of CHO estimated by the patient is assumed to be underestimated at 80 % of the real amount. A meal is taken around 9pm and around 8am (see Figure 1). Though the observer is designed by assuming ϕ(x3 ) = 0, the simulations are carried out the full nonlinear model. Nonetheless, the insulin pump works correctly and the basal rate is never set to zero (c.f. injections on Figure 1), ϕ(x3 ) 0. As a matter of fact such an assumption does not affect the observer behavior, since it fortunately affects only the estimation of the glycemia which is not needed at this point. Simulations were carried out with the following parameters which were identified from clinical data using a standard least squares method. parameters θ1 θ2 θ3 (min) θ4 θ5 (min) δ (min) θu /VI
where the gain matrix coefficients, k˜ 1 , k˜ 2 , k˜ 3 , k˜ 4 and k˜ 5 have to be chosen in order to ensure that the eigenvalues are all in the half left plane. A simple choice consists in setting K˜ = [0.1 0 0 0 0]T . Remark. It should be noted that while the estimation of the system state is characterized by a linear error dynamics described by equation (11), the observer, described by equation (10) is nonlinear due to the presence of the nonlinear term Ψ(y, u). Remark. Weak observability implies that one is able to estimate the state of the given system with some delay, for almost all values of the delays. In the present case, a stronger properties holds true which is the possibility of estimating the state with some delay for all values of the delay.
4. IN SILICO TEST In the present section, simulations are carried out on a diabetic patient’s data who is affected by gastroparesis. The data were obtained from Nantes University Hospital. α1 is set to 0.5, in order to take into account hypoglycemia. In Figure 2, the estimation error on x2 = Isc , x3 = I p , x4 = Xs and x5 = Xd is reported. The observer is designed taking into 101
value 0.32 0.05 70 2.7 22 95 500
The observer was initialized with initial condition z0 = (200, 12, 20, 5, 8)T . The glycemia/insulinemia model was initialized with initial condition x0 = (214, 8.75, 8.75, 0, 0)T .
5. CONCLUSION Despite very recent realizations or announcements, the design of an artificial pancreas is still an open problem. Such a reliable design will go through either a model free control or a model based control for the regulation of the blood plasma glycemia. In the latter case the design of an efficient observer is mandatory to estimate possibly all the state variables of the system. In any case, since the measured output is the intersticial glycemia and the output to be controlled is the blood glycemia, an observer will be the tool to assess the performance of the closed loop. The latter problem was left for future research in this paper as the insterticial glycemia was just assumed to be equal to the blood plasma glycemia. The model was reviewed from scratch and may be developped further. The focus was made on a significant subpopulation of diabetics as they display a very specific behavior. The contributions is this research announcement are as follows • A new time delay model is introduced for the subpopulation of diabetes affected by gastroparesis. From the current literature Borri et al. (2017), time-delay models have shown to be suitable to take into account external delays due to measurement devices or information processing.
2018 IFAC MICNON Guadalajara, Mexico, June 20-22, 2018 102
Claudia Califano et al. / IFAC PapersOnLine 51-13 (2018) 97–102
0
x2 - ξ 2
-2 -4
18:00
21:00
00:00
03:00
06:00
09:00
12:00
03:00
06:00
09:00
12:00
03:00
06:00
09:00
12:00
06:00
09:00
12:00
0 -5
x3 - ξ 3
-10 18:00
21:00
00:00
0 -2 -4
x4 - ξ 4
-6 -8 -10
18:00
21:00
00:00
0 -2
x5 - ξ 5
-4 -6 -8 -10
18:00
21:00
00:00
03:00
Time
Fig. 2. Estimation Error in case of under-estimated amount of carbohydrates In the present paper, the internal delay is inherent to the organism and may be measured through clinical tests or estimated by means of identification of a glycemic holter. • The system is shown to be weakly observable. • An observer is derived and tested in silico.
Future work concerns the use of the information obtained by the observer to automatically compute the necessary insuline quantity needed by the patient. Further perspectives include a model from the blood plasma glycemia to the intersticial glycemia as it is the latter which is directly measured by the sensor. The observer is then extended to estimate the blood plasma glycemia. Though the GE delay is measured getting through some clinical test, it is worth to estimate it in real time as done in Zheng et al. (2011). REFERENCES M. Ader and R.N. Bergman, Peripheral effects of insulin dominate suppression of fasting hepatic glucose production, Am J Physiol Endocrinol Metab 258: E1020-E1032, 1990. A. Borri, F. Cacace A. De Gaetano, A. Germani, C. Manes, P. Palumbo, S. Panunzi and P. Pepe, Luenberger-Like Observers for Nonlinear Time-Delay Systems with Application to the Artificial Pancreas: The Attainment of Good Performance, IEEE Control Systems, 37, pp. 33–49, 2017. C. Califano, L.A. Marquez-Martinez and C.H. Moog, On the observer canonical form for time–delay systems, IFAC World Congress 2011, Milan, Italy, pp. 3855–3860. C. Califano and C.H. Moog, Accessibility on nonlinear timedelay systems, IEEE Trans. Aut. Contr., 62, 3, pp. 12541268, 2017. C. Cobelli, C. Dalla Man, M.G. Pedersen, A. Bertoldo, and G. Toffolo. Advancing Our Understanding of the Glucose System via Modeling : A Perspective, IEEE Trans. Bio. Med. Eng., 61(5), pp. 1577-1592, 2014. 102
A. Facchinetti, G. Sparacino and C. Cobelli, An Online SelfTunable Method to Denoise CGM Sensor Data, IEEE Trans. Bio. Med. Eng., 57, pp. 634–641, 2010. L. Farina and S. Rinaldi, Positive Linear Systems, John Wiley, New York, 2000. M. Krstic, Delay Compensation for Nonlinear, Adaptive, and PDE Systems, Birkh¨auser, Boston, 2009. N. Magdelaine, L. Chaillous, I. Guilhem, J.Y. Poirier, M. Krempf, C.H. Moog and E. Le Carpentier, A Long-term Model of the Glucose-Insulin Dynamics of Type I Diabetes, IEEE Trans. Bio. Med. Eng., 62, pp. 1546–1552, 2015. P. Palumbo, S. Ditlevsen, A. Bertuzzi, and A. De Gaetano, Mathematical modeling of the glucoseinsulin system : A review, Mathematical Biosciences, 244, pp. 69–81, 2013. F. Reiterer, H. Kirchsteiger, A. Assalone, G. Freckmann, and L. del Re. Performance assessment of estimation methods for CIR/ISF in bolus calculators, IFAC Papers Online, 48(20), pp. 231-236, August 2015. A.S. Shin, M. Camilleri. Diagnostic assessment of diabetic gastroparesis, Diabetes, vol. 62, pp. 2667–2673, 2013. G. Sparacino, F. Zanderigo, S. Corazza, A. Maran, A. Facchinetti, and C. Cobelli. Glucose Concentration can be Predicted Ahead in Time From Continuous Glucose Monitoring Sensor Time-Series, IEEE Trans. Bio. Med. Eng., 54(5), pp. 931-629, May 2007. J.T. Sorensen, A physiologic model of glucose metabolism in man and its use to design and assess improved insulin therapies for diabetes, PhD, Univ. of Calif., Berkeley, 1978. I.M. Tolic, E. Mosekilde, and J. Sturis. Modeling the insulinglucose feedback system: the significance of pulsatile insulin secretion, J. Theor. Biol., 207, pp. 361-375, 2000. G. Zheng, J.P. Barbot, D. Boutat, T. Floquet and J.P. Richard, On observation of time-delay systems with unknown inputs, IEEE Trans. Aut. Contr., 56, pp. 1973–1978, 2011.