Optimal regulation of blood glucose level in Type I diabetes using insulin and glucagon.

The Glucose-Insulin-Glucagon nonlinear model accurately describes how the body responds to exogenously supplied insulin and glucagon in patients affected by Type I diabetes. Based on this model, we design infusion rates of either insulin (monotherapy) or insulin and glucagon (dual therapy) that can optimally maintain the blood glucose level within desired limits after consumption of a meal and prevent the onset of both hypoglycemia and hyperglycemia. This problem is formulated as a nonlinear optimal control problem, which we solve using the numerical optimal control package [Formula: see text]. Interestingly, in the case of monotherapy, we find the optimal solution is close to the standard method of insulin based glucose regulation, which is to assume a variable amount of insulin half an hour before each meal. We also find that the optimal dual therapy (that uses both insulin and glucagon) is better able to regulate glucose as compared to using insulin alone. We also propose an ad-hoc rule for both the dosage and the time of delivery of insulin and glucagon.

1 Introduction

Insulin and glucagon are pancreatic hormones that help regulate the levels of glucose in the blood [1-4]. Insulin is produced by the beta-cells in the pancreas and carries glucose from the bloodstream to the cells throughout the body. Glucagon releases glucose from the liver into the bloodstream in order to prevent hypoglycemia. In people affected by diabetes insulin is either absent (type I diabetes) or not produced in the proper amount (type II diabetes). In type I diabetes the body’s immune system attacks and destroys the beta cells. As a result, insulin is not produced and glucose accumulates in the blood which may cause serious harm to several organs. Type II diabetes is a metabolic disorder in which the beta cells are unable to properly regulate the blood glucose within limits. Common therapies for diabetes involve the administration of exogenous insulin. Currently glucagon is not typically included in therapies because it does not preserve its chemical properties at room temperature and also because diabetic patients are still able to produce it.

The control of glucose levels in diabetic patients is an active field of research [5-19]. The approval by the FDA of a simulator which replaces in-vivo with in-silico therapy testing has greatly benefited this area of research. This simulator implements a mathematical model, first proposed in [1] and updated in [2-4], and provides an alternative to often slow, dangerous and expensive human testing.

Typically, insulin is administered manually approximately half an hour before each meal where the amount is determined from the current glucose level (measured through a blood sugar test), the expected glucose intake, and the patient’s sensitivity to insulin. In what follows we will refer to this as the standard therapy. In 1992 the first insulin pumps were introduced to the market. They delivered both a consistent basal amount of insulin and an insulin bolus determined by the patients based on their glucose level. It was only in 2016 that the first autonomous system for glycemic control was approved by the FDA. The system consists of an insulin pump, a sensor that measures the blood glucose level continuously in time, and control software that is able to regulate the insulin level in the blood without needing any input from the patient.

Many control techniques have been proposed and tested to regulate blood glucose levels using insulin pumps including PID (proportional–integral–derivative) control [5, 6, 8–10, 20], fuzzy logic control [11-13] and bio-inspired techniques [14] which do not rely on a mathematical model. In [21] closed loop control has been used on a so called “minimal model” [22-24]. In [15–18, 25] a linear model predictive control (MPC) has been used in a model with fixed structure but for which parameters are constantly updated to adapt to the patient’s response. In [26] linear MPC has been used in silico. In [27] MPC has been applied to a system linearized around the operating points of a physically derived nonlinear model and in [28] multiple model probabilistic predictive control has been used. In [29] MPC has been applied together with a moving horizon estimation technique to a linear model. Most of the models used when designing the above controllers are simplified versions of the FDA approved model and all the control techniques considered only use insulin (but not glucagon) as control input.

Because insulin delivered exogenously is not subject to normal physiological feedback regulation, hypoglycemia is common in patients with Type 1 diabetes who undergo treatment [30]. For these patients it has been proposed that exogenous insulin can be used to lower their blood glucose level and exogenous glucagon can be used to prevent hypoglycemia [31, 32]. Currently, a commercial pump that delivers both insulin and glucagon is not available, and the development of a two-drug artificial pancreas is still the subject of clinical research [33-41]. An outstanding research question, which we address in this paper, is the determination of the temporal dosages of both insulin and glucagon, in the case of the dual therapy.

Following the study in [42] which optimized multi-drug therapies for autophagy regulation, here we seek to determine an optimal strategy for delivery of both insulin and glucagon. We consider the combined effects of insulin and glucagon in regulating blood glucose levels in patients with Type 1 diabetes, using the model in [3] and nonlinear optimal control theory. Additionally, the objective function that we seek to minimize is the Blood Glucose Index which is a well known tool to measure the risk for a patient to enter either hyperglycemia or hypoglycemia. To design the optimal control problem, we use the balance control technique of ref. [43], which introduces a trade-off between the error allowed with respect to a state based cost (Blood Glucose Index) and the control effort. Our goal is to evaluate the performance limits of a control algorithm in the blood glucose problem, and to discuss the advantages of the dual drug therapy compared to the single drug therapy. Note that even though we do not attempt to design a closed-loop control strategy that works without the patient’s intervention, the solution we propose can be adapted for that purpose.

From solving the optimal control problem for a family of objective functions derived from the balance control paradigm, we observe the emergence of a pattern, from which we propose a simple rule for the delivery of insulin and glucagon similar to the standard therapy, but for the case that both insulin and glucagon are used. While this therapy is suboptimal, we see that it still performs better than the optimal solution with insulin alone.

Finally, we test the robustness of the optimal solution. While optimal control does not guarantee robustness of the optimal solution with respect to model uncertainty or parameter mismatches, we see that our proposed solution still performs well in the presence of model parameter perturbations and variations affecting the time and glucose intake of the meal.

2 Materials and methods

2.1 Model and parameters

We consider the model in [3, 4] which is a system of nonlinear ordinary differential equations (ODEs). The equations are given in Eqs (S1)-(S9) in S1 Appendix. We write the ODEs in Eqs (S1)-(S9) in the form where the state vector is x = [x1(t), x2(t), …, x17(t)] and t is the physical time (in min). In Table 1 we tabulate all of the variables x and their names. The control input vector is u(t) = [u(t), u(t)], where u(t) ≥ 0 is the exogenous insulin infusion rate (in insulin Unit/min) and u(t) ≥ 0 is the exogenous glucagon infusion rate (mg/min). Both u(t) and u(t) are the external inputs to the system in Eq (1). A schematic network diagram has been presented in S1 Fig The scalar quantity D(t) represents the exogenous glucose input, that is, the glucose intake with a meal. The output of the system is the quantity G(t), which measures the density of glucose in the blood, obtained as the ratio between the plasma glucose and the distribution volume of glucose V.

Table 1

Variables and their physical meaning.

Variables	Names	Representing	Units
x1	Gp	Mass of glucose in plasma	mg/kg
x2	Gt	Mass of glucose in tissue	mg/kg
x3	Il	Mass of insulin in liver	pmol/kg
x4	Ip	Mass of insulin in plasma	pmol/kg
x5	I′	Mass of delayed in compartment 1	pmol/L
x6	XL	Amount of delayed insulin action on EGP (Endogenous glucose production)	pmol/L
x7	Qsto1	Amount of solid glucose in stomach	mg
x8	Qsto2	Amount of liquid glucose in stomach	mg
x9	Qgut	Amount of glucose in intestine	mg
x10	X	Amount of interstitial fluid	pmol/L
x11	SRHs	Amount of static glucagon	ng/L/min
x12	H	Amount plasma glucagon	ng/L
x13	XH	Amount of delayed glucagon action on EGP	ng/L
x14	Isc1	Amount of nonmonomeric insulin in the subcutaneous space	pmol/kg
x15	Isc2	Amount of monomeric insulin	pmol/kg
x16	Hsc1	Amount of glucagon in the subcutaneous space 1	ng/L
x17	Hsc2	Amount of glucagon in the subcutaneous space 2	ng/L

State variables and their physical meaning.

When u(t) = 0, u(t) = 0 and D(t) = 0, the model reaches (for physically meaningful parameters) a steady state, also known as the basal condition of a patient. The basal condition depends upon the parameters of the models Θ. We denote by a set of parameters for which the basal glucose level G is equal to G. The basal levels for the other states are found according to Eq (S10).

2.2 Problem formulation

We formulate a nonlinear optimal control problem with two control goals. The first goal is to regulate the glucose at levels corresponding to low clinical risk of either hyperglycemia or hypoglycemia during a time period over which a meal is consumed. We assume that a meal is ingested at time t = τ, which we assume to be modeled as a Dirac delta function D(t) = Dδ(t − τ). To evaluate the clinical risk of a particular glycemic value, Kovatchev et al. [44, 45] proposed the Blood Glucose Index (BGI), defined as where a small BGI value corresponds to low risk of either hyperglycemia or hypoglycemia. This metric also takes into account the fact that (i) the target blood glucose range as defined by the Diabetes Control and Complications Trial [46] (between 70 and 180 mg/dL) is not symmetric about the center of the range and (ii) hypoglycemia occurs at glucose levels closer to the basal level than hyperglycemia. The second goal is to limit the overall usage of insulin and/or glucagon over the period [t0, t].

We formulate the optimization problem according to these two goals, subject to the following constraints,

In Eqs (2) and (3), the insulin infusion rate u(t) and the glucagon infusion rate u(t) are the two control inputs. The three coefficients α, α and α in Eq (2) are tunable factors through which we may vary the weight associated with each of the three terms in the cost function J. The first coefficient, α is dimensionless while the units of α and α are (U/min)− and (mg/min)−, respectively. Note that by setting u = 0 in Eq (6), we have an optimal control problem in terms of insulin only.

The first term in the objective function (2) defines a regulation problem, i.e., we try to maintain the glucose at low risk levels. The second and third terms in the cost function are chosen to avoid using excess insulin or glucagon. For p = 1 in Eq (2), the second and third terms define a ‘minimum fuel’ problem, thus we call the optimization problem ReMF (Regulation and Minimum Fuel). In this case, we expect the optimal solution to consist of pulsatile inputs and [47, 48]. For p = 2, the second and third term inside the cost function define a ‘minimum energy’ problem, thus we call the optimization problem ReME (Regulation and Minimum Energy). In this case, we expect the optimal control inputs and to be continuous. The set of equations in (3) coincide with the ODEs in Eqs (S1)-(S9) of S1 Appendix. In Eq (4) G and G are the lower and upper bounds for G(t), they can be set in order to avoid undesired hypoglycemic or hyperglycemic states. In Eqs (5) and (6) and are upper bounds for the insulin and glucagon delivery rates, respectively. These constraints are set by the maximum infusion rates allowed by the insulin pump. In Eq (5) is the lower bound for u(t), i.e., a minimum insulin delivery rate that can be used to set a basal insulin infusion rate to counteract endogenous glucose production [49]. Finally, in Eqs (7) and (8)), and set limits to the total limits of insulin and glucagon that can be delivered over the time period [t0, t]. The initial condition in Eq (9) defines the patient’s condition before administration of the therapy. In the Results section, we discuss how we choose the bounds on G(t), u(t), u(t), ϕ, ϕ, the control time period [t0, t] and the initial condition .

Our goal is to find an optimal solution which satisfies the constraints in Eqs (3)–(9) and minimizes the objective function in Eq (2). Note that the BGI only depends upon G(t): we are making no attempt to control the states of the system, only its output. In the literature, such an approach is often referred to as target control [47, 50].

2.3 Method: Pseudo-spectral optimal control

Optimal control theory combines aspects of dynamical systems, optimization, and the calculus of variations [47] to solve the problem of finding a control law for a given dynamical system such that the prescribed optimality criteria are achieved. The Eqs (2) and (3)–(9) together form a constrained optimal control problem, which can generally be written as, In general, there exists no analytic framework that is able to provide the optimal time traces of the controls u*(t) and the states x*(t) in (10), and so we must resort to numerical techniques.

Pseudo-Spectral Optimal Control (PSOC) is a computational method for solving optimal control problems. Here we present a brief overview of the theory of pseudo-spectral optimal control. PSOC has become a popular tool in recent years [51, 52] that has let scientists and engineers solve optimal control problems like Eq (10) reliably and efficiently in applications such as guiding autonomous vehicles and maneuvering the international space station [52]. PSOC is an approach by which an OCP can be discretized by approximating the integrals by quadratures and the time-varying states and control inputs with interpolating polynomials. Here we summarize the main concept of the PSOC. We choose a set of N discrete times {τ} i = 0, 1, …, N where τ0 = −1 and τ = 1 with a mapping between t ∈ [t0, t] and τ ∈ [−1, 1]. The times {τ} are chosen as the roots of an (N + 1)th order orthogonal polynomial such as Legendre polynomials or Chebyshev polynomials. The choice of dicretization scheme is important to the convergence of the full discretized problem. For instance, if we choose the roots of a Legendre polynomial as the discretization scheme, the associated quadrature weights can be found in the typical way for Gauss quadrature. The time-varying states and control inputs are found by approximating them with Lagrange interpolating polynomials, where and are the approximations of x(τ) and u(τ), respectively, and L(τ) is the ith Lagrange interpolating polynomial. The dynamical system is approximated by differentiating the approximation with respect to time. Let which allows one to rewrite the original dynamical system constraints in (10) as the following set of algebraic constraints. The integral in the cost function is approximated as, The original time-varying states, control inputs, the dynamical equations constrained and the cost function are now discretized approximation of the continuous NLP problem. Thus the discretized approximation of the original OCP is compiled into the following nonlinear programming (NLP) problem. We have used [53], an open-source PSOC library, to perform the above PSOC discretization procedure. The NLP in (16) can be solved with a number of different techniques, but here we use an interior point algorithm [54] as implemented in the open-source software Ipopt [55].

2.3.1 Continuous approximation of non-differential function in ODEs

The optimization algorithms implemented in require the derivatives of the function exists. As there are terms that contain discontinuities in Eqs (S1)-(S9), we replace them with smooth approximations which are described in S3 Appendix.

3 Results

We now describe in more detail the optimal control problem in Eqs (2)–(9) by setting the constraint and parameter values. In Fig 1 we plot the function BGI(G) versus the glucose G. The minimum BGI(G) occurs at G = G = 112.51 mg/dL, which corresponds to a clinical target set for the glucose level [46]. Based on the data in [56], the average fasting plasma glucose level of patients with type I diabetes is G = 130 (mg/dL). Thus, we set the the basal glucose level G = 130 (mg/dL). The parameters are set so that the steady state glucose is 130 (mg/dL) in the absence of a meal and of exogenously supplied insulin, i.e., we compute Θ130.

Fig 1

A) The Blood Glucose Index (BGI(G(t))) as a function of the blood glucose G(t). The function is minimized at G(t) = G = 112.51 (mg/dL). B) The response of glucose (G(t)) to different time-constant basal insulin infusion rates in the absence of a meal. We see that as u increases, the glucose is further down regulated.

We set the upper and lower bounds for the glucose level, G and G in Eq (4), to satisfy the target blood glucose range, 90 ≤ G(t) ≤ 180 [46]. The control time period is [t0, t] = [0, 300] minutes, and we assume that a meal with 70 grams of glucose is consumed at time t = 60 min (i.e. D(t) = 70δ(t − 60)).

We consider a situation in which the patient’s glucose level is partially controlled by providing a constant but low insulin infusion rate u > 0 (which is common for patients who use an insulin pump) [49] and serves to compensate for the endogenous glucose production. In Fig 1 we show glucose response G(t) for different values of constant u in the absence of a meal. We observe that for u = 0.0024 (U/min), G(t) converges to the desired glucose level G. We thus set the lower bound of u(t) in Eq (5), , while its upper bound is set to (in U/min), the maximum insulin flow allowed in commercial pumps [57]. In the absence of commercially available glucagon pumps, we will assume that a pump mechanically similar to an insulin pump is used to deliver glucagon. Since the maximum flow rate for an insulin pump is 0.15 mL/min (1 mL of insulin solution contains 100 U of insulin), and normally 1 mg of glucagon is diluted in 1 mL of solution, the maximum glucagon flow rate in Eq (6) is set to mg/min.

The amount of insulin administered in a bolus to a patient with a basal glucose level lower than 150 mg/dL normally ranges between 0.12 and 0.2 U/kg [58]. As the body mass of the in-silico patient we consider is 78 kg, we set U in Eq (7). The maximum total amount of glucagon administered in one shot to a patient who is in a hypoglycemic state is 1 mg, and a second identical shot can be administered after thirty minutes. We thus choose the maximum total amount of glucagon used (as defined in Eq (8)) throughout the five hour therapy to be mg.

The choice of the initial condition in Eq (9) is critical. We select the initial condition so that the solution of our optimal control problem only attempts to regulate glucose in response to a meal. In the results presented we have set the initial condition equal to the values of the states when u(t) = u after a period of fasting (the final point of the blue curve in Fig 1). If we were to select any alternative initial condition then the solution to the optimal control problem would try to ‘correct’ the initial condition as well, making comparisons between solutions difficult.

Once the parameters, bounds, the control time period and the initial condition are set, we solve the nonlinear optimal control problem using . We first solve the optimal control problem without glucagon (i.e., ), and then we solve the optimal control problem using both insulin and glucagon.

To evaluate the effectiveness of the obtained results, we introduce the following measures.

The cumulative insulin r(t) and cumulative glucagon r(t) used up to time t,

The total amount of insulin ϕ = r(t) and the total amount of glucagon ϕ = r(t) used up to final time t.

The integral of BGI over the entire time period [t0, t], where a large Δ indicates that the patient is at higher risk of either hyperglycemia or hypoglycemia for a prolonged period of time.

The maximum and minimum values attained by the blood glucose level over the entire time period [t0, t], which measure the risk for either hyperglicemia or an hypoglicemia [45, 59], respectively.

3.1 Insulin as control input

In this section we use only insulin as control input, i.e., we set u = 0 in Eq (3). As the orders of magnitude of the terms BGI and in the objective function are different, it is important to find the appropriate values of the scaling factors α and α. In what follows, we use a Pareto-front analysis to determine these values. We first rewrite the objective function as where ε = α/α. In Fig 2( we plot Δ, Gmin, Gmax and ϕ as functions of the coefficient ε. By looking at these plots, we see that the four measures can be divided into two groups. On the one hand, Δ and Gmax (panels A and C), improve (decrease) as ε increases, with a sharp transition around ε = 10 for the ReMF problem and around ε = 103 for the ReME problem. On the other hand, Gmin and ϕ (panels B and D), behave in the opposite way, i.e., they improve (insulin decreases and the minimum glucose level increases) as ε decreases, again with a sharp transition around ε = 10 for the ReMF problem and around ε = 103 for the ReME problem. Because the four curves in Fig 2( are monotone, all the points are Pareto-efficient, i.e., it is not possible to improve one objective (e.g. Δ) without worsening the other one (e.g. ϕ). We notice that past a certain value of ε (10 in the ReMF case, 103 in the ReME case) Δ and Gmax do not further decrease and Gmin and ϕ remain unchanged. We choose as weights α = 10 and α = 1 for p = 1, while we choose α = 103 and α = 1 for p = 2 (these are highlighted by dashed circles in Fig 2). The reason for these choices (for both values of p) is that these values yield ϕ ∼ 10 units, which is equal to two thirds of the maximum amount of insulin that can be supplied (), and Gmin ∼ 93mg/dL, which is far from the hypoglycemic risk region.

Fig 2

Performance of the optimal control solution as a function of ε.

Large (small) values of ε correspond to a large (small) weight associated with the BGI index in the objective function, compared to the weight for insulin expenditure. The first four plots show our metrics as functions of the objective function coefficients: A) Δ vs. ε, B) Gmin vs. ε, C) Gmax vs. ε, and D) ϕ vs. ε. E) We also project the Pareto front into the Δ - ϕ plane. We see a clear trade-off between Δ and ϕ as we vary ε. By increasing ε we can decrease the values of Δ and Gmax. However, the values of Δ and Gmax do not further decrease for ε larger than 10 for the ReMF problem (p = 1) and the value of Δ does not further decrease for ε larger than 103 for the ReME problem (p = 2). We choose ε = 10 for p = 1 and ε = 103 for p = 2, which are indicated by dashed circles in the figure, for the remaining simulations.

In Fig 2, we plot a projection of the Pareto front in the Δ and ϕ plane. Looking at this plot, the trade-off between Δ and ϕ is evident; if the total amount of insulin expenditure increases, Δ decreases and vice-versa. The ReMF and the ReME therapies can also be compared in Fig 2. The ReMF Pareto front dominates the ReME one (both Δ and ϕ are lower on the blue curve (p = 1) compared to the magenta curve (p = 2)). This indicates that a shot of insulin (the optimal solution of a ReMF problem is typically a pulsatile function) performs slightly better in terms of Δ than a therapy in which the drug is delivered over a longer period of time while using less insulin.

Fig 3 shows the results of the optimal control problem for the selected values of α and α. The blue and magenta curves are the optimal solutions of the ReMF and of the ReME problem, respectively. The orange curve corresponds to the case that 10 U of insulin are injected 30 minutes before the time of the meal, i.e., the standard therapy.

Fig 3

A) The time evolution of glucose G(t) (in mg/dL). The blue curve corresponds to the pulsatile optimal insulin supply rate u(t) (shown in B) obtained by solving the ReMF problem. The magenta curve corresponds to the continuous optimal insulin supply rate u(t) (shown in B) obtained by solving the ReME problem. The orange curve is the time evolution of G(t) corresponding to the standard therapy (10 U of insulin injected 30 minutes before the time of the meal). B) Time evolution of the optimal insulin infusion rates u(t) (in U/min). Color code is consistent with A. C) Cumulative insulin supply r(t) (in U) as a function of t.

We observe that for p = 1 the optimal insulin infusion rate is pulsatile with a pulse appearing at t ∼ 20 minutes, which is 40 minutes before the time of the meal. We obtained qualitatively similar results for different choices of the model parameters, with the pulse typically appearing at a time in the interval t ∈ [20, 30] minutes. It is noteworthy that the optimal solution is close to the standard insulin based therapy for glucose regulation in diabetics. The optimal insulin infusion rate is continuous when we solve the ReME problem, also shown in the inset of Fig 3. Note that the ReMF and ReME therapies perform very similarly with respect to glucose as the peak insulin infusion rate occurs at approximately the same time and the total amount of insulin administered is nearly equal.

3.2 Insulin and glucagon as control inputs

In the previous section we tuned the weights α and α inside the objective function (2). We now consider the case that u > 0 and we tune α, the weight associated with the glucagon expenditure in the objective function (2), by keeping α = 10, α = 1 for p = 1 and α = 103, α = 1 for p = 2, as previously determined.

In Fig 4, we plot the optimal Δ, Gmin and Gmax as functions of the parameter α, respectively. A large value of α indicates that we are placing a large weight on the expenditure of glucagon within the objective function (2), i.e., the larger the value of α, the less glucagon we use. By looking at Fig 4, we observe that the values of Δ decrease as α decreases, i.e., we can obtain lower (improved) values of Δ if we allow for a larger expenditure of glucagon. We note that past a certain value of α (10−2 in the both the ReMF and ReME problems) no further reduction in Δ is observed. As in the previous case, the maximum glucose level Gmax, shown in Fig 4, improves (decreases) when Δ improves (decreases). Interestingly, different from the previous case, also the minimum glucose level Gmin (Fig 4) improves (increases) with Δ and Gmax: this is a consequence of the fact that we are using both insulin and glucagon as control inputs, which enables us to avoid both hypoglycemia and hyperglycemia.

Fig 4

Performance of the optimal control solution as a function of α.

A) Δ vs. α. B) Gmin vs. α. C) Gmax vs. α. D) ϕ vs. ϕ. E) Δ vs. ϕ. We select α = 10−2 for both of the REMF and REME problems, which are indicated by dashed circles in the figure.

In Fig 4 we plot the projection of the Pareto front in the (ϕ, ϕ) plane. By looking at the Fig 4, ϕ and ϕ appear to be positively correlated and related by an approximately linear relation. While the timing of administration of insulin and glucagon is different, we see that overall the more insulin is used in the optimal solution, the more glucagon is used as well. This is because the two hormones have opposite effects in the regulation problem and thus they work so as to balance each other. This is also consistent with the observation that with the dual drug therapy (insulin and glucagon) it becomes possible to simultaneously improve Δ, Gmin, and Gmax. From the data in Fig 4 we derive the following approximate linear relationship between ϕ and ϕ, Obviously, glucagon should be used only when ϕ(ϕ) > 0.

Panel Fig 4 shows a projection of the Pareto front on the (ϕ, Δ) plane. We see again that the ReMF front dominates the ReME one, i.e., a pulsatile therapy gives better results than a continuous therapy in terms of Δ and also uses lower amounts of the two drugs (smaller ϕ, and thus smaller ϕ due to the positive correlation found in Fig 4).

The Pareto front is monotonically decreasing in Fig 4 which indicates a trade-off between the total amount of drugs used and the achievable glucose control performance. We choose the value of α for which the ratio between the increase in Δ and the decrease in ϕ is minimized, i.e., α = 10−2 for both ReMF and ReME problems, which are indicated by dashed circles in the figure.

Fig 5 show the results of the optimal control problem for α = 10, α = 1 and α = 10−2 when p = 1; and α = 103, α = 1 and α = 10−2 when p = 2. In Fig 5 we plot the time evolution of glucose G(t) for the different optimal solutions. The blue curve corresponds to the solution of the ReMF problem when only insulin is used (the blue curve in Fig 3). The red and green curves correspond to the solution of the ReMF and the ReME problems for the dual drug therapy. We observe that G(t) reaches the desired level G faster if we use both insulin and glucagon as control inputs, compared to the case that only insulin is used. We also see that in this case both Gmax decreases and Gmin increases. We therefore conclude that the therapy with both insulin and glucagon performs better than the therapy with only insulin, as the risks for both hypoglycemia and hyperglycemia are reduced and glucose fluctuations are suppressed.

Fig 5

A) Time evolution of glucose G(t) (in mg/dL). The blue curve corresponds to u(t) obtained by solving the ReMF problem. The red curve corresponds to u(t) and u(t) obtained by solving the ReMF problem using the dual drug therapy. The green curve corresponds to u(t) and u(t) obtained by solving the ReME problem using the two-drug-therapy. B) Time evolution of the insulin infusion rate u(t) (in mg/dL). Color code is consistent with A. C) The cumulative insulin supply r(t) as a function of time t. D) Time evolution of the glucagon infusion rate u(t) (in mg/dL). E) The cumulative glucose supply r(t) as a function of time t.

In Fig 5 we plot the optimal insulin infusion rates and in Fig 5 we plot the cumulative insulin supply r(t) as a function of time t. We observe that for the ReMF problem, the pulse in insulin appears at t = 32 minutes in the case that both insulin and glucagon are used (28 minutes before the meal), whereas the pulse appears at t = 20 minutes when only insulin is used. From Fig 5, we see that, for the ReMF problem, the glucagon delivery function is pulsatile with a main pulse appearing at t = 145 min (one hour and 25 minutes after the meal) and a secondary pulse appearing at t = 203. The dual drug therapy shows a noticeable difference between the ReMF solution and the ReME solution. As expected, the solutions of the ReME problem are continuous. The glucose response to the ReME therapy is better than the glucose response to the ReMF solution. Specifically, the green curve has smaller oscillations (in panel A) at the cost of small increases in the total amounts of used insulin and glucagon (compare panels C and E)).

Based on the results in Fig 5, we propose a possible ad-hoc dual drug therapy to be used as an alternative to the standard therapy. Rather than administering insulin half an hour before the meal (standard therapy), better glucose regulation can be achieved with a slightly larger insulin injection half an hour before a meal followed by a glucagon injection one hour and thirty minutes after a meal. The insulin injection of the ad-hoc dual drug therapy is 25% larger than the one used in the standard therapy, which is consistent with the relation between ϕ for the monotherapy ReMF optimal solution and the one used in the dual drug therapy.

In Fig 6 we present a comparison between the glucose response to the standard insulin base therapy (orange curve) and the proposed ad-hoc dual therapy (cyan curve) for the case of a meal with 70 grams of glucose (for the particular patient considered this corresponds to 10 units of insulin half an hour before the meal) and the proposed ad-hoc dual drug therapy (which consists of 12.43 units of insulin thirty minutes before the meal and 0.40 mg of glucagon one hour and thirty minutes after the meal). We observe that the ad-hoc dual drug therapy performs better in terms of all of the proposed measures (Δ, Gmin, Gmax, ϕ and ϕ) as opposed to the standard insulin based therapy.

Fig 6

Comparison between the glucose response to the standard insulin base therapy (orange curve) and the proposed ad-hoc dual therapy (cyan curve).

3.3 Robustness analysis

We now analyze the robustness of the optimal control therapies we have proposed with respect to model parameter mismatches, which is a fundamental step for implementation of model based control. We consider two different types of mismatches. The first type accounts for variability in the patient’s behavior, in terms of both the time of the meal τ and the amount of glucose intake D. The second type accounts for deviations in the parameter estimation, as well as the temporal variability of the parameters that a patient may experience during the day [4].

3.3.1 Robustness against variability of the meal time and glucose intake

In this section we analyze the robustness of the optimal ReMF therapies (both monotherapy and dual therapy) with respect to the two “control” parameters the patient has. The first one is the variation in the meal time, () (min), where is the time of a meal and τ is the time of a meal we assumed in order to compute the optimal therapies. The second one is the variation of glucose in the meal, (), where is the glucose intake in a meal and D is the glucose intake we assumed to compute the optimal therapies. We consider variations in the meal time in the interval [30, 90] min and variations of the glucose intake in the interval [40, 100] g.

The results of this study are illustrated in Fig 7. Fig 7 provides a visual assessment of the quality of the optimal therapies in terms of the three proposed measures Δ, Gmax and Gmin (the over-bar stands for evaluation at the perturbed parameter values ). The color in Fig 7 varies according to the control performance from green (good) to red (dangerous). In the upper panels (A–C) we consider the optimal ReMF monotherapy, while in the lower panels (D-F) we consider the optimal ReMF dual thearpy. Cross symbols indicate the application of the optimal control therapies under ideal condition, i.e., when and . The black curves labeled by 180, 90 and 70 in Fig 7 are the curve level plots for , and , respectively. The black curves labeled by 180 in Fig 7, are the curve level plots for .

Fig 7

Robustness of the optimal control solution against variations in the meal timing and the amount of glucose in the meal.

A)–C) show the results obtained for the ReNF problem (p = 1) with only insulin provided, D)-F) ReMF (p = 1) problem with both insulin and glucagon provided. Cross symbols indicate the application of the optimal control therapies for and . The blue cross symbols correspond to the optimal therapies for the ReMF problem with only insulin. The red cross symbols correspond to the optimal therapies for the ReMF problem with both insulin and glucagon. A) and D) are plots of in the control parameters space . B) and E) are the plots of in the control parameters space . C) and F) are the plots of in the control parameters space .

We see from Fig 7 that the optimal therapies for the ReMF problem (using only insulin or both insulin and glucagon) are robust with respect to variations in the control parameters: remains well bounded in most of the considered parameter space. In particular we see from Fig 7 that the proposed optimal therapies are robust against hyperglicemic events: for example, even if exceeds D by 50% and exceeds τ by 30 minutes, the patient will not enter the hyperglycemic regime (Gmax > 300). Fig 7 reveal that the proposed therapies suffer from a certain lack of robustness with respect to hypoglycemic events (Gmin < 70), the most dangerous ones. The dangerous cases are, however, confined to extreme situations in which and minutes. Fig 7 show also that the optimal therapy for the ReMF problem with both insulin and glucagon is more robust (larger green region and smaller yellow region) than the optimal therapy for ReMF problem with only insulin (smaller green region and larger yellow region): thus the use of glucagon alleviates the risk of severe, life-threatening hypoglycemia.

We obtain qualitatively similar results when performing the same analysis for the other therapies we proposed (the ReME therapies and the ad-hoc dual drug therapy).

3.3.2 Robustness to parameter mismatches

We consider perturbation of the model parameters up to 20% of their nominal values, where φ is a random number from a normal distribution , Θ is a nominal parameter for a given patient with basal glucose level G and represents the associated perturbed parameter. We then apply the optimal insulin and glucagon dosing, calculated for the unperturbed system, to 100 perturbed systems. This is analogous to testing the computed optimal control therapy on a specific patient, but the patient’s parameters may vary due to imperfect knowledge or due to the parameter variability throughout the day. The results of this study are illustrated with a Control Variability Grid Analysis (CVGA), see Fig 8. The CVGA provides a simultaneous visual and numerical assessment of the overall performance of the glycemic control strategies in terms of the achieved minimum/maximum glucose values in the space of parameters mismatches. In Fig 8, points in the light green region indicate accurate blood glucose control while points in the dark green regions indicate the patient is not immediately at risk of either hypoglycemia or hyperglycemia. Points in the top two yellow/orange regions indicate an elevated risk of hyperglycemia and points in the the right two yellow/orange regions indicate an elevated risk of hypoglycemia. Finally, points in the red corner region indicate an elevated risk of both hyperglycemia and hypoglycemia. Each point reported in the figure is a plot of Gmax vs. Gmin. Here, the black dots correspond to the glucose response when a certain therapy is applied to a system with perturbed parameters. Cross symbols indicate application of the optimal control therapies to the unperturbed systems.

Fig 8

Robustness of the optimal control solution against parameter perturbations of the system and CVGA in the G, G plane.

The analysis is performed for A) ReMF (p = 1) problem with only insulin provided, B) ReME (p = 2) problem with only insulin provided, C) ReMF (p = 1) problem with both insulin and glucagon provided, D) ReME (p = 2) problem with both insulin and glucagon provided, E) the standard therapy, and F) the proposed ad-hoc dual drug therapy. Cross symbols indicate the application of the optimal control therapies to the unperturbed systems.

For the of monotherapy (ReMF in Fig 8 and ReME in Fig 8) we find that the control is 67% and 61% accurate, respectively. For the dual therapy case (ReMF in Fig 8 and ReME in Fig 8) we find the control is more accurate than for the case of monotherapy, 92% and 94% accurate, respectively. The least robust control is obtained with the standard therapy (shown in Fig 8), attaining only 37% accuracy. Note that the optimal dual drug therapies (Fig 8) are not only more robust than the optimal insulin therapies (Fig 8), but also than the standard therapy (Fig 8). We also see that the ad hoc therapy (Fig 8) is more robust than the standard therapy (Fig 8).

4 Discussion

In this paper we have used the Glucose-Insulin-Glucagon mathematical model proposed in [2-4], which describes how the body responds to exogenously supplied insulin and glucagon in patients affected by Type I diabetes and designed an optimal dosing schedule of either insulin or insulin and glucagon together to regulate the blood glucose index (BGI), while limiting the total amount of insulin and glucagon administered. The numerical optimal control software has been used to solve this optimal control problem. While the numerical solution requires knowledge of the set of model parameters, which are patient specific, the solutions we obtain provide insight into the best possible glucose regulation with insulin or with insulin and glucagon together. Our approach is in agreement with the results of references [60-62], in which simplified models are used to analytically establish general theoretical properties and control limitations for the glucose regulation problem.

Two distinct regulation problems have been considered: the minimum fuel problem (ReMF) which yields pulsatile (shot-like) type solutions and the minimum energy problem (ReME) which yields longer periods of time over which insulin is administered but with smaller delivery rates. This has allowed us to compare standard therapies which typically consist in shots of insulin with therapies in which insulin is delivered continuously. In [63, 64] it has been proven that the optimal control is pulsatile when the aim of the control is to minimize the variation in the maximum and minimum output response, the system is positive (like the one we are considering) and the disturbance (the meal, in our case) is pulsatile. Our work indicates that pulsatile control is still a good choice when more complex objective functions are chosen. Moreover, a pulsatile control appears to be optimal for alternative more realistic models of the meal (for example, a meal that is consumed over a window of 15 minutes). We also see that a continuous drug delivery can achieve better results in the case of the dual therapy, thus pointing out the importance of developing a commercial pump able to deliver both insulin and glucagon.

For both the ReMF and ReME problems, we compute the optimal drug dosing schedules when only insulin is available and when both insulin and glucagon are available. The solution of the insulin only ReMF problem, astoundingly, is nearly equal to the standard method of insulin based glucose regulation. Similarly, the solution of the ReMF problem when insulin and glucagon are used is also pulsatile, except that the amount of insulin administered is larger and the administration time is closer to the time of the meal, while the glucagon is mostly delivered in a shot about an hour and thirty minutes after the meal.

The solution for the ReME problem when insulin only is available as well as when both insulin and glucagon are available is different from the ReMF solution in that, rather than being pulsatile, insulin and possibly glucagon are delivered at a slower rate over longer periods of time. Nonetheless, the total amount of insulin and possibly glucagon is about the same, and the peak of the longer delivery time occurs approximately at the same time of the shot according to the solution of the ReMF problem. The obtained glucose profiles for the optimal ReMF and ReME problem solutions do not differ too much from each other: taken together, these results indicate that the amounts of insulin and glucagon, and the peak times of delivery, are the most important factors to determine when computing the optimal solutions.

Based on the above results, we have proposed the following ad hoc therapy when insulin and glucagon are used in combination: Administer a shot of insulin (with 5% more insulin than the amount required by the standard therapy based on the planned meal) 30 minutes before eating. Administer a shot of glucagon of an amount specified by Eq (18) one hour and thirty minutes after completing the meal. This therapy must be used with caution as the amount of insulin injected can lead to hypoglycemia if the shot of glucagon is not administered as well.

All optimal dosing schedules we computed were tested for robustness with respect to variations in the meal timing and size and with respect to variability of the parameters. The therapies we proposed typically maintain the patient in the healthy region even under variable conditions and patient behavior. Note, however, that the proposed therapies are open-loop (the drug schedule is computed only from the condition of the patient at the initial time), thus cannot compensate for unexpected behavior that can arise due to modeling simplifications (e.g. we do not consider how physical activity influences the blood glucose production and consumption [65, 66]), measurement noise or bias. A step towards the real application of our methodology is a real-time closed-loop strategy; this is possible, since the typical time needed to compute an optimal solution on a standard laptop (i7-8550U CPU with 16GB RAM) is around 2 minutes. Another main limitation of our study is that real life constraints and long term physiological effects may make a therapy based on exogenous administration of both insulin and glucagon impractical.

Our optimal control strategies require knowledge of the meal time and meal glucose amount. This is somewhat undesirable, as recent advances in diabetes therapy have moved towards devices that do not require the user to provide information about the meals. Our results emphasize the importance of knowing when the meals will occur, and that new dosing schedules would benefit from some knowledge about the meals. This seems to indicate that it would be beneficial to provide the pump with the ability to interpret the patient’s behavior.

GIGM model and parameters.

The nonlinear equations of the GIGM model are given in Eqs (S1)-(S9). The parameters are given in S1 Table. The equations are given in Eq (S10) and the basal steady states are given in Eq (S11). The basal values are given in S2 Table. The references [67-69] are appeared here.

(PDF)

Click here for additional data file.

Network representation of the GIGM model.

The network representation of the GIGM model. The references [70-85] are appeared in this section.

(PDF)

Click here for additional data file.

Smooth approximation of non-differential functions in ODEEs.

There are some non-differential functions in Eqs (S1)-(S9). We present the smooth approximations of these non-differential functions in Eqs (S12)-(S15).

(PDF)

Click here for additional data file.

Average parameters.

Average Parameters of the GIGM model.

(PDF)

Click here for additional data file.

Basal values.

Basal values of the GIGM model.

(PDF)

Click here for additional data file.

(EPS)

Click here for additional data file.