Discrete Blood Glucose Control in Diabetic Göttingen Minipigs

Abstract: Despite continuous research effort, patients with type 1 diabetes mellitus (T1D) experience difficulties in daily adjustments of their blood glucose concentrations. New technological developments in the form of implanted intravenous infusion pumps and continuous blood glucose sensors might alleviate obstacles for the automatic adjustment of blood glucose concentration. These obstacles consist, for example, of large time-delays and insulin storage effects for the subcutaneous/interstitial route. Towards the goal of an artificial pancreas, we present a novel feedback controller approach that combines classical loop-shaping techniques with gain-scheduling and modern H∞-robust control approaches. A disturbance rejection design is proposed in discrete frequency domain based on the detailed model of the diabetic Göttingen minipig. The model is trimmed and linearised over a large operating range of blood glucose concentrations and insulin sensitivity values. Controller parameters are determined for each of these operating points. A discrete H∞ loop-shaping compensator is designed to increase robustness of the artificial pancreas against general coprime factor uncertainty. The gain scheduled controller uses subcutaneous insulin injection as a control input and determines the controller input error from intravenous blood glucose concentration measurements, where parameter scheduling is achieved by an estimator of the insulin sensitivity parameter. Thus, only one controller stabilises a family of animal models. The controller is validated in silico with a total number of five Göttingen Minipig models, which were previously obtained by experimental identification procedures. Its performance is compared with an experimentally tested switching PI controller.


Introduction
Diabetes mellitus collectively denotes a group of metabolic diseases with a worldwide increasing prevalence over the last decades. Diabetes is characterised by abnormally high blood glucose concentrations (hyperglycaemia) resulting from a dysfunction of insulin secretion, insulin action, or both [1]. In diabetes, hyperglycaemia is typically characterised by blood glucose concentrations beyond 126 mg/dL (7 mmol/L) in fasting state or beyond 200 mg/dL (11.1 mmol/L) in oral glucose tolerance tests and by increased HbA1c fraction (glycated haemoglobin >6.5%). According to a study on global estimates of diabetes prevalence [2], 382 million people in 130 countries had diabetes in 2013. This number is estimated to increase to 592 million people by 2035. A subgroup of diabetes mellitus is type 1 diabetes mellitus (T1D), which is characterised by the destruction of pancreatic β-cells (in healthy state about 65% to 80% of the islets of Langerhans) and a subsequent lack of the hormone insulin, which is important in the blood glucose regulation. About 10% of all diabetic patients are patients with T1D [3]. Diabetes might induce long term secondary diseases that may affect heart, vascular system, peripheral nervous system, and can lead to blindness, heart disease or limb amputation.
To this day, manual control has to be used for patients with T1D to reduce increased blood glucose levels to the normophysiological range. This is achieved by (a) a measurement of the blood glucose concentration; (b) the determination of an appropriate size of subcutaneously injected insulin bolus; and (c) the subcutaneous insulin injection. An overview of the feedback control system is presented in Figure 1. The patient can be seen as a feedback controller acting on its own metabolic dynamics. However, manual control of blood glucose bears several severe disadvantages. The glucose metabolic system is subject to external disturbance and internal dynamics, which affect the blood glucose level after insulin bolus application. Blood glucose predictions that are made by the patients rely heavily on expert knowledge and on sparsely available blood glucose measurements. In practice, the manual blood glucose control often leads to so-called hypo-or hyperglycaemic events. In both conditions, too low or too high blood glucose concentrations, respectively, can lead to acute problems due to too low or too high blood glucose levels.
A continuous research effort towards an automatic control of blood glucose, the so-called artificial pancreas (AP), was initiated in the early 1970s. The AP consists of three main components: (i) an insulin infusion pump; (ii) a blood glucose sensor; and (iii) a feedback control algorithm that computes insulin infusion rates based on glucose sensor measurements. Note that other additional measurements such as meal size or physical activity might be measured and integrated into the control algorithm. An overview of the history of the AP is, for example, given in [4]. First approaches relied on intravenous blood glucose measurements and insulin application, where the main challenge is the necessary catheterisation. More recently, the subcutaneous-to-subcutaneous AP have led to a tremendous increase in research. However, this approach is associated with multiple challenges, where it is commonly agreed that a large additional time-delay in the process dynamics is most prominent. Moreover, insulin storage in the subcutaneous compartment and sudden release is associated with a variable time-delay and the risk of severe hyperinsulinaemia by some researchers [4]. A first promising scenario for the application of AP might therefore be the intensive care unit with respect to central venous catheterised patients with blood glucose imbalances [5].
Over the last three decades, a number of AP approaches were proposed by different research groups. A recent review on the state-of-the-art T1D control strategies is presented in [6]. Early controllers used, for example, proportional-integral-derivative (PID) strategies. An example is a nonlinear PD-type controller presented by Clemens [7], which was extended by moving average filtering of the measured blood glucose. However, Hernjak and Doyle [8] concluded that a simple PD controller structure is not able to stabilise the blood glucose level around a tight operating point as the nonlinear, time-varying metabolic system is subject to external disturbances. More recent approaches to AP rather focus on modern optimal or robust control methodologies [9], model predictive control (MPC) [10] or the combination of feedback and feedforward control [11,12] (which is sometimes called two-degrees-of-freedom control).
An AP should reliably stabilise the blood glucose concentration while rejecting physical activity and meal uptake metabolic disturbances with satisfactory performance. The controlled system is nonlinear, has time-varying dynamics and shows intra-individual and inter-individual parameter uncertainties. In this article, we present a novel controller design procedure that is fully discrete. In contrast to [9], our controller design procedure combines the advantages of classical frequency domain techniques with fixed structure controllers with modern internal model-based procedures for robust control, for example H 2 -or H ∞ -design procedures. The central controller design is basically a discrete frequency domain disturbance rejection procedure. The design is conducted at a large number of operating points of the linearised process models. Gain scheduling of the discrete controllers is proposed on the measured blood glucose concentration and the estimated insulin sensitivity, thus guaranteeing an adaption to changing process conditions. The classical controllers are finally robustified by the discrete version of the Glover-McFarlane H ∞ -loop-shaping procedure [13].
The AP approach in this paper is model-based and tested in an extensive in silico validation procedure with n = 5 Göttingen Minipig models that were obtained from previously conducted animal experimental studies: • In Section 2, we briefly describe the mathematical model of a Göttingen Minipig, which is based on the Sorensen model [14] and has already been published [15,16]. • The controller design procedure is based on the Göttingen Minipig model and is introduced in Section 3. • Section 4 describes the animal experimental study and the results of the model identification procedure that is the basis for the controller validation. • In Section 5, the controller performance is evaluated in silico by means of the identified animal models and compared to the performance of a structural switching P/PI-controller that was validated in vivo [17]. • Discussion of the results and a summary are given in Section 6.

Göttingen Minipig Model
The metabolic model of a Göttingen Minipig was adopted from the Sorensen model [14,16]. An overview of the model is presented in Figure 2. The model consists of three main subsystems: the blood circulation, the interstitium and the gastro-intestinal tract. The models of the interstitium and the gasto-intestinal tract are interconnected via the subcutaneous insulin appearance rate r ISC (t) and the intravenous glucose appearance rate r GGA in a forward manner, respectively. External inputs to the model can be grouped into disturbances and manipulated variables. Disturbances that enter the model are the pancreatic glucagon infusion rate S Γ (t) = r N PΓP (t)r B PΓP = r PΓP (t) and orally and intravenously applied glucose rates, denoted by D oral (t) and D iv (t), respectively. The glucagon production r N PΓP (t) is normalised to the basal rate by r B PΓP . Manipulated variables are the subcutaneous U sc (t) or the intravenous U iv (t) insulin infusion rates. G p (t) denotes the blood glucose concentration that is measured in the plasma.

Nonlinear State Space Model
The Göttingen Minipig model consists of n x = 16 first-order ordinary differential equations. The control and disturbance inputs are subsequently grouped as u(t) = [U sc (t), U iv (t)] T ∈ R 2 and d(t) = [D oral (t), D iv (t), D Γ (t)] T ∈ R 3 , with external glucagon infusion rate D Γ (t). With the state vector x(t) ∈ R n x and the blood glucose concentration y(t) = G P (t) ∈ R as system output, the nonlinear state space model is given bẏ where f(·) : R n x × R + → R n x is a nonlinear smooth vector field. c, B and E are linear mappings with according dimensions. As further simplifications, we assume that only subcutaneous insulin is available as a control input, that is U iv (t) = 0. Moreover, only oral glucose disturbance D oral (t) is regarded in the controller design; therefore, D iv (t) = 0 and D Γ (t) = 0.

Model Structure
The equations of the model are derived on the basis of a compartmental modelling approach. Figure 3 gives an overview depicting an interaction diagram of a general compartment i of the model. Figure 3. Basic compartment i of the blood glucose model illustrating mass exchange of substance X.
The mass exchange takes place between a vascular (blood) compartment V and an interstitial (bloodless) compartment I. The mass exchange in compartment i can be described by two coupled first-order differential equations where X(t) denotes the mass distribution of a particular substance in the corresponding space. In Equation (2), Q X i denotes the blood stream, T X i is a time constant associated with mass diffusion and additional mass inflow and outflow are denoted by r in (t) and r out (t), respectively. The compartment i is coupled to other compartments by the mass concentration in the inflowing blood stream X in (t) and the mass concentration in the outflowing blood stream X iV (t). The blood circulation equations of the Göttingen Minipig model and corresponding parameters are not provided in this paper. Instead, the reader is referred to [15].
The dynamics of orally uptaken glucose D oral (t) are described by a second-order model as also presented by Hovorka et al. [18]. By denotingṀ G ing,1 andṀ G ing,2 the solid and the liquid glucose mass flow in the stomach and the intestine, respectively, the equations are given by In Equation (3), T ing,1 and T ing,2 are time constants and f G is the bioavailability of glucose in the blood. Similar to the model of the gastro-intestinal tract, the interstitium model consists of two coupled linear differential equations of the first-order. The mass flows of nonmonomeric (inactive) insulin and monomeric (active) insulin are denoted byṀ I sc,1 (t) andṀ I sc,2 (t), respectively. Equations for these mass flows, that have been applied subcutaneously, are given by where T sc,1 and T sc,1 are suitable time constants and x U is a fractional parameter that is used to describe the partially activated insulin. Note that model nonlinearities are mainly due to hyperbolic functions and bilinearities with respect to certain states [14]. In the case of the time-varying parameter insulin sensitivity k IS (t), the glucose concentration in the interstitium G MI (t) is influenced by the glucose consumption rate r MGC (t), which is given by where V I MI denotes the interstitial space for insulin and I MI (t) denotes the interstitial insulin concentration. The mass balance for G MI (t) is given by with V G MI the interstitial balance volume for glucose, T G M the diffusion time constant and G MV (t) the glucose concentration in the muscle and adipose tissue. A change in k IS (t) affects the rate of glucose consumption and increases or decreases the process gain of the blood glucose controlled plant (from U sc (t) to G p (t)).

Controller Design
The controller design presented in this section is based on a model linearisation, a model reduction and a discretisation. The novel design procedure consists of a classical loop-shaping design and a subsequent H ∞ -loop shaping robustification. To form a gain scheduling controller by lookup table interpolation, the classical loop-shaping design is presented as an automatic procedure over a large number of operating points.

Model Linearisation
The nonlinear state space model given in Equation (1) is trimmed and linearised with a first-order Taylor approximation. External disturbances and the full input vector are used in the trimming procedure. However, corresponding entries are neglected in the linearisation procedure. The trimmed dynamics result in where the linearised state space model with state vector ∆x, no external disturbances d = 0 and an input consisting only of subcutaneous insulin ∆u is given by with the initial state vector ∆x 0 and (A eq , b, c) given by the corresponding Jacobians. The model is trimmed over a range of blood glucose operating points G p and insulin sensitivity values k IS . The set of blood glucose operating points is given by X and Y denote the integer sets of blood glucose and insulin sensitivity operating points, respectively. A total of 256 models are linearised for n G = 16 and n k = 16, which are collected in set of plants From the trimming procedure, the disturbance dynamics are obtained by linearisation around the same operating points. The linearised state space model with state ∆x d is given by which are collected in the set of disturbance dynamics

Discrete Controller Design Prerequisites
In the first step, the set of plants P and P d are discretised at a sampling time of T s = 15 min. A block diagram showing the loop structure is given in Figure 4.
Note that T s is chosen in accordance to the animal experimental study presented in Section 4. With respect to the linearised state space models given in minimum realisation the discrete realisation of the plant is given by where H 0 = H denotes a zero-order hold element ( Figure 4) and Z {·} is the z-transform Z {·} of the δ-sampled impulse response series with the inverse Laplace-transform L −1 {·}. The disturbance transfer function is discretised by applying the z-transform to the impulse response of the disturbance dynamics which corresponds to an impulse-invariant discretisation. After obtaining the discretised state space realisations G(z) and G d (z), the orders are subsequently reduced. Towards this end, G(z) and G d (z) are brought to a balanced realisation by employing the state transformationx = Tx [19], with respect to the discete-time controllability W c and observability W o gramians where T is chosen such thatW In Equation (17), Λ denotes a diagonal matrix filled with (positive) square roots of the eigenvalues Consequently, balanced realised systemsḠ(z) andḠ d (z) are reduced toĜ r (z) andĜ r d (z), respectively, by truncation of all states that do not correspond to the maximum eigenvalue in Λ. The resulting discrete transfer functions are given bŷ with positive constants b 0 , a 0 , b d,0 and a d,0 . Figure 5 shows a comparison of the linearised and discretised full order systemsḠ(z) andḠ d (z) and the reduced order systemsĜ r (z andĜ r d (z). The systems correspond to a blood glucose concentration of G p (t) = 50 mg/dL and an insulin sensitivity of k IS (t) = 0.4 mg/(min mU).  Full order G(z)

Disturbance Rejection Design
This section presents the classical loop-shaping controller for disturbance rejection design. The presented procedure consists of three design stages, which are automatically computed over a large operating range of linearised and discretised models. A series interconnection of controllers C 1 (z), C 2 (z) and C 3 (z), as a result of the three design stages, forms the classical loop-shaping controller.

Initial Controller Gain for Disturbance Rejection
The controller is designed for disturbance rejection. Therefore, the reference input is set to zero, r(kT s ) = 0 (see Figure 4). Remaining control loop inputs and outputs are scaled with respect to their maximum possible (expected) magnitudes [20]. Scaling factors are introduced for the error d e = e max , the process input d u = u max and the disturbance input d d = d max . The scaled and discretised plant of reduced order G r d (z) and the scaled and discretised disturbance process of reduced order G r d (z) are then given by With respect to the discrete time disturbance signal d(z) and the error of the control loop e(z) = y(z), the objective of the controller is e(e jωT s ) < 1, assuming d(e jωT s ) with ω s = 2π/T s . Assuming a worst case disturbance of |d(e jωT s )| = 1 for some frequency ω, one can obtain S(e jωT s )G r d (e jωT s ) < 1, ∀ ω ∈ 0, which can be rewritten by introducing the sensitivity transfer function as S(z) = 1/(1 + L(z)), with L(z) = G r (z)C(z). The requirement for disturbance rejection thus follows as For frequencies ω with |G r d (e jωT s ))| > 1, it follows that |L(e jωT s ))| |G r d (e jωT s ))|. The initial controller should, therefore, satisfy the following gain requirement It follows from Equation (19) that the initial controller is given by where the superscript s denotes the parameters of the scaled transfer functions.

Integral Action
In the second design phase, the initial controller is extended with integral action to reduce the remaining steady-state error of the disturbance rejection. A PI-controller C(s) = K p + K i /s with proportional gain of K p = 1 is discretised by using the Tustin approximation. With s ≈ 2 T s z−1 z+1 , the discrete PI-controller is given by The value of K i is determined by employing the phase margin. Candidate integral gains K i n are obtained in the log space by the set with the initial, possibly negative exponent n 0 ∈ Z, the maximum exponent n 1 ∈ Z and the number N ∈ N + of candidate integral gains to be evaluated. The phase margin PM n is then determined over the set of candidate gains K for the full order open-loop compensated plant PM n = ∠ G(e jω PM T s )C 1 (e jω PM T s )C 2 (e jω PM T s ) where the largest gain K i n with a minimum phase margin of PM n > 30 • was chosen. The resulting controller has a relatively large bandwidth but an insufficiently low gain margin.

Lead-Lag Compensator
The third and final step of the classical loop-shaping procedure therefore consists of the design of a lead-lag compensator to increase the phase margin. A Tustin approximated lead-lag discrete compensator is used that has the form Time constants τ LL and τ LG correspond to the lead and the lag pole, respectively. By denoting the frequency of the previously determined phase margin of Equation (28) ω PM , the time-constant of the lead part is set τ LL = 1/(0.7ω PM ). Furthermore, the time-constant of the lag part is set to τ LG = 2 * τ LL . The maximum phase increase occurs at ω PM with a positive phase of ≈ 19.4 • . Note that the phase margin is increased by this value. However, the GM is reduced by the additional increase in gain. The classical loop-shaping procedure is repeatedly executed over the whole set of plants P and P d . Hence, the total number of 256 controllers for each C 1 (z), C 2 (z) and C 3 (z) are determined automatically over the whole range of sampled operating points.

Robust Loop-Shaping H ∞ -Controller
The robust loop-shaping H ∞ -controller follows the approach of Walker [13] and increases the robustness of the shaped loop with regard to discrete, not further specified, right coprime factor uncertainty. Given the discrete compensated plant in a normalised right coprime factorisation a family of coprime factor uncertain plants can be described by where [∆ M ∆ N ] provide a fairly general type of dynamic uncertainty, which allows both poles and zeros to cross out of the unit circle in the z-domain. For the generalised plant P(z) given in the form of a normalised left coprime factorisation the application of the small gain theorem (see for example [21]) ensures robust stability with respect to the normalised coprime factor perturbed family of plants P . Here, the optimal achieved performance by the H ∞ -norm γ = 1/ provides a measure of tolerable right coprime factor with respect to Equation (31). An internally stabilising controller is obtained given a minimal discrete state space realisation of P(z) and the solution of an indefinite discrete time algebraic Riccati equation. The central loop-shaping H ∞ -controller K(z) is employed in the overall mixed classical-robust controller This corresponds to a reduction of peak sensitivity measure M S = S(z) ∞ , with S(z) = 1/(1 + L(z)). The peak sensitivity measure is reduced from M S | L(z) = 4.7 dB to M S | L ∞ (z) = 2.28 dB. An = 0.5235 was accomplished by the the robustification procedure.
Note that only one robust H ∞ -compensator is designed. For the design, a model at an operating point corresponding to a blood glucose concentration of G p (t) = 130 mg/dL and an insulin sensitivity of k IS (t) = 2.0 mg/(min mU) is chosen. Stability is determined by analysing classical gain and phase margin, as well as sensitivity peak measure M S at every operating point of P. The resulting minimum gain and phase margin are computed as GM| L ∞ (z),min = 14.2 dB and PM| L ∞ (z),min = 63.6 • . The maximum sensitivity peak was determined as M S | L ∞ (z) = 3.1 dB. As a consequence, it is concluded that only one H ∞ -compensator is enough to robustly stabilise the set of plants P. Note that an additional sensor time-delay of about 2 min was neglected in the control design prodecure, as it was considered to be sufficiently small in comparison to the sampling time of 15 min. The additional time-delay is due to the time it takes from withdrawing venous blood until a sensor value is obtained. At the crossover frequency of the compensated open-loop of 0.02 rad/min, the time-delay leads to an additional phase lag of −2.3 deg.

Experimental Animal Study
To identify a realistic process and for the first test of a switching PI controller, an experimental animal study was conducted (the experimental study was approved by the State Agency for Nature, Environment and Consumer Protection, North Rhine-Westphalia, Germany) with Göttingen minipigs. A total number of eight animals were included in the study that were aged 2-3 years and weighted 40-70 kg upon arrival [17]. The animals were acclimatised to standardised conditions of the animal housing in the first three weeks. Each animal was equipped with two central venous lines (Cavafix Certo 475, B.Braun Melsungen AG, Melsungen, Germany), which were implanted into the right jugular vein. An acute diabetes was induced according to the protocol of Strauss et al. [22]. During this procedure, the β-cell toxin streptozotocin (#S0130, Sigma-Aldrich, St. Louis, MO, USA or #572201 Streptozotocin, Merck Millipore, Darmstadt, Germany) was applied to the animals followed by an intensive care phase of 36 h. After 10 additional days of convalescence, glucose and insulin tolerance tests were conducted. During glucose control tests, blood glucose was measured intravenously with a minimum time interval of T s . For this, 1 mL blood with no pre-treatment was withdrawn from the venous line and applied to a blood gas analyser (ABL 810 Flex, Radiometer Medical ApS, Herlev, Denmark). Short acting insulin (Humalog 100 U/mL, Lilly, Bad Homburg, Germany) was applied subcutaneously with a commercially available insulin pump (AccuCheck Combo Spirit, Roche Diagnostics, Rotkreutz, Switzerland). After euthanization, the pancreas was removed from the animals and histologically analysed by means of standard haematoxylin-eosin immunohistochemistry. An anti-insulin antibody (#IR002, Polyclonal Guinea Pig Anti-Insulin, Dako, Germany) was applied showing that only a small number of β-cells survived (5-50 cells per standardised area), compared to the reference pig (600-2800 cells per standardised area).

Model Identification
For each of the five minipigs, one parameter set could be identified, resulting in five individualised mathematical minipig models. Details on the model adaption are given in [16]. Parameters of the minipig model are partially transferred from the Sorensen model, in cases where human and animal metabolism is similar. However, a remaining set of 14 parameters was determined in a mathematical identification procedure. These parameters were determined from the minimisation of the cost function J where the error e(p, t) = G P,sim (p, t) − G P,meas (t) is computed from a comparision between in silico G P,sim (p, t) and in vivo G P,sim (p, t) blood glucose values. p ∈ R 14 denotes the parameter vector.
To reduce the parameter space, the parameter vector is divided into five subgroups p = [p T 1 , . . . , p T 5 ] T , the values of which were determined according to Equation (35). The corresponding experimental tests included intravenous and oral glucose applications, intravenous and subcutaneous insulin injections, and glucose concentration measurements. The identification process was implemented in MATLAB using the Levenberg-Marquardt optimisation algorithm (lsqnonlin, Optimisation Toolbox, The MathWorks, Natick, MA, USA). Figure 7 shows a comparison of in vivo data with results of the model identification procedure. In the experiment, intravenous glucose was administered at t = 0 min followed by two consecutive intravenous insulin impulses at t = 210 min and t = 334 min. In silico results evidently capture the real-life dynamics with high accuracy with respect to ±10% relative sensor error (red area in the figure). Similar results were obtained with other minipig models [15].

In Silico Feedback Control Study
In this section, the robust disturbance rejection control is implemented as a gain scheduled controller and tested in silico with the five minipig models. The controller performance is compared against the performance of a classical switched PI-controller that has already been tested in vivo with the diabetic Göttingen Minipigs [17].

Controller Implementation
The blood glucose controller is implemented as a gain scheduling controller. For each of the scheduled controller parameters, a lookup-table procedure with linear interpolation and clipping at table endpoints is used. The scheduling variables are computed by a state observer that is based on a nonlinear reduced-order version of the Göttingen Minipig model. The state observer model is subsequently augmented with the insulin sensitivity as a state. For this, an Extended Kalman Filter (EKF) is used that is based on the augmented nonlinear minimal realisation of the Göttingen Minipig model [23]. The gain scheduled controller is presented in Figure 8. Based on the control input signal and the measured blood glucose concentration, discrete time updates of the estimated parameter vector θ(kT s ) = Ĝ p (kT s )k IS (kT s ) T are computed for scheduling. It should be mentioned that the EKF does not include the disturbance dynamics of a meal uptake. Instead, the EKF is driven by forward dynamics of the gastro-intestinal tract, which is not shown in Figure 8. To simulate blood-gas sensor behaviour, normally distributed (Gaussian) noise was added at the blood glucose concentration output, given by N (0, 0.01). Figure 8. Block diagram of the discrete disturbance rejection controller C(z, θ) implementation with gain scheduling depending on estimated process parameters insulin sensitivityk IS (kT s ) and blood glucose concentrationĜ p (kT s ).

Disturbance Rejection Performance
To validate the performance of the controller, disturbance rejection tests are conducted with the full-order nonlinear minipig model and the reduced-order EKF estimator. The model is initialised with glucose concentrations of 327 mg/dL in all glucose compartments, a glucagon concentration and W HGP,0 = 2.237 [-] (the reader is referred to [15] for more details). All other states of the model are initialised to zero. An example of a meal disturbance rejection of the gain scheduled controller is shown in Figure 9. The controller and the EKF are switched on simultaneously at time t = 0 h and receive intravenous blood glucose samples every 15 min. After a transient response of ∆t = 1 h, the EKF reliably estimates the insulin sensitivityk IS (kT s ).
The controller is able to reduce the high initial blood glucose level within less than two hours with an undershoot to absolute glucose concentration of 77 mg/dL. In addition, the controller is able to keep the blood glucose below 150 mg/dL and above 95 mg/dL after a oral glucose uptake of D oral (t) = 50 g at t = 15 h.
The disturbance rejection performance is furthermore compared with the PI-controller that has been successfully implemented under in vivo test conditions in the Göttingen Minipig study. Figure 10 shows the comparison of the disturbance rejection performance to a meal consisting of 50 g oral glucose update. Orally uptaken glucose is applied to the model over the course of 1 min simulation time. This disturbance is given to the closed-loop control system with previously converged dynamics. A lower constant insulin sensitivity of k IS ≈ 1.175 mg/(min mU) is set-up in the model. Note that this PI-controller had already been tested successfully in experiments with diabetic Göttingen Minipigs. From Figure 10, it can be seen that both controllers provide good disturbance rejection. However, the system response to the PI-controller action is more oscillatory.
To investigate the effects of changing process dynamics, the insulin sensitivity k IS (t) is changed in a sinusoidal manner in the full order model. Figure 11 presents the results of the in silico test with the tested PI-controller and the new H ∞ -controller. The disturbance rejection response of the classical PI-controller is oscillatory, as the controller is not adapted to the current process dynamics. On the contrary, the gain scheduled H ∞ -controller is adapted online to the process. Since the increase in insulin sensitivity leads to an increase in process gain, the resulting controller, obtained from the automatic tuning procedure presented in Section 3, has a reduced gain. Moreover, the PI-controller shows severe undershoot after disturbance rejection with minimum blood glucose G p (t) = 59 mg/dL. On the contrary, the undershoot of the gain scheduled H ∞ -controller is bounded by a minimum value of G p (t) = 80.25 mg/dL. The disturbance rejection performance is not oscillatory. The EKF is able to track the insulin sensitivity over a large range of insulin sensitivity k IS (t) = [1.35 . . . 3.4] mg min·mU . The gain scheduled H ∞ -controller was tested for all five animal models that were obtained from the model identification procedure of Section 4.2. Only one gain scheduling H ∞ -controller was used. The controller provided stable results with sufficient performance in disturbance rejection, while time-varying effects of changing insulin sensitivity were included in the process model. A comparison with the classical PI-controllers, which were obtained for each of the Göttingen Minipigs individually, underlined the advantages of the proposed strategy. In some of the in silico tests, the classical PI-controller showed strong oscillations while rejecting meal disturbances containing cases with dangerously low blood glucose concentrations of below 60 mg/dL. Table 1 shows an overview of the percentage of time in hyper-or hypoglycaemic conditions of the robust gain-scheduled controller over a 36 h in silico trial with 50 mg carbohydrates meal uptake disturbance rejection and a highly time-varying insulin sensitivity. Only one gain-scheduled controller was employed for the family of in vivo identified plant models. The insulin sensitivity was varied in a sinusoidal manner with a period of T IS = 24 h and a variation in amplitude of ±30%.
The percentage of time in hyperglycaemia or hypoglycaemia was determined by employing bounds on the blood glucose concentration of 120 mg/dL and 70 mg/dL, respectively. For all measurements, a reference value of 100 mg/dL blood glucose concentration was used. For only one of the Göttingen Minipig models, an episode of hypoglycaemia was obtained, where the minimum blood glucose concentation was kept above the critical threshold of 50 mg/dL. 6.1% 0% #08 6.0% 1.3% A disturbance rejection statistics is provided in Figure 12. A single robust gain-scheduled controller was used to stabilise the family of plants, under the assumption of fixed insulin sensitivity. Disturbance in the experiments were 50 mg of glucose that were uptaken orally. Given in the figure are the mean and the standard deviation that were calculated over time for all conducted experiments. It can be seen that the controller quite well adapts to the models with individual parameter variability, while rejecting the disturbance sufficiently fast. Furthermore, it can be seen that the blood glucose response does not fall below 80 mg/dL.

Conclusions and Discussion
A novel discrete blood glucose control design procedure was developed that is focused on disturbance rejection of meal uptake. The procedure combines the advantages of classical loop-shaping control design with modern robust H ∞ control design. Based on the linearisation of the detailed and accurate minipig model, the controller design is automatically conducted for a large operating range of trimmed models at varying blood glucose concentrations and insulin sensitivities. The family of gain scheduled controllers was then robustified with a single H ∞ discrete loop-shaping controller. Classical gain and phase margin, as well as H ∞ -norm peak sensitivity measures computed for the family of linear models, indicate strong stability properties.
The gain scheduled controller was tested in silico by using the five different full order nonlinear minipig models. Controller implementation was realised with the help of an augmented EKF that estimates insulin sensitivity as a system state. The controller was able to robustly stabilise the family of animal experimental models where additional uncertainty in terms of time-varying insulin sensitivity was included. A comparison with a previously designed classical PI-controller showed superior performance of the new control strategy.
The controller in its current form was implemented without anti-windup protection of the discrete integrator. Moreover, the performance of the gain scheduled controller depends on the accuracy of the insulin sensitivity estimator. As meal uptake and insulin sensitivity cannot be estimated at the same time, due to observability problems, the size of a meal uptake is considered in the estimator by regarding the forward dynamics. Since the subcutaneous application of insulin is considered as uncertain by different authors, a first suitable application of the proposed discrete blood glucose controller seems to be an intensive care unit glucose control scenario. Here, the intravenous-to-intravenous route minimises control input uncertainties, time-delay, and blood glucose sampling time. In its current form, the presented controller design procedure can be readily applied to this scenario, if a suitable model of the process dynamics is available, which can be linearised over the range of blood glucose operating points and insulin sensitivity values. Furthermore, the controller extension with a forward term regarding meal uptake disturbance compensation like, for example, in [12] is possible. The application of the controller in a clinical scenario will be the subject of future work.
Author Contributions: Berno J.E. Misgeld was responsible for the controller and state observer design and the evaluation of the controller performance. Philipp G. Tenbrock was involved in the state observer design and the in silico validation. Katrin Lunze conducted the experimental study with the Göttingen Minipgs and developed the Göttingen Minipig models based on the Sorensen model. Steffen Leonhardt provided technical support and guidance throughout the project and carefully edited the paper.

Conflicts of Interest:
The authors declare no conflict of interest.