Optimality of Unconstrained Pulse Inputs to the Bergman Minimal Model

TL;DR

This paper characterizes the optimality of rectangular bolus insulin inputs in the Bergman minimal model, ensuring plasma glucose levels stay above a minimum, with numerical illustrations using the Hovorka model.

Contribution

It provides a theoretical framework for determining optimal insulin bolus inputs under specific constraints in glucose regulation models.

Findings

Optimal bolus inputs are rectangular under given constraints.

The model ensures plasma glucose remains above a minimum level.

Numerical examples demonstrate the theoretical results.

Abstract

We characterise optimality of bolus insulin inputs, to the Bergman minimal model, by the predicted behaviour of the plasma glucose concentration for a given disturbance. The result is derived subject to the constraints that the plasma glucose concentration must attain but not go below a specified minimum value and the bolus input is rectangular. We give numerical examples of the results for the Hovorka model.

Equations

$$\begin{array}{rl}\dot{z}& =-dz+dku\\ \dot{y}& =-cy+cz\\ \dot{x}& =-ax+aby\\ \dot{g}& =-hg+w\end{array}$$

$$\begin{array}{r}h=x+G\\ w=r+E\end{array}$$

$$u(t,A)=\bar{u}+\hat{u}{\chi }_A(t)$$

$$\overline{u}=\frac{1}{kb}\left(\frac{E}{g(0)}-G\right)$$

$$L_i:=$$

$$\cdots ,\frac{k_i{\sigma }_{i-1}+(n-k_i){\alpha }_{i-1}}{n},\cdots ,{\alpha }_{i-1}\}$$

$$g\left(\frac{k_i{\sigma }_{i-1}+(n-k_i){\alpha }_{i-1}}{n}\right)$$

$$g\left(\frac{k_j{\sigma }_{i-1}+(n-k_j){\alpha }_{i-1}}{n}\right)$$

$${\sigma }_{i+1}=\frac{k_i{\sigma }_i+(n-k_i){\alpha }_i}{n}$$

$${\alpha }_{i+1}=\frac{(k_i-1){\sigma }_i+(n-k_i+1){\alpha }_i}{n}$$

$${\alpha }_{i+1}-{\sigma }_{i+1}=\frac{1}{n}\left({\alpha }_i-{\sigma }_i\right)=\dots =\frac{1}{n^{i+1}}\left({\alpha }_0-{\sigma }_0\right)$$

$$g(u_{i+1}){|}_B\ge g(u_i){|}_B$$

$$\left(\begin{array}{c}{\dot{f}}_1(t)\\ {\dot{f}}_2(t)\end{array}\right)=\left(\frac{1}{60}\right)\left(\begin{array}{cc}-1& 1\\ 0& -1\end{array}\right)\left(\begin{array}{c}{f}_1(t)\\ f_2(t)\end{array}\right)+\left(\begin{array}{c}0\\ \rho (t)\end{array}\right)$$

$$g=f(g,u,h,w,t)$$

$$\begin{array}{rl}\dot{z}& =-dz+u\\ \dot{y}& =-dy+dz\\ \dot{x}& =-kx+cdy\\ {\dot{x}}_1& =-a_1{x}_1+a_1{b}_{1}x\\ {\dot{x}}_2& =-a_2{x}_2+a_2{b}_{2}x\\ {\dot{x}}_3& =-a_3{x}_3+a_3{b}_{3}x\\ {\dot{q}}_1& =-h_1{q}_1+l{q}_2+w\\ {\dot{q}}_2& =-h_2{q}_2+x_1{q}_1\end{array}$$

$$\begin{array}{rl}w& =E+r\\ h_1& =V^{-1}\left(f_c+f_r+\frac{x_{3}E}{g}\right)+x_1\\ h_2& =l+x_2\end{array}$$

$$\begin{array}{rl}{f}_r& :=\{\begin{array}{ll}VR\left(1-\frac{{\overline{g}}_r}{g}\right),& g\ge {\overline{g}}_r\\ 0,& \text{otherwise}\end{array}\\ f_c& :=\{\begin{array}{ll}f{g}^{-1},& g\ge {\overline{g}}_c\\ f{\overline{g}}_c^{-1},& \text{otherwise}\end{array}\end{array}$$

$$\begin{array}{r}g=V^{-1}{q}_1\end{array}$$

$$\left(\begin{array}{c}{\dot{f}}_1(t)\\ {\dot{f}}_2(t)\end{array}\right)=\left(\frac{1}{55}\right)\left(\begin{array}{cc}-1& 1\\ 0& -1\end{array}\right)\left(\begin{array}{c}{f}_1(t)\\ f_2(t)\end{array}\right)+\left(\begin{array}{c}0\\ 0.8\rho (t)\end{array}\right)$$

Full text

I Introduction

Type one diabetes is a chronic disease affecting over thirty-eight million people [1]. Diabetics, typically, require the subcutaneuous administration of insulin to minimise plasma glucose concentrations whilst avoiding hypoglycaemia. Current treatment is invasive and often provides poor control. Hence, much recent effort has been devoted to developing an artificial pancreas [2] to automate treatment and better control plasma glucose concentrations.

Understanding and modelling the dynamics of glucose regulation assists the development of such systems and further treatment improvements. A number of models of glucose regulation have been proposed ([3, 4, 5]). Each is typically comprised of sub-systems describing different physiological processes such as insulin kinetics and glucose absorption.

Recently, research has focused on comprehensive models of glucose dynamics which are generally preferred to test treatment policies and control algorithms, for example [6]. Typically, these models are high order dynamic system with many parameters to ensure robustness to inter-individual variability. However, simpler models are useful to establish general theoretical properties that would otherwise be difficult to investigate analytically. Indeed, most models of glucose dynamics share certain analytic properties – such as positivity the of the plasma glucose. Thus analytic results obtained for simpler models can give insights into the behaviour of more comprehensive models.

We focus here on the Bergman (Khandarian) Minimal Model ([7, 8, 9]) which is a simplified model of glucose metabolism frequently used for virtual patient simulations and as the basis of more comprehensive models such as the Fabietti model ([10]) and the extensions of [11] and [12]. The model (1) is a non-linear continuous-time model comprising a set of first order linear ordinary differential equations which govern the subcutaneous, plasma and interstitial concentrations and effectiveness of insulin and a non-linear ordinary differential equation which governs the plasma glucose concentration $g(t)$:

[TABLE]

where all variables and constants are positive and $u(t)$ is the input function. The functions $h$ and $w$ in (1) are:

[TABLE]

where the function $r$ is a given bounded function. Specifically, the terms in (1) and (2) represent:

• •

$u(t),z(t),y(t)$ and $x(t)$ – the delivery, subcutaneuos concentration, plasma concentration and insulin effectiveness, respectively.

• •

$c$ and $d$ – inverse time constants.

• •

$a,b$ and $k$ – the insulin motility [12], insulin sensitivity and the clearance rate.

• •

$g(t)$ – the plasma glucose concentration.

• •

$E$ and $G$ – the endogenous glucose production and the effect of glucose on the uptake of plasma glucose and suppression of endogenous glucose production.

• •

$r(t)$ – the glucose absorption from meals.

Physiological values for the above are derived from [9] and given in Table 1 of [8].

We contribute to the theoretical understanding of this model by characterising the magnitude, delivery time and duration of insulin bolus inputs that are optimal in the sense that they give the lowest maximum glucose concentration whilst avoiding hypoglycaemia, see Sections III, III and III. Specifically, we impose a fixed constraint on the minimum glucose concentration and focus on lowering the maximum glucose concentration. We show that this fixed constraint induces a fundamental limitation on the controllability of the maximum glucose concentration when the control input is a pulse. We constrain the minimum glucose concentration because the risks associated with hypoglycaemia are, generally, far greater than those associated with hyperglycaemia. To ensure robustness against uncertainties this constraint could be set above the hypoglycaemic threshold.

The effect of a fixed constraint on the control of plasma glucose concentrations has been investigated in [8]. The authors consider a discretised non-linear model, derived from the Bergman model, which they use to derive a non-linear insulin bolus dosing algorithm. However, the bolus is constrained to be an impulse applied contemporaneously with an impulsive food input. This is a specific example of the cases considered here. For instance, in [8], the duration $\tau=0$ is fixed and $w$ is assumed to be the response of a second order system to a single impulse.

In [13] a pulse input $u(t)$ of fixed duration was shown to be optimal, i.e. it minimised the global maximum glucose concentration, if and only if either the fixed minimum glucose concentration occured between two global maxima of $g(t)$ or the global maximum occured between two fixed minima of $g(t)$. Here, we present the counterpart of these results by giving conditions on inputs of varying durations but fixed delivery time to minimise the global maximum glucose concentration. Furthermore, our main contribution is to generalise the results to pulse inputs of any duration and delivery time. This fully characterises optimality of arbitrary pulse inputs in the sense of minimising the maximum glucose concentration subject to a fixed constraint on the minimum glucose concentration.

The observations of this work are that: firstly, distinct optimal inputs – in the sense of [13], see Sections III and III – must intersect at least twice if one has a lower global maximum glucose concentration and, secondly, that pulse inputs of varying duration can intersect at most twice and will only do so if one input is nested inside the other. Our results confirm the intuition that responses with a maximum between two minima result from longer pulses than responses with a minimum between two maxima. Finally, decreasing the duration for the first type of response or increasing the duration the second type of response will lower the global maximum. As $g$ is a continuous function of the duration the lengthening and shortening of the duration converges.

The presented results reveal a fundamental limit on the controllability of the plasma glucose concentration achievable from a bolus input to the Bergman minimal model and allow the optimality, of the input, to be determined from the shape of the glucose response. They also specify the effect of changes to the parameters of a bolus input on the maximum plasma glucose concentration. This may, for example, act as metric for the optimality of control algorithms designed for artificial pancreas systems and assist in the determination of bolus guidelines. Regardless of our focus on the Bergman model, other models may be analysed mutatis mutandis, see Section V.

Notation:

We adopt the following notation throughout: $\overline{u}$ and $\hat{u}$ are the basal input and the magnitude of the bolus input; $\lambda$ and $\gamma$ are the global minimum glucose concentration and the global maximum glucose concentration; $t^{\prime},t_{i,\max},t_{i,\min}$ and $\tau$ are the delivery time, the $i^{\text{th}}$ time when the glucose concentration is at its global maximum, the $i^{\text{th}}$ time when the glucose concentration is at its minimum and the duration of the interval over which the bolus is delivered; $u(t,\tau)=u(t,A)$ is the input $u(t)$ applied over the interval $A:=[t^{\prime},t^{\prime}+\tau]$; $g(h(u),w)=g(t,\tau)$ is the reponse of $g$ to the functions $h$ and $w$, where $h(u)$ is the response of $h$ to the input $u(t,\tau)$; $t_{i}$ and $t_{g,i}$ are intersection points of the responses $h(u)$ and $h(v)$ for distinct inputs $u\neq v$ and the $i^{\text{th}}$ intersection point of the resulting $g(h(u),w)$ and $g(h(v),w)$ and, lastly, $\gamma(u)$ is the global maximum of $g(h(u),w)$.

II Assumptions and Preliminaries

Regardless of the nominal defintions given above, we do not require $r$ to be a positive bounded function corresponding to the glucose absorption from meals nor $E$ to be the endogenous glucose production. Instead, we require that $w$ is a positive function bounded below by any positive real $\overline{E}\geq G$. This allows, for example, $r$ to be negative if $\overline{E}<E$. By abuse of notation we denote $\overline{E}$ by $E$.

Throughout we impose the following initial conditions: $z(0)=y(0)=ku(0)$, $x(0)=bku(0)$ and $g(0)>0$. We assume the function $w$ is positive and bounded. We also assume the input $u$ is positive and bounded and of the form:

[TABLE]

where the constant $\bar{u}$ is the basal input, $\hat{u}$ is the magnitude of the bolus input applied at some time $t^{\prime}$, known as the delivery time and $\chi_{A}$ is the characteristic function of $A$. The bolus input is held constant over $A=[t^{\prime},t^{\prime}+\tau]$, where $\tau\in\mathbb{R}_{+}$. When $\tau=0$ we define $u(t):=\hat{u}\delta_{t^{\prime}}(t)$, where $\delta(t)$ is the Kronecker delta. The boundedness and positivity of $u(t)$ imply that $h$, given by (1) and (2), is a continuous, positive and bounded function. We desire that there exist $\lambda>0$ such that $g(t)\geq\lambda$ for all $t$. This is achieved if $\lambda$ is a global minimum of $g(t)$. We denote by $t_{min}\in\mathbb{R}_{+}$ a point such that $g(t_{min})=\lambda$.

Finally, unless otherwise stated we assume that $t_{\max}:=\operatorname*{arg\,max}_{t}{g(t)}<\infty$. The maximal time $t_{\max}$ exists if $w$ is assumed to vanish to its lower bound at infinity. 1 summarises a number of useful results from [13].

Definition \thedefn (Steady-State)

The steady-state of $g$ is $g(\infty):=\lim_{t\to\infty}g(t)$, when $\lim_{t\to\infty}u(t)=\overline{u}$ and $\lim_{t\to\infty}w(t)=E$ i.e. it is the limit of the response of $g(t)$ when the only input is the constant input $\overline{u}$.

Theorem 1** (Portmanteau)**

Suppose $h$ and $w$ are bounded positive real-valued functionals, $g$ is as in (1), $u(t)$ is as in (3) and choose $\lambda\leq g(0)$ and $\tau\geq 0$. Then:

Under the assumed initial conditions, $x(\overline{u})=bk\overline{u}$ for all $t$. Furthermore, $\lim_{t\to\infty}x(t)=bk\overline{u}$. 2. 2.

$g(t)$* is a strictly monotone function of $u(t,\tau)$.* 3. 3.

Setting:

[TABLE]

gives $g(\infty)=g(0)$.

Definition \thedefn (Proper Input)

For some $\lambda\leq g(0)$, an input function, $u(t,\tau)$, is proper, if there exists $t_{\min}$ such that $g(h(u(t_{\min})),w)=\lambda$ and $g(t)\geq\lambda$ for all $t$.

Theorem 2 (Theorem 7 [13]) proves the existence of a bolus input delivered at any $t^{\prime}$ and $\tau$ which achieves a specified minimum $\lambda>0$ and thus proves the existence of proper inputs of the form (3).

Theorem 2** (Insulin Bolus, (Theorem 7 [13]))**

Suppose $u(t)$ is of the form (3). Fix $\tau$ and $t^{\prime}$ – the input time i.e. $A:=[t^{\prime},t^{\prime}+\tau]$, choose $\lambda\in(0,g(t^{\prime})]$ and suppose $\bar{u}$ is as in (4). Then there exists $\hat{u}$ such that $u(t)$ is proper.

III Optimal Duration

Two necessary and sufficient conditions were given in [13] for the response $g(t)$ to an input of the form (3) with a fixed duration $\tau$ to be optimal.111As in [13] we say an input is optimal if it results in the lowest maximum of $g(t)$ for all inputs of the same duration. These conditions are summarised in Definitions III and III.

Definition \thedefn ($\lambda$–optimal)

An optimal input $u$ is $\lambda$–optimal if the global maximum of $g(u)$ occurs between two minima.

Definition \thedefn ($\gamma$–optimal)

An optimal input $u$ is $\gamma$–optimal if all minima of $g(u)$ occur between two global maxima.

For a fixed $w$, if the input $u$ is $\lambda$–optimal, respectively if $u$ is $\gamma$–optimal we say the response $g(u)$ is $\lambda$–optimal, respectively $\gamma$–optimal. We extend the results of [13] to inputs of the form (3) which may have any duration $\tau$. In Section III we define global optimality of an input. In this section, firstly, we consider global optimality over the class of $\lambda$–optimal inputs and secondly global optimality over the class of $\gamma$–optimal inputs. Finally, we characterise global optimality over all proper inputs of the form (3). The following defines some useful notation.

**Definition **\thedefn

Let $\gamma(u,w):=\max_{t}\{g(u(t),w(t))\}$. This is denoted $\gamma(u)$ for fixed $w$.

Definition \thedefn (Globally Optimal)

An input $u$ is globally optimal if $\gamma(u)<\gamma(v)$ for all $v\neq u$.

Lemma III specifies the maximum number of intersection points of responses to distinct inputs of the form (3).

Definition \thedefn (Nested Inputs)

Two inputs $u$ and $v$ are nested if $A\subset B$, where $A$ and $B$ are the intervals over which the boluses $\hat{u}$ and $\hat{v}$ are applied.

Throughout we adopt the convention that for two inputs $u$ and $v$ times related to $u$ are denoted by $t$ and times related to $v$ are denoted by $s$.

**Lemma **\thelem

Suppose $u$ and $v$ are distinct inputs with delivery times $t^{\prime}$ and $s^{\prime}$ respectively. Then, for each solution $\phi\in\{x(t),y(t),z(t),h(t),g(t)\}$ to (1)–(2), there are at most two $t_{i}>\min\{t^{\prime},s^{\prime}\}$ such that $\phi(u,t_{i})=\phi(v,t_{i})$ and these $t_{i}$ are distinct for all $\phi$ only if $u$ and $v$ are nested.

Proof:

Observe, for $t>\min\{t^{\prime},s^{\prime}\}$ that $z(u,t)=z(v,t)$ only if $u-v$ changes sign. As $u-\overline{u}$ and $v-\overline{u}$ are rectangular $u-v$ can change sign at most twice. Implying that $z(u,t)=z(v,t)$ at most twice and thus $z(u,t)-z(v,t)$ may change sign at most twice. We proceed similarly for all solutions, $\phi$. ∎

III-A $\lambda$–Optimal Inputs

**Lemma **\thelem

Suppose $g(t,\tau)$ and $g(t,\sigma)$ are the respective responses to the distinct $\lambda$–optimal inputs $u(t)$ and $v(t)$ with durations $\tau$ and $\sigma$, respectively. Then $\gamma(u)>\gamma(v)$ if and only if $[s^{\prime},s^{\prime}+\sigma]\subset[t^{\prime},t^{\prime}+\tau]$, where $t^{\prime}$ and $s^{\prime}$ are the respective delivery times of the inputs $u(t)$ and $v(t)$.

Proof:

As $u$ and $v$ are pulse inputs there are at most two points $t_{i}$ at which the response $h(u)$ and $h(v)$ intersect, and similarly at most two $t_{g,i}$. We denote these $t_{i,1}$ and $t_{i,2}$, similarly $t_{g,1}$ and $t_{g,2}$.

Suppose $\gamma(u)>\gamma(v)$. Then $g(v)<g(u)$ for all $t$ in some interval $I\ni t_{1,\max}$. Additionally, $g(v)\geq g(u)$ at both $t_{1,\min}$ and $t_{2,\min}$. Thus $g(v)>g(u)$ for all $t\in(\min\{s^{\prime},t^{\prime}\},t_{g,1})$, where $t_{g,1}\geq t_{1,\min}$. Otherwise there would exist more than two intersection points or $v$ would be non-proper. Therefore, $t^{\prime}<s^{\prime}$. Lastly, observe that should $s^{\prime}+\sigma>t^{\prime}+\tau$ then by Section III there is at most one $t_{i}>t^{\prime}$ such that $h(u)=h(v)$ which implies that either $u$ or $v$ is not $\lambda$–optimal.

Suppose instead that $[s^{\prime},s^{\prime}+\sigma]\subset[t^{\prime},t^{\prime}+\tau]$. As both $u$ and $v$ are $\lambda$–optimal, then by the above, $[s_{1,\min},s_{2,\min}]\subseteq[t_{1,\min},t_{2,\min}]$ or the converse. By assumption on the inputs we have that $g(u)<g(v)$ for all $t\in(t^{\prime},t_{g,1})$. If $s_{1,\min}<t_{1,\min}$. Then we have that $t_{g,1}<t_{g,2}\leq t_{1,\min}$. Thus $g(v)$ would not be $\lambda$–optimal. As $g(v,t_{1,\max)}>\gamma(u)$ and $g(v)>\lambda$ for all $t>t_{g,2}$ which must occur before $t_{1,\min}$. Thus $[s_{1,\min},s_{2,\min}]\subseteq[t_{1,\min},t_{2,\min}]$. Finally, as $g(v)$ and $g(u)$ intersect at most twice we have that $\gamma(u)>\gamma(v)$. ∎

**Corollary **\thecor

Suppose $g(v)$ is $\lambda$–optimal for all $\tau>0$. Then $u$ is globally optimal if and only if $A$ is a singleton.

III-B $\gamma$–Optimal Inputs

**Lemma **\thelem

Suppose $u$ and $v$ are distinct inputs for which there exists unique $t_{i}$ such that $h(u)=h(v)$ and $h(u)>h(v)$ for all $t\neq t_{i}$ and $t>\min\{s^{\prime},t^{\prime}\}$. Then there are two distinct $p_{i}$ such that $y(u,p_{i})=y(v,p_{i})$, where $y(u,t)$ is the response of $y(t)$, from (1) to the input $u$.

Proof:

Note that $t_{i}$ must be a minimum of the non-negative function $f(t):=x(u)-x(v)$. Thus there exists $\varepsilon>0$ such that $x^{\prime}(u)<x^{\prime}(v)$ for $t\in(t_{i}-\varepsilon,t_{i})$ and $x^{\prime}(v)<x^{\prime}(u)$ for $t\in(t_{i},t_{i}+\varepsilon)$. By assumption $x(v)\leq x(u)$ for all $t$ thus, from (1), $f^{\prime}(t)$ can only change sign about $t_{i}$ if $y(v)-y(u)$ changes sign about $t_{i}$. Hence, there is some $t_{y}<t_{i}$ such that $y(u)<y(v)$ and similarly there is some $s_{y}>t_{i}$ such that $y(v)<y(u)$. By continuity of $y(t)$ we see that there is a $p_{2}\in(t_{y},s_{y})$ such that $y(u)=y(v)$. Additionally, $f^{\prime}(t)>0$ in some non-empty interval $[\min\{t^{\prime},s^{\prime}\},\min\{t^{\prime},s^{\prime}\}+\delta)$ as $h(u)>h(v)$, for almost all $t>\min\{t^{\prime},s^{\prime}\}$. We have that $y(v)<y(u)$ on this interval. This implies that there must exist $p_{1}<t_{y}$ such that $y(v)=y(u)$, again by continuity. ∎

**Corollary **\thecor

Suppose there are at most countably many $t_{i}$ such that $h(u)=h(v)$ and that $h(u)>h(v)$ for all $t\neq t_{i}$ and $t>\min\{t^{\prime},s^{\prime}\}$. Then for each $t_{i}$ there are two $p_{i}$ such that $y(u,p_{i})=y(v,p_{i})$.

**Lemma **\thelem

Suppose $g(t,\tau)$ and $g(t,\sigma)$ are the respective responses to the $\gamma$–optimal inputs $u(t)$ with duration $\tau$ and $v(t)$ with duration $\sigma$. Then $\gamma(u)>\gamma(v)$ if and only if $[s^{\prime},s^{\prime}+\sigma]\supset[t^{\prime},t^{\prime}+\tau]$, where $t^{\prime}$ and $s^{\prime}$ are the delivery times of the inputs $u(t)$ and $v(t)$, respectively.

Proof:

Assume that $\gamma(u)>\gamma(v)$ and suppose $s^{\prime}>t^{\prime}$. This implies that $t_{g,1}<t_{1,\max}$ as $g(u)<g(v)$ for all $t\in(t^{\prime},t_{g,1})$. Hence $t_{g,2}\leq t_{\min}$. This would imply either that $g(v)\geq g(u)=\gamma(u)$ at $t_{2,\max}$ or by Section III-B that there are two additional $p_{i}$ such that $y(v,p_{i})=y(u,p_{i})$. This contradicts Section III. Hence $s^{\prime}\leq t^{\prime}$. Now suppose either $s^{\prime}=t^{\prime}$ or $s^{\prime}+\sigma=t^{\prime}+\tau$. This implies that there is at most one intersection point, and that $\hat{v}<\hat{u}$, which implies that either $g(u)$ or $g(v)$ is non-optimal.

Suppose that $[s^{\prime},s^{\prime}+\sigma]\supset[t^{\prime},t^{\prime}+\tau]$. This implies that $g(u)>g(v)$ for all $t\in(t^{\prime},t_{g,1})$. If $t_{g,1}\leq t_{1,\max}$. Then $t_{g,2}<s_{\min}$. Should $s_{\min}<t_{\min}$ then $g(v)$ would be non-optimal. Thus $s_{\min}\geq t_{\min}$ after which $g(v)\leq g(u)$ which implies $g(v)$ is not $\gamma$–optimal as $g(v,s_{1,\max})\neq\max\{g(v,t):t\geq s_{\min}\}$. Therefore $t_{g,1}>t_{1,\max}$. This together with the assumption that $g(v)$ is $\gamma$–optimal implies that $\gamma(v)<\gamma(u)$. ∎

III-C Amalgamation

**Lemma **\thelem

Suppose $u$ is $\lambda$–optimal and $v$ is $\gamma$–optimal. Then $v$ is nested in $u$.

Proof:

For the sake of contradiction suppose $u$ is nested in $v$. We know that $t^{\prime}$ must occur before $t_{1,\min}$ which, by the assumption that $u$ is nested in $v$ implies that $s^{\prime}$ occurs before $t_{1,\min}$. Hence, there exists $t_{g,i}$ at or before each minimum of both $g(u)$ and $g(v)$. As there are at least three minima and at most two possible $t_{g,i}$ we see that $u$ cannot be nested in $v$. Instead, suppose, $u$ and $v$ are not nested. From Section III this implies that there is at most one intersection point of $g(u)$ and $g(v)$ contradicting optimality of $u$ or $v$. ∎

**Lemma **\thelem

Suppose $g(u)$ is $\lambda$–optimal and $g(v)$ is $\gamma$–optimal such that $\gamma(u)=\gamma(v)$. Then there exists an input $m$ with duration $\tau_{\lambda}>\sigma>\tau_{\gamma}$ such that $\gamma(m)<\gamma$.

Proof:

Choose $\sigma$ as in the statement of the Theorem. We know the input $m$ is either $\gamma$–optimal or $\lambda$–optimal. In either case as it satisfies the conditions of Lemmas III-A and III-B we have that $\gamma(m)<\gamma$. ∎

Theorem 3

Suppose there exists $\tau>0$ such that $g(v)$ is $\gamma$–optimal and $t$ such that $g(t)>g(0)$. Then an input $u$ of the form (3) is globally optimal if and only if the response has at least two global maxima interlaced between two minima.

Proof:

Let $\tau$ be the duration as in the statement of the Theorem. By Lemma III-B so long as $\sigma>\tau$ produces a $\gamma$–optimal input $v$ then $\gamma(v)<\gamma(u)$. As $w$ vanishes at infinity, there exists a duration $\alpha$ such that the input $u$ is $\lambda$–optimal.

Denote by $g(\alpha)$ the response of $g(t)$ to the input $u(t,\alpha)$ and by $g(\sigma)$ the response of $g(t)$ to the input $v(t,\sigma)$. We now construct a globally optimal $g$ and show that its shape is as in the statement of the Theorem. Recursively define the sequences $\overline{\alpha}:=(\alpha_{i})_{i=0}^{\infty}$ and $\overline{\sigma}:=(\sigma_{i})_{i=0}^{\infty}$ by $\alpha_{0}:=\alpha$ and $\sigma_{0}:=\sigma$ and $\alpha_{i}$ the least element of the following finite ordered partition of the interval $[\sigma_{i-1},\alpha_{i-1}]$:

[TABLE]

where $n\in\mathbb{N}$ is arbitary and $k_{i}\leq n$, such that the response:

[TABLE]

is $\lambda$–optimal. Similarly, $\sigma_{i}$ is defined to be the greatest element of $L_{i}$ such that:

[TABLE]

is $\gamma$–optimal. The sequence $\overline{\sigma}$ is a monotone increasing sequence bounded above by $\alpha_{i}$ for all $\alpha_{i}\in\overline{\alpha}$ and therefore has a limit $\tau_{-}$. Similarly, $\overline{\alpha}$ is a monotone decreasing sequence bounded below by $\sigma_{i}$ for all $\sigma_{i}\in\overline{\sigma}$ and thus has a limit $\tau_{+}$. It remains to show that $\tau_{-}=\tau_{+}$. Suppose, for all $i$, that $\sigma_{i}<\alpha_{i}$. We see that if:

[TABLE]

Then $\alpha_{i+1}$ must be the next element of $L_{i}$, as if were not the next element of $L_{i}$ would be $\gamma$–optimal contradicting our choice of $\sigma_{i+1}$, that is:

[TABLE]

Thus:

[TABLE]

i.e. $\lim_{i\to\infty}\left(\alpha_{i+1}-\sigma_{i+1}\right)=0$ i.e. $\tau_{-}=\tau_{+}$. Set $\tau:=\tau_{+}$.

Thus for all $\varepsilon\in(0,\tau)$ the optimal input with duration $\tau+\varepsilon$ must be $\lambda$–optimal an the optimal input with duration $\tau-\varepsilon$ must be $\gamma$–optimal. By continuity of $g$ there must be at least two equal maxima and two minima.

As the limits, $\tau_{-}$ and $\tau_{+}$, are equal and $g$ is a continuous function of the duration we may consider the sequence $(g(\alpha_{i}))_{i=0}^{\infty}$ to determine the shape of $g(\tau)$. Since the durations $\alpha_{i}$ decrease we have, as in the proof of Section III-A that $[t_{1,\min,i+1},t_{2,\min,i+1}]\subseteq[t_{1,\min,i},t_{2,\min,i}]$ for each $i$. Thus:

[TABLE]

where $B:=[t_{1,\min,i},t_{2,\min,i}]^{c}$ – the complement in $\mathbb{R}_{+}$. Indeed: $\gamma(u_{i+1})|_{B}\geq\gamma(u_{i})|_{B}$. Therefore the response $g(u_{i+1})$ has no additional minima outside the interval $[t_{1,\min,i},t_{2,\min,i}]$. By (5), the function $G(i):=\gamma(u_{i})-\gamma(u_{i})|_{B}$ is monotone decreasing, as the global maximum, $\gamma(u_{i})$, is decreasing by Section III-A, and bounded below by [math]. Suppose $\lim G(i)>0$. This is true only if $\tilde{u}:=\lim u_{i}$ is $\lambda$–optimal. If $g(\tilde{u})$ is $\lambda$–optimal then there exists strictly positive $\varepsilon<\gamma(\tilde{u})-\max\{g(\tilde{u}):g(\tilde{u})\neq\gamma(\tilde{u})\wedge\dot{g}(\tilde{u})=0\}$. For all such $\varepsilon$ there is $\delta>0$ such that $u(\tau-\delta)$ is $\lambda$–optimal as $g(\tau)$ is continuous. By Section III $\gamma(u(\tau-\delta))<\gamma(\tilde{u})$. This implies that $\tilde{u}$ is not the limit of the sequence of $\lambda$–optimal inputs $u_{i}$ with durations $\alpha_{i}$. Lastly, $\lim G(i)=0$ implies that $\gamma(\hat{u})=\gamma(\hat{u})|_{B}$ i.e. the maxima outside the interval $[t_{1,\min},t_{2,\min}]$ are equal to the maxima inside the interval.

Suppose $g(u)$ is as in the Theorem but there exists proper $v\neq u$ such that $\gamma(v)<\gamma(u)$. Thus, for $j\in\{1,2\}$, we require that $g(v)<g(u)$ at $t_{j,\max}$ and $g(v)\geq g(u)$ at $t_{j,\min}$. Therefore there must be more than two $t_{i}$ such that $h(u)=h(v)$, unless both $t_{g,i}$ occur at $t_{j,\min}$. In this case, by III-B, there must be at least four distinct points at which $y(u)=y(v)$ contradicting Section III. ∎

It is not possible to lower the global maximum of $g$ if $g(t)\leq g(0)$ for all $t$. However, in this case the shape of the reponse, $g$, specified in 3 does minimise the maximum of $h$. This is shown by Section III-C.

**Proposition **\theprop

Suppose $g(v)\leq g(0)$ for all $t$ and proper inputs $v$. Additionally, suppose, there exists $u$ for which there are two minima and there is a $t\in(t_{1,\min},t_{2,\min})$ such that the response $g(u(t))=g(0)$. Then $\max\{h(u)\}<\max\{h(v)\}$ for all proper $v\neq u$.

Proof:

The proof follows if all such $v$ are nested in $u$ as $h$ is a monotonic function of $\hat{u}$. As $g(u)$ and $g(v)$ may intersect at most twice we have that $t_{g,1}\in[t_{1,\min},t_{1,\max}]$ and $t_{g,2}\in[t_{1,\max},t_{2,\min}]$ and that $g(v)<g(u)$ for all $t\in(t_{g,1},t_{g,2})$. Thus $g(v)>g(u)$ for all $t>t^{\prime}$ such that $t\not\in[t_{g,1},t_{g,2}]$. This occurs only if $v$ is nested in $u$. ∎

III-D Optimal Duration for a Fixed Delivery Time

We conclude this Section by characterising the optimality of an input with varying duration and fixed delivery time. 4 is the analogous result for durations to the results for delivery times derived in [13]. Part 1 of 4 is a generalisation of the main result, (Theorem 16), of [14], when restricted to the class of rectangular inputs. This generalisation stems from only assuming that $w$ is bounded and vanishes at infinity and as we do not require that the global maximum, of $g$, occurs before its global minimum222This assumption holds if case A from 4 holds for $\tau=0$. As, if case A holds for $\tau=0$. Then it holds for all $\tau\geq 0$. for all inputs $u$, of the form (3). The results of [14] hold for more general inputs of the form $u(t)=\overline{u}+\hat{u}(t^{\prime},t)$ where $\hat{u}(t^{\prime},t)$ is a positive bounded function such that $\hat{u}(t^{\prime},t)=0$ for all $t<t^{\prime}$ and $\overline{u}$ is as in (3).

Theorem 4

Consider the following two cases for $g(t)$:

1. A.

no global maximum occurs after a global minimum. 2. B.

no global maximum occurs before a global minimum.

Fix $t^{\prime}$. Let $u(t^{\prime},\tau)$ and $v(t^{\prime},\sigma)$ be two distinct inputs delivered at $t^{\prime}$. Suppose either: $u$ and $v$ satisfy A, $u$ and $v$ satisfy B. or $u$ satisfies A and $v$ satisfies B. Then, for each respective case:

$\gamma(u)<\gamma(v)$* if and only if $\tau<\sigma$.* 2. 2.

$\gamma(u)<\gamma(v)$* if and only if $\tau>\sigma$.* 3. 3.

there exists $\alpha\in(\sigma,\tau)$ such that $\gamma(m(t^{\prime},\alpha))<\min\{\gamma(v),\gamma(u)\}$. Furthermore, for all $\alpha\not\in(\sigma,\tau)$ the maximum $\gamma(m(t^{\prime},\alpha))\geq\min\{\gamma(v),\gamma(u)\}$.

Proof:

Throughout this proof we say $g(t,\tau)<g(t,\sigma)$ initially if there exists $\varepsilon>0$ such that $g(t,\tau)<g(t,\sigma)$ for all $t\in(t^{\prime},t^{\prime}+\varepsilon)$.

Part 1

Suppose $\tau<\sigma$. Then, by 1 Part 2 and as $u$ and $v$ are proper, $\hat{u}>\hat{v}$. Thus, initially $g(t,\tau)<g(t,\sigma)$. By Section III there is at most one $t_{g}$ at which $g(t_{g},\tau)=g(t_{g},\sigma)$. As $u$ and $v$ are proper this $t_{g}$ must exist and $t_{g}\leq s_{1,\min}$. Otherwise $g(s_{1,\min},\tau)<\lambda$. Note that Section III implies that $g(u)-g(v)$ must change sign at $t_{g}$ even if $t_{g}=s_{1,\min}$. If $t_{g}\in(\max\{s_{\max}\},s_{1,\min}]$, where $\max\{s_{max}\}$ is the greatest time at which $g(v)$ is maximised. Then $\gamma(u)<\gamma(v)$ and $\{t:g(t,\tau)=\lambda\}\subset(\max\{s_{\max}\},s_{1,\min}]$ i.e. all minima of $g(u)$ occur between the last maximum of $g(v)$ and the first minimum of $g(v)$.

Instead, suppose $t_{g}\leq\max\{s_{\max}\}$. This implies that $\gamma(u)$ occurs after $t_{g}$. Thus, as $g(t,\tau)>\lambda$ for all $t\leq\max\{t_{\max}\}$ and $g(u)>g(v)$ for all $t>t_{g}$, $u$ is not proper.

Suppose $\gamma(u)<\gamma(v)$. If initially $g(t,\tau)<g(t,\sigma)$ then $\tau<\sigma$. Suppose, initially $g(t,\tau)>g(t,\sigma)$. As $g(v)>g(u)$ for all $t>t_{g}$ we have that $\max\{s_{\min}\}\leq t_{g}$. Otherwise $u$ would not be proper. By assumption $\max\{s_{\max}\}<\min\{s_{\min}\}$. Therefore $\max\{s_{\max}\}<t_{g}$. Then as there is at most one intersection point, of $g(u)$ and $g(v)$, we see that $t_{g}>\max\{t_{\max}\}$. Contradicting our assumption that $\gamma(u)<\gamma(v)$.

Part 2

Suppose $\tau>\sigma$. This implies that initially $g(t,\tau)>g(t,\sigma)$. Hence, similarly to Part 1 above, we see that $t_{g}\geq\max\{s_{\min}\}$. Thus all minima of $g(v)$ must occur before $t_{g}$ and all minima of $g(u)$ must occur after $t_{g}$. Thus $t_{g}\in[\max\{s_{\min}\},\min\{t_{\min}\}]$. In particular, this implies that $t_{g}<\min\{t_{\max}\}$. As $g(v)>g(u)$ for all $t>t_{g}$ we have, by the assumed shape of $g(u)$, that $g(t_{\max},\sigma)>\gamma(u)$ i.e. $\gamma(v)>\gamma(u)$.

Suppose $\gamma(u)<\gamma(v)$. If $t_{g}>\min\{t_{\min}\}$. Then $g(v)>g(u)\geq\lambda$ for all $t<t_{g}$. Otherwise there would exist $s$ such that $g(v)<\lambda$. Additionally, $g(v)<g(u)\leq\gamma(u)$ for all $t>t_{g}$. Thus all $s_{\max}<t_{g}$ and as no minimum of $g(v)$ exists after $\min\{s_{\max}\}$, we see that $v$ is not proper. Hence, $t_{g}\leq\min\{t_{\min}\}$. If $g(v)>g(u)$ initially. Then $\max\{g(v,s)\}<\max\{g(v,t)\}<\gamma(u)$ for $s\in[t^{\prime},t_{g}]$ and $t>t_{g}$. This is because no global maximum may occur before the last minimum of $g(v)$. Hence $g(v)<g(u)$ initially. This implies that $\sigma<\tau$.

Part 3

We note that the problem is well-posed as if $\sigma>\tau$ for $u$ satisfying A and $v$ satisfying B. Then initially $g(v)>g(u)$. Thus $t_{g}\leq\min\{s_{\min}\}$, after which $g(u)>g(v)\geq\lambda$. This implies either $u$ is not proper or $t_{g}=\min\{s_{\min}\}$. In which case there exists $s<t_{g}$ such that $g(v)>\gamma(u)$. As $v$ satisfies B there must exist $s>t_{g}$ such that $g(v)>\gamma(u)>g(u)$. Contradicting the uniqueness of $t_{g}$.

Suppose $\alpha<\sigma$. Then $m(t^{\prime},\alpha)$ satisfies B as $g(m)<g(v)$ initially. Thus, by the above, $\gamma(m)>\gamma(v)$. Similarly for $\tau$.

As $g$ is a continuous function of the duration there exists $\alpha<\tau$ such that $m(t^{\prime},\alpha)$ satisfies A and by the above $\gamma(m)<\gamma(u)$ and $\alpha>\sigma$. Similarly for $v$.

∎

Section III-D characterises when it is better to optimise the delivery time instead of the duration of an input.

**Corollary **\thecor

Fix $t^{\prime}$. There exists delivery time $s^{\prime}\neq t^{\prime}$ such that $\gamma(s^{\prime},\tau)<\gamma(t^{\prime},\tau)$ for all $\tau$ if and only if $u(t^{\prime},\tau)$ satisfies either A or B of 1 for all $\tau$.

We conclude this section with Section III-D which extends Sections III-A and III-B to the case of a fixed input time.

**Corollary **\thecor

Fix $t^{\prime}$. Suppose $u(t^{\prime},\tau)$ is proper and either $\lambda$ or $\gamma$–optimal. Then $\gamma(v(t^{\prime},\sigma))>\gamma(u(t^{\prime},\tau))$ for any $\sigma\neq\tau$ where $v(t^{\prime},\sigma)$ is proper. Futhermore, such $v$ is neither $\lambda$ nor $\gamma$–optimal.

Proof:

This follows by similar argument to the proofs of Sections III-A and III-B. ∎

IV Algorithm for Optimal Duration

As the duration, $\tau$, is bounded below by [math] the following algorithm may be used to locate the optimal duration:

Algorithm:

Set $\tau=0$ i.e. $\chi_{A}=\delta$. If the response $g$ is $\lambda$–optimal then $\tau$ is globally optimal. Otherwise: 2. 2.

Choose $\tau>0$:

1. (a)

if $g$ is $\lambda$–optimal then proceed to step 3 2. (b)

otherwise increase $\tau$ until $g$ is $\lambda$–optimal 3. 3.

Recursively bifurcate the interval $[\sigma,\alpha]$, where $\sigma$ is the largest known $\tau$ such that $g$ is $\gamma$–optimal and $\alpha$ is the least known $\tau$ such that $g$ is $\lambda$–optimal.

**Remark **\therem

If the condition that $g$ is $\lambda$–optimal is replaced by condition B from 4 and the condition that $g$ is $\gamma$–optimal is replaced by condition A from 4, this algorithm may be adapted to find the optimal duration for a fixed delivery time $t^{\prime}$.

Numerical Example:

In the example presented in Figures 1–3, the algorithm to locate the optimal duration was applied to a system where the parameters of (1) and (2) were chosen to be: $d=0.025$, $k=1806^{-1}$, $c=0.025$, $a=0.0101$, $b=8.16\times 10^{-4}$, $G=0.0023$, $E=1.0$, and $r(t)=263^{-1}f_{1}(t)$, where $f_{1}(t)$ is the solution to the system of linear differential equations:

[TABLE]

where: $\rho(t):=5\chi_{[300,800]}(t)+100\chi_{[450,460]}(t)$. We take the initial conditions to be as in Section II and set $g(\infty)=g(0)=5.0\mathrm{mmolL}^{-1}\,(90\mathrm{mgdl}^{-1})$. The minimum glucose concentration $\lambda$ is chosen to be $4.0\mathrm{mmolL}^{-1}\,(72\mathrm{mgdl}^{-1})$.

For computational reasons the smallest duration tested was $\sigma=2$. The longest duration considered was $\alpha=1000$. In Figures 1 and 3 the blue, green, dashed black, red and cyan lines correspond to the durations $\tau=100,250,370,550$ and $600$ respectively.

As the duration approaches $\tau=370$, which corresponds to the dashed black profile of Figure 1, the maximum glucose concentration decreases. This is shown in Figure 2 Indeed $\gamma(u)$ is monotonic as $\tau\to^{-}370$ and monotonically decreasing as $\tau\to^{+}370$. The small deviations are an artifact of the numerical precision. With no optimisation: $\gamma(u(100))=7.2\mathrm{mmolL}^{-1}\,(129\mathrm{mgdl}^{-1})$ and $\gamma(u(700))=8.3\mathrm{mmolL}^{-1}\,(149\mathrm{mgdl}^{-1})$. Whilst $\gamma(u(370))=5.54\mathrm{mmolL}^{-1}\,(100\mathrm{mgdl}^{-1})$.

Lastly, Figure 3 shows $\hat{u}\chi_{A}$ for the $A$ yielding glucose profiles shown in Figure 1. Each interval over which the input $u(t,\tau)\neq\overline{u}$ is nested in the next larger interval.

**Remark **\therem

Figure 2 indicates that the rate of decrease in $\gamma(u)$ drops about the optimal duration i.e. $\frac{d}{d\tau}\gamma(u,\tau)\to 0$. Therefore, it seems that there is little benefit in over-optimising the duration.

V Application to Other Models

For models of the form:

[TABLE]

where $g$ is the plasma glucose concentration, $u$ is the insulin input and $w$ and $h$ are some bounded positive functions, the results of [13] and those presented here require that:

$g$ is a continuous function of $u,w$ and $t$ that decreases monotonically with respect to $u$ 2. 2.

$h$ is continuous function of $u$ and $t$ and is monotone in $u$. 3. 3.

$w$ and $h$ decay to their respective lower bounds. 4. 4.

$g\geq g(\infty)\geq\lambda$ if $u(t)=0$, for all $t$, where $g(\infty)$ is the desired steady-state glucose concentration.

The Hovorka model ([15]) is another dynamic model of glucose metabolism which explicitly includes a number of physiological factors, for example a renal excretion term. We give numerical examples of our results for the Hovorka model which suggest it may satisfy our assumptions. The Hovorka model is:

[TABLE]

where:

[TABLE]

and:

[TABLE]

where $a_{i},b_{i},c,d,E,f,\overline{g}_{c},\overline{g}_{r},k,l,V$ and $R$ are positive constants, physiological values for which may be found in [15] and $r$ is a positive bounded function. The plasma glucose $g$ is a scalar mulitple of $q_{1}$:

[TABLE]

As for the Bergman model, we assume $u$ is a positive bounded function of the form (3) such that $\lim_{t\to\infty}u(t)=\overline{u}$. We also assume that $\lim_{t\to\infty}r(t)=0$. The steady-state value has a positive upper bound, $g>0$, for $\overline{u}=0$. As $g$ is a continuous function of $u$, we may choose $\lambda$ and the steady-state value, $g(\infty)=g(0)$, to be any value less than this upper bound.

V-A Numerical Examples

In the example presented in Figures 4–6 the algorithm to locate the optimal duration was applied to the Hovorka model with parameters as in [15] and a body weight, on which the values in [15] depend, of $70\mathrm{kg}$. The function $r(t):=\left(\frac{1}{55}\right)f_{1}(t)$ where $f_{1}(t)$ is the solution to the differential equations:

[TABLE]

where: $\rho(t):=0.2\chi_{[300,800]}(t)+5\chi_{[450,460]}(t)$. The initial conditions where set such that the system was in steady-state at $t=0$ for a steady-state value $g(\infty)=g(0)=5.0\mathrm{mmolL}^{-1}\,(90\mathrm{mgdl}^{-1})$. The minimum glucose concentration $\lambda$ was chosen to be $4.0\mathrm{mmolL}^{-1}\,(72\mathrm{mgdl}^{-1})$.

The duration tested ranged from $\sigma=2$ to $\alpha=1000$. In Figures 4 and 6 the blue, green, dashed black, red and cyan lines correspond to the durations $\tau=150,250,305,450$ and $550$ respectively.

As the duration approaches $\tau=305$, which corresponds to the dashed black profile of Figure 4, the maximum glucose concentration decreases monotonically from above and below. With no optimisation: $\gamma(u(2))=9.8\mathrm{mmolL}^{-1}\,(176\mathrm{mgdl}^{-1})$ and $\gamma(u(1000))=17.5\mathrm{mmolL}^{-1}\,(315\mathrm{mgdl}^{-1})$. Whilst, the optimal duration $\gamma(u(305))=7.7\mathrm{mmolL}^{-1}\,(139\mathrm{mgdl}^{-1})$.

Figure 6 shows $\hat{u}\chi_{A}$ for the $A$ yielding the glucose profiles shown in Figure 4. Each interval over which the input $u(t,\tau)\neq\overline{u}$ is nested in the next larger interval.

In the example presented in Figures 7–8, we demonstrate that the results of [13] apply to the Hovorka model. All values and functions were taken to be as in the previous example with a fixed duration $\tau=200$. Figure 7 shows three responses of the Hovorka model to a proper pulse delivered at $300,339$ and $337$, these correspond to the blue, green and red responses, respectively. The green response has two equal maxima bounding the minimum and for which $\gamma=8.7\mathrm{mmolL}^{-1}\,(157\mathrm{mgdl}^{-1})$. The dashed black line is the optimal glucose concentration achieved for this system in the previous example i.e. when both the input time and duration where optimised.

In Figure 8 the maximum plasma glucose and magnitude of the proper input bolus $\hat{u}$ is shown as a function of the input time $t^{\prime}$. The lowest maximum occurs at $t^{\prime}=339$, corresponding to the green response in Figure 7.

VI Conclusions and Further Work

We have given necessary and sufficient characterisations of the optimality of pulse inputs to the Bergman minimal and Hovorka models in terms of the shape of the predicted plasma glucose concentration. This paper, in conjunction with [13], determines the magnitude of the maximum glucose concentration in response to changes in the parameters of a pulse input. These results demonstrate the possibility of rejecting disturbances by tuning the duration and delivery time of a bolus input of some shape.

Current research aims to generalise the presented results to any bounded input function $u(t)$. We are also interested in characterising the behaviour of $\frac{d}{d\tau}\gamma(u(\tau))$ – the rate of change of the maximum of the response $g$ as a function of the duration $\tau$. This may provide conditions which guarantee the existence of $g(t)>g(0)$, for all durations, or a $\gamma$–optimal input, which are required for 3.

Given the general nature of the proofs of the current results we believe it is likely that similar results hold for other models of glucose metabolism.