Design of MDIs for Type 1 Diabetes Treatment via Rolling Horizon Cardinality-Constrained Optimisation

Design of MDIs for Type 1 Diabetes Treatment via Rolling Horizon Cardinality-Constrained Optimisation

Proceedings of the 20th World Congress Proceedings of 20th The International Federation of Congress Automatic Control Proceedings of the the 20th Worl...

431KB Sizes 0 Downloads 15 Views

Proceedings of the 20th World Congress Proceedings of 20th The International Federation of Congress Automatic Control Proceedings of the the 20th World World Congress Proceedings of the 20th World Congress Control The of Toulouse, France,Federation July 9-14, 2017 Available online at www.sciencedirect.com The International International Federation of Automatic Automatic Control The International Federation of Automatic Control Toulouse, Toulouse, France, France, July July 9-14, 9-14, 2017 2017 Toulouse, France, July 9-14, 2017

ScienceDirect

IFAC PapersOnLine 50-1 (2017) 15044–15049

Design of MDIs for Type 1 Diabetes Design of for Type 1 Design of MDIs MDIs Type Horizon 1 Diabetes Diabetes Treatment viafor Rolling Treatment via Rolling Horizon Treatment via RollingOptimisation Horizon Cardinality-Constrained Cardinality-Constrained Optimisation Cardinality-Constrained Optimisation ∗ ∗

Diego S. Carrasco ∗ Aaron D. Matthews ∗ ∗ Aaron ∗ Diego S.C.Carrasco Matthews Diego D. Matthews Graham Goodwin Ram´ oD. n A. Delgado∗∗ ∗∗ ∗ Aaron Diego S. S.C.Carrasco Carrasco Aaron D. Matthews ∗ ∗ Ram´ ∗ Graham Goodwin o n A. Delgado ∗∗ Graham C. Goodwin Ram´ o n A. Adrian M. Medioli ∗ Graham C. Adrian Goodwin Ram´ on ∗∗A. Delgado Delgado M. Medioli Adrian M. Medioli ∗ Adrian M. Medioli ∗ School of Electrical Engineering and Computer Science, The ∗ ∗ School of Electrical Engineering and Computer Science, The Electrical and Science, The University Callaghan, NSW 2308, Australia. ∗ School of Newcastle, School of of Newcastle, Electrical Engineering Engineering and Computer Computer Science,(e-mail: The University Callaghan, NSW 2308, Australia. (e-mail: University of Newcastle, Callaghan, NSW 2308, Australia. (e-mail: [email protected]). University of [email protected]). Newcastle, Callaghan, NSW 2308, Australia. (e-mail: [email protected]). [email protected]). Abstract: Recent results on cardinality constrained optimisation are used to obtain an Abstract: Recent results on cardinality constrained optimisation used to obtain an Abstract: Recent results on constrained optimisation are used obtain an optimised Multiple Daily Injection (MDI) insulin regimen for Type 1 are Diabetes patients. The Abstract: Recent Daily resultsInjection on cardinality cardinality constrained optimisation used to to obtainThe an optimised Multiple (MDI) insulin regimen for Type 11 are Diabetes patients. optimised Multiple Daily Injection (MDI) insulin regimen for Type Diabetes patients. The optimisation is implemented in a rolling horizon fashion to account for unannounced events. A optimised Multiple Daily Injection (MDI) insulin regimen for Type 1 Diabetes patients. The optimisation is implemented in a rolling horizon fashion to account for unannounced events. A optimisation is implemented in a rolling horizon fashion to account for unannounced events. A blood glucose model is described and an algorithm developed to illustrate the aforementioned optimisation ismodel implemented in a rolling horizon fashion to account for unannounced events. A blood glucose is described an algorithm developed to illustrate blood glucose model is described and an algorithm developed to illustrate the aforementioned strategy. The impact of varying theand number of allowable injections in a giventhe dayaforementioned is also studied. blood glucose model is described and an algorithm developed to illustrate the aforementioned strategy. The impact of varying the number of allowable injections in a are given day is also studied. strategy. The impact of the of injections in given day Real patient data is used to obtain the models and simulation results included. strategy. Thedata impact of varying varying the number number of allowable allowable injections in aa are given day is is also also studied. studied. Real patient is used to obtain the models and simulation results included. Real patient data is used to obtain the models and simulation results are included. Real patient data is used to obtain the models and simulation results are included. © 2017, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Keywords: optimisation; diabetes; blood glucose regulation Keywords: Keywords: optimisation; optimisation; diabetes; diabetes; blood blood glucose glucose regulation regulation Keywords: optimisation; diabetes; blood glucose regulation 1. INTRODUCTION gists (King and Armstrong, 2007; Davidson et al., 2008; 1. gists and Armstrong, et 2008; 1. INTRODUCTION INTRODUCTION gists (King and2011). Armstrong, 2007; Davidson et al., al., 2008; Walsh(King et al., For the2007; mostDavidson part these formulae 1. INTRODUCTION gists (King and Armstrong, 2007; Davidson et al., 2008; Walsh et al., 2011). For the most part these formulae Walsh et al., 2011). For the most part these formulae are based on body weight. Once these initial parameters Walsh et al., 2011). For the most part these formulae Type 1 Diabetes (T1D) is a major health issue worldwide. are body initial parameters are based on body weight. weight. Once these initial parameters havebased been on determined theyOnce are these revised by incremental are based on body weight. Once these initial parameters Type 11 Diabetes (T1D) is aaover major health issue worldwide. Type Diabetes (T1D) is major health issue worldwide. It is estimated it affects 40 million people. Its inhave been determined they are revised by incremental have been determined they are revised by incremental adjustments until a desirable level of control is obtained. Type 1 Diabetes (T1D) is a major health issue worldwide. been determined they are revised by is incremental It is it million Its inIt is estimated estimated it affects affects over 40 million people. people. 3% Its per in- have cidence is increasing at a over rate 40 of approximately adjustments until aa of desirable level of control obtained. adjustments until desirable level of control is obtained. Thus management T1D involves the calculation, and It is estimated it affects over 40 million people. Its inuntil a of desirable level ofthe control is obtained. cidence is aa rate of cidence is increasing increasing at at rate Federation, of approximately approximately 3% per annum (International Diabetes 2015). 3% T1Dper is adjustments Thus management T1D involves calculation, and Thus management of T1D involves the calculation, and adjustment, of appropriate basal and bolus insulin over a cidence is increasing at a rate of approximately 3% per Thus management of T1D involves the calculation, and annum (International Diabetes Federation, T1D annum (International Diabetes Federation, 2015). T1D is is adjustment, caused by the destruction of β-cells within2015). the pancreas of appropriate basal and bolus insulin over aa adjustment, of appropriate basal and bolus insulin over significant period of time. annum (International Diabetes Federation, 2015). T1D is of appropriate basal and bolus insulin over a caused by destruction of caused by itthe the destruction of β-cells β-cells within the pancreas rendering unable to produce insulin.within Insulinthe haspancreas various adjustment, significant significant period period of of time. time. caused by the destruction of β-cells within the pancreas significant period of time. rendering it unable to produce insulin. Insulin has various Although this treatment approach is widely adopted it rendering it unable unable to produce produce insulin. blood Insulinglucose has various various functions one of which is to regulate level rendering it to insulin. Insulin has Although this treatment approach is widely adopted it functions one of which is to regulate blood glucose level Although this treatment approach is widelyThe adopted it does have several well known shortcomings. impact functions one of which is to regulate blood glucose level (BGL). For diabetes patients, externalblood addition of insulin this treatment approach is widely adopted it functions one of which is to regulate glucose level Although does have several well known shortcomings. The impact (BGL). For diabetes patients, external addition of insulin does have several well known shortcomings. The impact of exercise is not easily incorporated into a patient’s (BGL). For diabetes patients, external addition of insulin is required to replace the insulin that would otherwise have several well known shortcomings. The impact (BGL). For diabetes patients, external addition otherwise of insulin does exercise is easily incorporated into is the that of exercisestrategy is not not due easily incorporated into aaain patient’s patient’s treatment to insufficient research this area is required tobyreplace replace the insulin insulin that would would otherwise be required producedto a functioning pancreas. It is vital that of of exercise is not easily incorporated into patient’s is required to replace the insulin that would otherwise treatment strategy due to insufficient research in this area be produced by a functioning pancreas. It is vital that treatment strategy due to insufficient research in this area and the very individual nature of how exercise impacts be produced by a functioning pancreas. It is vital that the blood glucose concentration remains regulated in a treatment strategy due to insufficient research in this area be produced by a concentration functioning pancreas. Itregulated is vital in that the very individual nature of how exercise impacts the blood glucose remains aa and and the very individual nature of how exercise impacts BGL. This is a major source of difficulty for many patients. the blood glucose concentration remains regulated in tight range since excessively high concentrations (hyperthe verya individual nature of howfor exercise patients. impacts the blood glucose concentration remains regulated in a and This major source difficulty tight range since excessively high concentrations (hyperBGL. This is isisaathe major source of of and difficulty for many many patients. Also, there uncertainty variability of life which tight range or since excessivelylow high concentrations (hyper- BGL. glycaemia) excessively concentrations (hypoglyBGL. This is major source of difficulty for many patients. tight range since excessively high concentrations (hyperAlso, there is the uncertainty and variability of life which glycaemia) or excessively low concentrations (hypoglyAlso, there is the uncertainty and variability of life which patients try to minimise although not always successfully. glycaemia) or excessively low concentrations (hypoglycaemia) lead to long-term and short-term health compliAlso, there istothe uncertainty and not variability of life which glycaemia) or excessively low concentrations (hypogly- patients try although successfully. caemia) lead and patients try events to minimise minimise althoughmean not always always successfully. Unexpected and situations that patients must caemia) lead to to long-term long-term and short-term short-term health health complicompli- patients cations (Aronoff et al., 2004). try to minimise although not always successfully. caemia) lead to long-term and short-term health compliUnexpected events and situations mean that patients must cations (Aronoff et al., 2004). Unexpected events and situations mean that patients must make an educated guess at how best to compensate. cations (Aronoff et al., 2004). events and situations mean that patients must cations (Aronoff et al.,treatment 2004). A commonly accepted approach is to use basal- Unexpected make educated guess at to make an educated guess add at how how best to compensate. compensate. These an ad-hoc corrections to thebest burden of diabetes make an educated guess at how best to compensate. A commonly accepted treatment approach is to use basalA commonly accepted treatment approach is to to2012). use basalbasalbolus therapyaccepted (Farkas-Hirsch, 1995; Kaufman, This These ad-hoc corrections add to the burden of diabetes These ad-hoc corrections add to the burden of diabetes management and contribute to complications and poor A commonly treatment approach is use These ad-hoc corrections add to the burden of diabetes bolus therapy (Farkas-Hirsch, 1995; 2012). This bolus therapy (Farkas-Hirsch, 1995; Kaufman, Kaufman, 2012). This management approach divides the insulin requirements into two compoand contribute to complications and poor management and contribute to complications and poor health outcomes. bolus therapy (Farkas-Hirsch, 1995; Kaufman, 2012). This and contribute to complications and poor approach the requirements into two approach divides the insulin insulin requirements into background two compocompo- management nents: (i) divides Basal insulin provides the constant health health outcomes. outcomes. approach divides the insulin requirements into two compohealth outcomes. nents: (i) Basal insulin provides the constant background Around 80% of type 1 diabetics rely upon multiple daily nents: (i) Basalby insulin provides the constant constant background supply (i) needed insulin dependent cells. It balances the nents: Basal insulin provides the background Around 80% of type 11 diabetics rely supply needed by insulin dependent cells. It balances the Around 80% of type diabetics rely upon upon multiple multiple daily insulin injections (MDIs) and self-monitoring blooddaily glusupply needed by insulin insulinand dependent cells. It Itofbalances balances the endogenous production consumption glucose the to Around 80% of type 1 diabetics rely upon multiple daily supply needed by dependent cells. insulin injections (MDIs) and self-monitoring blood gluendogenous production and consumption of glucose to insulin injections (MDIs) and self-monitoring blood glucose (SMBG) measurements to administer basal-bolus endogenous production and consumption of glucose to maintain normoglycaemia during fasting. (ii) Bolus insulin injections (MDIs) and to self-monitoring blood gluendogenous production and consumption of glucose to insulin (SMBG) measurements administer basal-bolus maintain normoglycaemia (ii) cose (SMBG) measurements to administer basal-bolus therapy (Pickup, 2011). This to treatment is invasive and maintain normoglycaemia during fasting. (ii) Bolus Bolus insulin insulin balances exogenous glucoseduring inputsfasting. (or disturbances), such cose cose (SMBG) measurements administer basal-bolus maintain normoglycaemia during fasting. (ii) Bolus insulin therapy 2011). and balances exogenous glucose inputs (or disturbances), such therapy (Pickup, 2011). This This treatment is invasive and disruptive(Pickup, to the patient since treatment they need is to invasive receive mulbalances exogenous glucose inputs (or disturbances), such as food ingestion. If food is consumed without extra insulin therapy (Pickup, 2011). This treatment is invasive and balances exogenous glucose inputs (orwithout disturbances), such disruptive to since they to mulas food ingestion. If is extra insulin disruptive to the the patient sinceSMBG they need need to receive receivebeing multiple injections perpatient day, with measurements as foodsupplied ingestion. If food food is consumed consumed without extrahigh insulin being then BGL will rise and remain for disruptive to the patient since they need to receive mulas food ingestion. If food is consumed without extra insulin tiple injections per day, with SMBG measurements being being supplied then BGL will rise and remain high for tiple injections per day, with SMBG measurements being obtained through periodic finger prick testing. Although being supplied then BGL will rise and remain high for an extended period of time. The treatment for this kind injections per day, with SMBG measurements being being supplied then ofBGL will rise and remain highkind for tiple through finger testing. Although an extended The for obtained through periodic finger prick testing. treatment Although invasive, this is stillperiodic one of the mostprick cost effective an extended period period of time. time.provided The treatment treatment for this this kind of disturbance is typically by a bolus (dose) of obtained obtained through periodic finger prick testing. Although an extended period of time. The treatment for this kind of disturbance is provided by aa bolus of invasive, this is is still still one one of of the the most most cost cost effective effective treatment treatment options. this of disturbance is typically typically provided byprior bolus (dose) of invasive, insulin delivered subcutaneously just to, (dose) or with, invasive, this is still one of the most cost effective treatment of disturbance is typically provided by a bolus (dose) of options. insulin delivered subcutaneously just prior to, or with, options. insulin delivered subcutaneously just prior to, or with, meal consumption. insulin delivered subcutaneously just prior to, or with, options. In recent years there has been a large effort aimed at meal meal consumption. consumption. recent there has aa large effort aimed In recent years years there Artificial has been beenPancreas large (AP) effort(Bequette, aimed at at meal consumption. developing a, so called, Current clinical practice determines initial basal and bolus In In recent years there has been a large effort aimed at developing a, so called, Artificial Pancreas (AP) (Bequette, Current clinical practice determines initial basal and bolus developing a, so called, Artificial Pancreas (AP) (Bequette, 2005; Lee et al., 2009; Klonoff et al., 2009; Harvey et al., Current clinical practice determines initial basal and bolus doses based on practice formulae provided initial by clinicians. Many a, so called, Artificial Pancreas (AP) (Bequette, Current clinical determines basal and Many bolus developing et Klonoff al., Harvey et doses based on by 2005; Lee et al., al., 2009; Klonoff et et2014; al., 2009; 2009; Harvey et et al., al., 2010; Lee Cefalu and2009; Tamborlane, Kovatchev doses basedformulae on formulae formulae provided by clinicians. clinicians. Many alternative have provided been developed, along Many with 2005; 2005; Lee et al., 2009; Klonoff et al., 2009; Harvey et al., doses based on formulae provided by clinicians. 2010; Cefalu and Tamborlane, 2014; Kovatchev et al., alternative formulae have been developed, along with 2010; Cefalu and Tamborlane, 2014; Kovatchev et al., 2016). An AP is intended to be a fully operational aualternative formulae have been developed, along with proprietary schemes proposed by individual endocrinoloCefalu and Tamborlane, 2014; Kovatchev et aual., alternative formulae have been developed, along with 2010; 2016). An AP is intended to be a fully operational proprietary schemes proposed by individual endocrinolo2016). An AP is intended to be a fully operational auproprietary schemes schemes proposed proposed by by individual individual endocrinoloendocrinolo- 2016). An AP is intended to be a fully operational auproprietary

Copyright 15609Hosting by Elsevier Ltd. All rights reserved. 2405-8963 © © 2017 2017, IFAC IFAC (International Federation of Automatic Control) Copyright © 2017 15609 Copyright © under 2017 IFAC IFAC 15609Control. Peer review responsibility of International Federation of Automatic Copyright © 2017 IFAC 15609 10.1016/j.ifacol.2017.08.2516

Proceedings of the 20th IFAC World Congress Diego S. Carrasco et al. / IFAC PapersOnLine 50-1 (2017) 15044–15049 Toulouse, France, July 9-14, 2017

tonomous closed loop system for regulating blood glucose by interconnecting a continuous glucose monitor (CGM) and an insulin infusion system (IIS). Such systems could potentially reduce the inconveniences associated with MDI treatment. However, the high cost of such devices, and the limited access to them, reduces the availability of such an approach. Indeed, only about 20% of diabetics use an IIS (Pickup, 2011) (which is needed for an AP). Therefore, there are compelling reasons to study optimal injection algorithms for MDI patients, even more so since the algorithms developed for pumps (including recent closed loop algorithms) are not translatable to MDI therapy since the algorithms developed for pumps rely upon continuous delivery of insulin. The current paper uses optimisation tools to obtain a treatment regimen for MDI patients. To account for unexpected disturbances, the optimisation is implemented in a rolling horizon fashion (Goodwin et al., 2006). The problem of choosing the optimal delivery times and dose for MDIs can be considered as a specific instance of cardinality-constrained optimisation (Markovsky, 2011; Lee and Zou, 2014). In particular, cardinality-constrained optimisation allows one to answer the following question: “Given r allowable injections, when is the optimal time and what is the optimal dose to maintain a suitable blood glucose level (BGL) throughout a day?”. A priori, it is not obvious when the injections should be delivered. To illustrate this point, say that only 3 injections are allowed within the period 6:00am to 10:00pm, in any one of 5 minute intervals. Under these conditions there are over one million possible combinations for the insulin delivery times. Moreover, for each of the insulin delivery time patterns, the dosage would need to be determined. This would require massive computational effort if solved by brute force. In the current paper, we will use a recently developed algorithm for cardinality-constrained optimisation which allows one to solve general problem scenarios. Although initially the optimisation could be performed offline for a 24 hour period, many factors could change during the day, making the initial injection strategy inappropriate. A rolling horizon implementation allows one to update the optimal delivery times and doses as new information arrives in the form of BGL data and food consumed. The treatment strategy described above raises the question, “What is the quantitative trade-off between the number of injections and blood glucose regulation?”. This paper studies the impact of varying the number of allowable MDIs in a given day. More specifically, we address the question, “What is the improvement in performance achieved by allowing one extra insulin injection?”. The results will be exemplified by a model obtained from real data from a patient. The layout for the remainder of the paper is as follows. Section 2 outlines insulin regimens for MDI patients. Section 3 introduces a blood glucose regulation model that accounts for macro-nutrients. Section 4 illustrates the use of cardinality-constrained optimisation when considering a fixed horizon. Section 5 presents algorithms for

15045

implementing cardinality-constrained optimisation with a rolling horizon. Results from simulations with a real patient model are included in Section 6. Conclusions are drawn in Section 7. 2. INSULIN REGIMEN FOR MDI PATIENTS Currently there are two main types of insulin that are available for use: short acting and long acting. Short acting insulin has a quick release and relatively short duration once injected into the body. Long acting insulin has a slow release and long duration once injected into the body. Most treatment regimes use a combination of both types of insulin. Long acting insulin is used to provide the “basal rate” and short acting is used for the “bolus”. The insulin regimen is tailored to the individual. However each regimen is typified by multiple insulin injections and blood glucose monitoring. Type 1 diabetics aim to maintain a constant basal rate, and bolus whenever a meal/food is consumed. The goal of treatment is to regulate a patient’s BGL within a range of 4 − 8 [mmol/l] (70 − 140 [mg/dl]) (Kenny, 2014).

The unit of measurement for doses of insulin is the Unit (U), where 1 [U ] is equivalent to 0.01 [ml]. An Exchange (EX) measures the amount of carbohydrate content within food, where 1 [EX] is equivalent 15 [g]. Note that this does not account for the glycaemic index (GI) content of the food.

Throughout the remainder on the paper the term long acting insulin and basal will be used interchangeably. Similarly, the term short acting insulin and bolus, will be used interchangeably. In the sequel, when referring to MDIs, these only include the bolus component of injections and not the basal injections. 3. LINEAR BLOOD GLUCOSE MODEL WITH MACRO-NUTRIENTS There are many blood glucose models available in the literature, ranging from simple linear impulse response models to non-linear state-space models. For simplicity a linear blood glucose model will be used in the sequel. The results can, in principle, be extended, mutatis mutandis, to any model. The novelty is that the model differentiates the responses to different types of food depending on the GI, i.e. high, moderate, or low. We show that GI content has a significant impact on the BGL response. The model used takes the form: y = Gba (s)uba + Gbo (s)ubo + Gh (s)dh + Gm (s)dm + Gl (s)dl + C (1) where uba , ubo , dh , dm , and dl denote the basal insulin flow rate, bolus insulin flow rate, high GI carbohydrate consumption rate, moderate GI carbohydrate consumption rate, and low GI carbohydrate consumption rate. The above inputs and disturbances are respectively linked to the BGL response by the linear transfer functions Gba , Gbo , Gh , Gm , and Gl . C is a constant which captures the steady-state balance between endogenous glucose production and its consumption. The output y is the resultant BGL.

15610

Proceedings of the 20th IFAC World Congress 15046 Diego S. Carrasco et al. / IFAC PapersOnLine 50-1 (2017) 15044–15049 Toulouse, France, July 9-14, 2017

GI carbohydrates, and low GI carbohydrates consumed within the given time interval (k − 1)∆ ≤ t < k∆, ¯ ba , G ¯ bo , where k ∈ N>0 and ∆ is the sample period. G ¯h, G ¯ m , and G ¯ l are discrete transfer functions linked G to the aforementioned inputs and disturbances. yk is the resultant BGL at time sample k.

118 High GI Moderate GI Low Gi

116 114 BGL [mg/dl]

112 110 108 106 104

4. CARDINALITY CONSTRAINED OPTIMISATION

102 100 0

50

100

150 Time [min]

200

250

300

Fig. 1. Patient data. BGL response to 1 [EX] of High, Moderate and Low GI carbs ingested on a short time interval 100 98

In the case of MDI therapy, the number of injections is limited. Therefore, using the model (8), the problem of selecting the delivery times, and the associated doses, can be formulated as follows: k  P: ui = argmin (yi − y ∗ )2

BGL [mg/dl]

96 94 92 90 Bolus Basal

88

u

86 0

200

400

600

800 Time [min]

1000

1200

Cardinality Constrained Optimisation is a class of optimisation problems where the cardinality of the free variable is constrained. This class of optimisation problems is generally non-convex. Computing the optimal solution is NP-hard. This is because the number of combinations to be tested rises exponentially with the size of the problem.

1400

Fig. 2. Patient data. BGL response to 1 [U ] of Bolus and Basal Insulin injected on a short time interval The parameters for this model were obtained by collecting data from a patient (see Fig. 1 and 2) and using standard system identification tools. Figure 1 shows the blood glucose concentration response to 1 [EX] of high, moderate and low GI carbohydrates applied as a discrete impulse at t = 0 with no insulin injection. Figure 2 shows the blood glucose response to 1 [U ] of bolus and basal insulin applied as discrete impulses. Note that there is approximately a 2 : 1 difference between responses to a high GI and low GI meal.

where ymin

i=1

u = u0 , . . . , uN −1 ui ∈ [0, umax ] yi ≥ ymin card (ui ) ≤ r is a lower bound and y ∗ is the reference value.

This problem is a quadratic program that includes a cardinality constraint. There are several existing tools to solve this problem. For example, by using Mixed Integer Quadratic Programming or cardinality-constraint optimisation. In this paper we adopt the latter approach. In particular we follow the approach described in Dattorro (2005); Delgado et al. (2016) that reformulates a problem with a cardinality constraint as an optimisation problem subject to bilinear constraints. The approach in Delgado et al. (2016) is based on the fact that for x ∈ Rn , then n card(x) n ≤ r is equivalent to ∃w ∈ {w ∈ R |0 ≤ wi ≤ 1, i=1 wi = n − r}, such that xi wi = 0 for i = 1, . . . , n.

The resulting transfer functions obtained from the impulse responses shown in Fig. 1 and 2 are: −5.76 · 10−4 (2) Gba = 2 The main advantage of this approach is that the resulting s + 0.0084s + 1.764 · 10−5 −2 problem can be solved with standard tools in non-linear −3.6 · 10 (3) programming. Based on this idea, the solution of problem Gbo = 2 s + 0.0336s + 0.0002822 P is equivalent to the solution of the following optimisa4.277 · 10−5 tion problem: Gh = 5 s + 0.27s4 + 0.029s3 + 0.0016s2 + 4.25 · 10−5 s + 4.59 · 10−7 k  (4) P : min (yi − y ∗ )2 equiv u,w 4.752 · 10−3 i=1 (5) Gm = 3 s + 0.12s2 + 0.004752s + 6.221 · 10−5 ui ∈ [0, umax ] 3.24 · 10−3 yi ≥ ymin (6) Gl = 3 2 s + 0.12s + 0.0045s + 5.4 · 10−5 k  C = 108 [mg/dl] (7) wi = k − r

Although the identified system is continuous, it is necessary to work in discrete time due to the optimisation algorithm, as shown in Section 4. Therefore we introduce the discretised version of the model as follows: ¯ ba (z)ubak + G ¯ bok (z)ubok + G ¯ h (z)dhk yk = G ¯ ¯ + Gm (z)dmk + Gl (z)dlk + C (8) where ubak , ubok , dhk , dmk , and dlk denote the basal insulin, bolus insulin, high GI carbohydrates, moderate

i=1

for i = 1, . . . , k 0 ≤ wi ≤ 1 ui wi = 0 for i = 1, . . . , k This problem can be solved using known techniques in non-linear programming (Bertsekas, 1999).

−1 In summary, given a food profile {dhk , dmk , dlk }N k=0 , a N −1 basal insulin profile {ubak }k=0 , and restricting the number of allowable MDIs over a specific time period to r, the

15611

Proceedings of the 20th IFAC World Congress Diego S. Carrasco et al. / IFAC PapersOnLine 50-1 (2017) 15044–15049 Toulouse, France, July 9-14, 2017

optimal delivery times and dose for the bolus insulin are obtained by solving Pequiv .

BGL (mmol/l)

Blood Glucose Response with Two MDIs

5. ROLLING HORIZON IMPLEMENTATION

12 10 8 6 4 6

Exchanges (EX)

5.1 Event Triggered Rolling Horizon Implementing an optimisation-based control algorithm in a rolling horizon fashion allows the system to account for any update of relevant information. In the case of MDI therapy, the main source of state information are SMBG measurements. However, rolling horizon optimisation also allows the system to respond to unexpected changes in the daily routine, such as missing a meal, missing an injection or performing unanticipated exercise. A key point is that a re-calculation of the optimal insulin delivery times is carried out only when triggered by an external event. In addition, a special strategy is needed to ensure the cardinality-constraint does not become redundant. External events provide an update of information to the system at infrequent times. The new information, which may include actual BGL, food, and insulin data, is used to obtain a new state estimate of the system. Such an implementation is non-standard since the information update does not occur at regular times. This can be overcome with a non-equal sampling Kalman filter. Another matter of importance is the window of the horizon. The only control action in this system is insulin, which acts to lower blood glucose. Insulin cannot be removed nor can a control action be used to increase the blood glucose. The use of glucagon has been considered as a control action to raise BGL, however it is not currently used in practice due to complications with its storage at room temperature. For this reason it is vital that the rolling horizon window looks far enough into the future to ensure any delivered insulin does not cause a hypoglycaemic event at a later time. Here a 24 hour window is used. 5.2 Enforcing Cardinality-Constraint Care must be taken to ensure that the cardinalityconstraint is enforced, i.e. taking into account injections that have already been delivered. Thus, the following additional constraints are added to Pequiv : for t =[ko , ki ] (9) ud ≤ r for t =[ki+1 , kN ] (10) ua = r − u d up = r for t =[kN +1 , kN +i ] (11) where ki is the current sample, N is the number of samples within a day, ud is the number of MDIs delivered so far within the current day, ua is the number of MDIs remaining within the current day, and up is the number of MDIs allowed for the next day. 6. SIMULATIONS In this section we illustrate how the Cardinality Constrained Optimisation algorithm described above can be applied to MDI therapy. We study the incremental performance gain due to allowing an extra insulin injection. We then illustrate the impact of rolling horizon implementation with several daily scenarios.

15047

8

10

12

14

16

18

20

22

24

18

20

22

24

14 16 18 Time of Day (hours)

20

22

24

Food Profile 5 0

6

8

10

12

14

16

Units (U)

Insulin Profile 30

Basal

20 10 0

6

8

10

12

Fig. 3. Blood Glucose Response and Optimisation Results for Two MDIs r = 2 (solid). No injections (dashed) 6.1 Fixed Horizon Using the algorithm described in Section 4 and the model described in Section 3, we obtain the resulting insulin injection schedule for a specific daily food profile, i.e. a combination of low, moderate and high GI food throughout the day. Figure 3 shows the result for r = 2 (2 MDIs) The solid line in subplot 1 represents the blood glucose response of the model when the given MDIs are allowed and the dashed line represents the blood glucose response of the model without any injections (r = 0). The latter is included for comparison purposes. The food profile is shown in subplot 2. The resulting insulin injection schedule is shown in subplot 3. 6.2 Incremental Performance Gain The next question studied is what is the quantitative gain achieved by increasing the number of MDIs? To study this question, offline optimisations were performed for r = 0, . . . , 10. The performance of each optimisation was determined by calculating the blood glucose response deviation from the ideal (but unattainable) result of a constant blood glucose response (y = 6.0) throughout the whole day. These were then normalised against the blood glucose response for no MDIs. Hence an error of 1 corresponds to the worst blood glucose response (when disturbances are not dealt with) and an error of 0 corresponds to the best blood glucose response (when disturbances are completely mitigated). The results are shown in Fig. 4. Figure 4 shows that each additional allowable MDI improves the blood glucose response. With 4 MDIs a 90% improvement has been achieved with respect to blood glucose regulation. Interestingly, no significant benefit is obtained by increasing the MDIs beyond 4 (for this particular food and basal insulin profile). Of course this result could change if different food and basal insulin profiles were adopted. However, it raises the question whether there is an acceptable performance trade-off between the number of insulin injections and continuous insulin supply.

15612

Proceedings of the 20th IFAC World Congress 15048 Diego S. Carrasco et al. / IFAC PapersOnLine 50-1 (2017) 15044–15049 Toulouse, France, July 9-14, 2017

Improvement in Blood Glucose Regulation 1

Food Profile 100

CHO (g)

0.9

0.7

50

0

0.6

6

8

10

12

0.5

16

18

20

22

24

6

0.4 0.3 0.2 0.1

Original Updated

4 2 0 6

0

2

4 6 Allowable MDIs

8

8

10

12

14

16

18

20

22

24

Modelled BGLs

10

BGL (mmol/l)

0

14

Insulin Profile Units (U)

Normalised Deviation Error

0.8

Fig. 4. Effect of Increasing MDIs on the BGL

10 Original Updated

8 6 4

6.3 Rolling Horizon

6

8

10

12

14

16

18

20

22

24

18

20

22

24

18

20

22

24

18

20

22

24

Time (hours)

Fig. 5. Traces for Scenario 1 Food Profile

CHO (g)

In this section we illustrate the benefits of a rolling horizon implementation. Consider a typical day for a patient having breakfast at 08:00 (90g CHO), lunch at 13:00 (75g CHO), an afternoon snack at 16:00 (15g CHO) and dinner at 19:00 (75g CHO). All moderate GI meals. We will consider three scenarios. In Figures 5, 6 and 7 the top, mid and bottom panels show the food profile, the insulin profile and the simulated BGL trace.

100 50

Scheduled Unscheduled

0 6

Scenario 1 The patient has a normal lunch at 13:00 but forgets to inject the corresponding insulin bolus. At 14:00 the algorithm is updated. The results are shown in Figure 5.

Scenario 2 The patient has an unscheduled snack of 120 [g] at 18:30, right before dinner. The algorithm is updated immediately. The results are shown in Figure 6. The top panel in Figure 6 shows the additional unscheduled meal. The mid panel shows the updated insulin profile, with the original injection schedule shown with a dashed line. The bottom panel shows the BGL traces had the original injection schedule occurred (dashed) despite the new meal, i.e. assuming the algorithm was not updated regarding the unscheduled meal, and the BGL trace corresponding to the updated injection schedule (solid). Note that the rolling horizon optimisation has significantly improved the BGL response. In particular, it can be seen that informing the algorithm immediately about the unscheduled snack shifts the third injection to the current time and increases the dosing by almost double. This

10

12

14

16

Insulin Profile

Units (U)

10

The mid panel in Figure 5 shows the original injection schedule with a dashed line. The bottom panel shows the BGL trace for the updated injection schedule (solid line) and the BGL trace if the original injection schedule had been used as planned.

Original Updated

5

0 6

8

10

12

14

16

Modelled BGLs

BGL (mmol/l)

It can be seen that, given the update at 14:00, the algorithm recalculates the new optimal time and injection dose, which is to bolus 1.9 units of insulin at 14:00. Most importantly, note that this bolus is not the same size as the originally planned bolus. Had the same bolus been administered, the patient would have had a hypoglycaemic event. Note that the delayed injection of insulin has led to an unavoidable decrease in performance.

8

12 10 8 6 4

Original Updated

6

8

10

12

14

16

Time (hours)

Fig. 6. Traces for Scenario 2 avoids a considerable hyperglycaemic event had the event not been informed. Scenario 3 The patient has an unscheduled snack of 120 [g] at 18:30, right before dinner. The original dinner and the third insulin injection occur as planned. However, importantly the algorithm is not updated until 19:30 with the new information about the unscheduled snack. The results are shown in Figure 7. Since the three allowed daily injections have already occurred, the two options for the algorithm are to do nothing or to suggest a fourth injection. The mid panel in Figure 7 shows the timing and dosage of the additional injection, if allowed. The bottom panel shows the BGL traces corresponding to the original three injection schedule (dashed) and the four injection schedule (solid). It can be seen that, although the unscheduled snack is much bigger in comparison to the other meals, the extra injection is smaller than the other ones. The algorithm

15613

Proceedings of the 20th IFAC World Congress Toulouse, France, July 9-14, 2017 Diego S. Carrasco et al. / IFAC PapersOnLine 50-1 (2017) 15044–15049

CHO (g)

Food Profile 100 50

Scheduled Unscheduled

0 6

8

10

12

14

16

18

20

22

24

18

20

22

24

18

20

22

24

Insulin Profile

Units (U)

6 4 Original Additional

2 0 6

8

10

12

14

16

BGL (mmol/l)

Modelled BGLs 12 10 8 6 4

Original Additional

6

8

10

12

14

16

Time (hours)

Fig. 7. Traces for Scenario 3 takes into account the difference in time and the fact that the meal was compensated when the algorithm was updated. Had a normal dose for that snack been injected at 19:30 the patient would have had a severe hypoglycaemic event. Indeed, this exact scenario was taken from a real life experience. The patient, a 9 year old child in this case, did not immediately inform his parents of two chocolate bars he ate. When he told his parents several hours later the parents injected the same dose they would have used had the snacks been ingested at that time. He ended up having an hypoglycaemic event and taken to the hospital. These scenarios provide typical complex scenarios that a patient may have to address in the course of their daily activity. Current treatment strategies do not provide an objective patient response to these scenarios as there are no stipulated rules for the actions required in such situations. We believe that they would provide a significant challenge to even the most attentive patient or health care professional. 7. CONCLUSION In this paper we have shown how Cardinality Constrained Optimisation can be applied to the problem of Multiple Daily insulin Injections in Type 1 Diabetes therapy. We have also introduced a rolling horizon scheme that allows for updates as the patient’s day evolves. Simulations have shown the advantages of the algorithm as well as the versatility to deal with unexpected situations. We have shown that the proposed algorithm is capable of diminishing the complexity of decision making by providing systematic, quantitative advice. Moreover, the receding horizon implementation addresses changing daily situations. REFERENCES Aronoff, S.L., Berkowitz, K., Shreiner, B., and Want, L. (2004). Glucose Metabolism and Regulation: Beyond Insulin and Glucagon. Diabetes Spectrum, 17(3), 183– 190. Bequette, B. (2005). A critical assessment of algorithms and challenges in the development of a closed-loop

15049

artificial pancreas. Diabetes technology & therapeutics, 7(1), 28–47. Bertsekas, D.P. (1999). Nonlinear programming. Athena scientific Belmont. Cefalu, W.T. and Tamborlane, W.V. (2014). The artificial pancreas: are we there yet? Diabetes Care, 37(5), 1182– 1183. Dattorro, J. (2005). Convex optimization & Euclidean distance geometry. Meboo Publishing, USA. Davidson, P., Hebblewhite, H., Steed, R., and Bode, B. (2008). Analysis of guidelines for basal-bolus insulin dosing: basal insulin, correction factor, and carbohydrateto-insulin ratio. Endocrine Practice, 14(9), 1095–1101. Delgado, R.A., Aguero, J.C., and Goodwin, G.C. (2016). A novel representation of rank constraints for real matrices. Linear Algebra and its Applications, 496, 452 – 462. Farkas-Hirsch, R. (1995). Intensive diabetes management. American Diabetes Association. Goodwin, G.C., Seron, M.M., and De Don´a, J.A. (2006). Constrained control and estimation: an optimisation approach. Springer Science & Business Media. Harvey, R.A., Wang, Y., Grosman, B., and et al. (2010). Quest for the artificial pancreas: Combining technology with treatment. IEEE Engineering in Medicine and Biology Magazine, 29(2), 53–62. International Diabetes Federation (2015). IDF Diabetes Atlas, 7th edn. Brussels, Belgium. Kaufman, F. (2012). Medical Management of Type 1 Diabetes, 6th Ed. American Diabetes Association. Kenny, C. (2014). When hypoglycemia is not obvious: Diagnosing and treating under-recognized and undisclosed hypoglycemia. Primary care diabetes, 8(1), 3–11. King, A.B. and Armstrong, D.U. (2007). A prospective evaluation of insulin dosing recommendations in patients with type 1 diabetes at near normal glucose control: basal dosing. Journal of diabetes science and technology, 1(1), 36–41. Klonoff, D.C., Cobelli, C., Kovatchev, B., and Zisser, H.C. (2009). Progress in development of an artificial pancreas. Journal of diabetes science and technology, 3(5), 1002–4. Kovatchev, B., Tamborlane, W.V., Cefalu, W.T., and Cobelli, C. (2016). The artificial pancreas in 2016: a digital treatment ecosystem for diabetes. Diabetes Care, 39(7), 1123–1126. Lee, H., Buckingham, B.A., Wilson, D.M., and Bequette, B.W. (2009). A closed-loop artificial pancreas using model predictive control and a sliding meal size estimator. Journal of diabetes science and technology, 3(5), 1082–90. Lee, J. and Zou, B. (2014). Optimal rank-sparsity decomposition. Journal of Global Optimization, 60(2), 307– 315. Markovsky, I. (2011). Low rank approximation: algorithms, implementation, applications. Springer Science & Business Media. Pickup, J. (2011). Insulin pumps. International Journal of Clinical Practice, 65, 16–19. Walsh, J., Roberts, R., and Bailey, T. (2011). Guidelines for optimal bolus calculator settings in adults. Journal of diabetes science and technology, 5(1), 129–135.

15614