Investigation of stability in a two-delay model of the ultradian oscillations in glucose-insulin regulation

Abstract In this paper, a two-delay model for the ultradian oscillatory behaviour of the glucose–insulin regulation system is studied. Hill functions are introduced to model nonlinear physiological interactions within this system and ranges on parameters reproducing biological oscillations are determined on the basis of analytical and numerical considerations. Local and global stability are investigated and delay-dependent conditions are obtained through the construction of Lyapunov–Krasovskii functionals. The effect of Hill parameters on these conditions, as well as the boundary of the stability region in the delay domain, are established for the first time. Numerical simulations demonstrate that the model with Hill functions represents well the oscillatory behaviour of the system with the advantage of incorporating new meaningful parameters. The influence of the time delays on the period of oscillations and the sensitivity of the latter to model parameters, in particular glucose infusion, are investigated. The model can contribute to the better understanding and treatment of diabetes.


Introduction
Diabetes mellitus is a prominent problem worldwide with cases rising year on year. A 2011 Lancet report gave figures that there are 347 million diabetic sufferers worldwide [1]. Further to this, the World Health Organisation predicts that diabetes will be the 7th highest leading cause of death in the world by the year 2030 [2]. The disease itself is a defect of the metabolic system causing high blood glucose levels and has two main forms. Type I diabetes is a genetic disease which cannot be prevented where the pancreatic b-cells do not produce the regulatory hormone insulin. Type II diabetes is developed later in life and is much more common. It is caused by a large number of risk factors causing the glucose-insulin regulatory system to stop functioning correctly. The glucose-insulin regulatory system is an important part of the metabolic system. The pancreas and the liver regulate the production of insulin and glucose respectively, in order to keep the blood glucose at a normal level. Insulin is produced by the pancreas which is needed by the cells in the body to allow glucose to be used as fuel. Failure to do this will result in high blood sugar and hence diabetes which is related to many other long term problems. The system operates with ultradian oscillations which were first discovered by Hansen in 1923 [3] and have been studied in a variety of ways since.
The work in this paper to model the ultradian oscillations of glucose and insulin stems from the work of Sturis et al. [4][5][6]. The model [4] is based on a negative feedback loop, typically found in most biological systems, and incorporates two insulin compartments, plasma and remote, and one glucose compartment. An important feature of the system is the inclusion of a time delay for the effect of insulin on glucose production. The model is created as a system of six ordinary differential equations including a three tiered system representing the time delay. Our model includes five functions, f 1 À f 5 , which represent key features of the system, as illustrated by the modelling framework in Fig. 1.
The functions f i are commonly used in most models related to this problem and the model [4] has been used as a base for many models since its publication up to the recent paper [7]. The model [8] investigates the advantages of periodic versus constant insulin supply through a simplification of the original model [4] by replacing the functions f i with constants and Taylor expansions.
One main development from the model [4] is to incorporate explicitly the time delay in glucose production to function f 5 , creating a system of three delay differential equations as can be seen in [9]. The model is further simplified to a system of two delay differential equations by combining the two insulin compartments [10,11]. These models also explicitly introduce the second time delay, to represent the delay in insulin production. However, the analytical analysis of the two-delay models is not completely understood.
The analysis of the two-delay model focuses on the apparition of oscillations and the effect of varying certain biological parameters numerically. The work of Engelborghs et al. [10] carries out a numerical analysis of the steady state for their two equation model with two delays by finding branches of periodic solutions via a Hopf bifurcation. The paper of Li et al. [11] uses numerical bifurcation analysis for the two time delays and two selected parameters to give a threshold values leading to ultradian oscillations. A follow-up paper [12] provides an analytical study of local stability in the model [11] for both the one-and two-delay case. The global stability of a similar model was considered by Bennett and Gourley [13] in the one-delay case using the comparison principle approach and a Lyapunov functional.
Our work is based on the two-delay model of Li et al. [11]. One contribution of the current paper is the usage of Hill functions for each of the functions f 1 ; f 2 ; f 4 and f 5 which allows for more adequate modelling of underlying physiological dynamics. The Hill function was first introduced in 1910 as a mechanism for describing the interaction between oxygen and haemoglobin within the blood [14]. Since its creation it has been widely used in clinical and pharmacological models to great effect, specifically for modelling nonlinear relations. For reviews see for example [15]. Previously, the Hill function has been applied to the glucose-insulin regulation system but only as a single function within a simplified model [16][17][18]. One such example comes from the work of Huang et al. [17] which models an artificial pancreas. They consider a delay-free simplified version of the model [11] with one Hill function where the remaining functions are linear. This shows that these models can be used to aid the treatment of diabetes mellitus, which is one of the motivations of this work.
Furthermore, this paper for the first time provides a stability analysis both locally and globally of the two-delay system of glucose and insulin regulation. The linearisation of the system is studied analytically and a linear matrix inequality (LMI) approach is employed to numerically simulate the stability region. A Lyapunov functional approach is used to analytically give conditions for global stability, while providing numerical simulations to show the relation to the required physiological range.
This paper is organised as follows: the model with new functions f i is presented with a general analysis of parameter effect in Section 2. An analytical local stability analysis of the two-delay model is presented in Section 3 along with a numerical study based on LMI simulations. Section 4 provides the analysis of global stability for the two-delay case including an investigation of optimal ranges for the Hill parameters. Finally, concluding remarks are made in Section 5.

Mathematical model and Hill functions
Ultradian oscillations in the glucose-insulin regulatory system are studied by considering the following two-delay model _ GðtÞ ¼ G in À f 2 ðGðtÞÞ À f 3 ðGðtÞÞf 4 ðIðtÞÞ þ f 5 ðIðt À s 2 ÞÞ; _ IðtÞ ¼ f 1 ðGðt À s 1 ÞÞ À d i IðtÞ; ð2:1Þ as introduced by Li, Kuang and Mason in [11]. The functions GðtÞ and IðtÞ represent the plasma glucose (in mg=dl) and insulin (in mU=l) respectively while the two time delays s 1 and s 2 are assumed to be positive and constant. The positive constant d i represents the insulin degradation and clearance rate. In order to give a more accurate account of the oscillatory behaviour, we construct the functions f 1 ; f 2 ; f 4 and f 5 using Hill functions. This modification enables us to introduce new meaningful parameters for which we obtain analytic conditions guaranteeing asymptotic stability. The precise form of these functions is given by ð2:2Þ Function f 5 is the only decreasing function, f 1 ; f 2 and f 4 are all increasing. This implies that h 5 is negative while h 1 ; h 2 and h 4 are all greater than zero. The values for the physiological parameters are chosen as in [11] (see Table 1). We use the following criteria for selecting the Hill coefficients h i ; k i : The steady state of the system must be located in a physiologically relevant range. This is discussed in Section 2.1. The model should be capable of rendering oscillations which match physiological behaviour. An example of an oscillatory regime using this model is shown in Fig. 2. Our numerical simulations demonstrate that our model is also capable of reproducing sustainable oscillations for different values of the rate of insulin degradation d i within its physiological range (Fig. 3).
We note that the maximal value for the derivative of a Hill function can be written explicitly in terms of the Hill parameters. Indeed, for a general Hill function where c > 0; h; k are constants. A direct calculation shows that the maximal value for its first derivative is given by where the principal value has to be used whenever hÀ1 h or hþ1 h is not an integer.

Sensitivity of the steady state on Hill parameters
Although an explicit expression cannot be obtained, system (2.1) possesses a unique steady state ðG Ã ; I Ã Þ, defined by the equations _ G ¼ _ I ¼ 0, which is always positive as was shown in [13]. Since its value is influenced by the choice of the model parameters, the Hill coefficients h i are constrained by the requirement that ðG Ã ; I Ã Þ be located in a physiologically feasible range. In order for the system to fit physiological values the fixed point for glucose, G Ã , should be within the range, 90 < G Ã < 120, and the fixed point for insulin, I Ã , should fit the range, 25 < I Ã < 40 [4,11]. Here, we will investigate the general effect of the Hill parameters on both fixed points by altering the coefficients h 1 and h 4 as the corresponding nonlinear functions f 1 and f 4 are linked to physiological mechanisms related to Type I and Type II diabetes respectively. As function f 2 represents glucose utilisation independent of insulin, the parameters for this function are assumed to remain relatively constant. Fig. 4 shows the effect on the fixed points of parameter h 1 while also altering h 4 , which is represented by the dots.  It is clear from Fig. 4 that there is a large effect of the value of h 1 on the fixed points, particularly for G Ã , and that there is a relatively consistent and small range of values for h 4 for all values of h 1 . A similar result can be seen when comparing h 1 and h 5 . This numerical investigation indicates that choosing the Hill parameter h 1 in a neighbourhood around 2 is necessary in order to keep the steady state within the physiological range. A value of h 1 ¼ 2 had already been used in a number of simplified models in the literature, for example in [18].

Local stability and oscillations
We now focus on describing analytically the boundary of the local stability region in the case when both delays are nonzero. This extends the studies in [10][11][12][13] where several aspects of the local stability of model (2.1) were investigated with very generic assumptions on functions f 1 À f 5 . In matrix form, the linearisation of (2.1) at its unique steady state ðG Ã ; I Ã Þ is given by where we have introduced the following positive quantities The characteristic equation of (3.1) is The analysis of the exponential polynomial (3.2) when s 1 s 2 ¼ 0 can be performed using the circle method, as done in [10,11], leading to individual thresholds s 2 > s Ã 2 (s 1 > s Ã 1 ) for the apparition of oscillatory solutions when s 1 ¼ 0 (resp. s 2 ¼ 0). In the case when s 1 s 2 > 0, it was shown in [12] that if d i A P DðB þ CÞ, then for any given s 1 P 0, there exists s 2 ðs 1 Þ P 0 such that a Hopf bifurcation occurs at ðs 1 ; s 2 Þ. The separatrix curve can be described by seeking imaginary solutions k ¼ ix of (3.2), where x is assumed to be strictly greater than 0. This gives the system Using trigonometric identities, system (3.3) implies that ð3:4Þ so that x satisfies the following transcendental equation Solving (3.5) for sinðs 1 xÞ and introducing into (3.4), we obtain the expression which depends on s 1 only implicitly through x. Therefore, for a fixed value of s 1 , the transcendental Eq. (3.5) can be numerically solved for x whilst the critical value for s 2 is thereafter recovered from ð3:7Þ Using this approach, we obtain a graphical representation of the curve delimiting zones of asymptotic stability and oscillatory regime. Our study shows that the boundary is not a straight line (Fig. 5 Given a fixed value of s 1 , the threshold curve provides the value s 2 ðs 1 Þ such that a Hopf bifurcation occurs at ðs 1 ; s 2 Þ. This is illustrated in Fig. 6, where the distribution of eigenvalues for a specific choice on the threshold curve clearly shows a pair of eigenvalues sitting on the imaginary axis. It is interesting to note that within the physiological ranges for the delays, there is only one pair of eigenvalues that crosses the imaginary axis. Mathematically, model (2.1) does allow to have interaction of unstable modes, but this seems to occur only for very large values for the delays. For instance, when s 1 ¼ 6, a numerical investigation indicates that one needs to have s 2 > 182 in order to have more than one pair of eigenvalues with positive real parts.
It is known that the period of the oscillations depends on several physiological parameters, especially the amount of glucose infused continuously [4,11]. Under continuous glucose infusion, the period of oscillations is generally accepted to lie within the range ½90; 150 minutes. On the threshold curve (3.5), the influence of delays on the period of oscillations can be computed numerically for physiological values for s 1 and s 2 . When G in ¼ 1:35 mg/ml, this gives periods within the range ½134; 155 minutes, which is in the upper part of the physiological range mentioned above. With a reduced amount of infused  glucose, G in ¼ 0:54 mg/ml, the oscillations on the threshold curve have period between 96 and 100 min. The effect of the delays on the period of oscillations is shown in Fig. 7.
Moreover, we can show that the frequency x is maximal when s 2 ¼ 0. The derivative of x with respect to s 1 is given by Since x is assumed to be positive, this shows that x attains locally maximal values at points where sinðs 2 xÞ ¼ 0. The first positive value of s 2 for which this occurs is s 2 ¼ 0. This leads to the expression which coincides with the one-delay case. In that case, the threshold value s Ã 1 for the apparition of oscillations is obtained through the following system The fact that x is maximal when s 2 is zero implies that the delay in the glucose production by the liver leads to a longer period in the regulation loop. Similarly when s 1 ¼ 0, the threshold s Ã 2 is given by

ð3:15Þ
This recovers the results for the one-delay subcases which have been thoroughly investigated in [10,12].

Local asymptotic stability
The boundary of the asymptotic stability region can also be approached through the construction of an appropriate Lyapunov functional for the linear system (3.1) which is further studied numerically by making use of Linear Matrix Inequalities (LMIs). The approach adopted here uses a method for systems of delayed differential equations with two delays developed by He et al. [20]. The basis of the method is to use free weighting matrices to convey the relationships of Leibniz-Newton form. This allows to take into account the relationship between the two delays. If there exist the free weighting matrices, necessary to satisfy the criteria for the set of LMIs, then the linear system with delays is asymptotically stable.
Following [20], we denote x ¼ ½u v T and consider the following Lyapunov functional It was shown in [20] that the requirement that V be a Lyapunov-Krasovskii functional for the linear system (3.1) is fulfilled if the following linear matrix inequalities are satisfied The free weighting matrices P and Q i , i ¼ 1; 2; 3, are symmetric definite positive, W i X ii ; Y ii and Z ii , i ¼ 1; 2; 3, are symmetric positive semi-definite and the free matrices N i ; S i ; T i ; X ij ; Y ij and Z ij , ij ¼ 1; 2; 3, have no restrictions. This method was applied for the linearised system (3.1) using the parameter values as given in (3.8). The solution scheme, which relies on interior fixed point methods, was implemented using the Robust Control Toolbox in Matlab [21]. As a result, a plot of the region of asymptotic stability in the ðs 1 ; s 2 Þ plane has been obtained. The results are depicted in Fig. 8 where an almost linear pattern can be seen distinguishing between the area of asymptotic stability and the area of oscillations. It shows that the LMI optimisation based on the Lyapunov functional (3.16) allows to obtain better results in the two-delay setting than in its onedelay limiting subcases. It appears to be remarkable that the numerical scheme approximates the stability region very accurately in the physiological range when 5 < s 1 < 20.

Global stability and Hill parameters
We now proceed to the construction of a Lyapunov functional for the nonlinear two-delay model (2.1) which allows us to obtain delay-dependent stability conditions. Following the work of Bennett and Gourley [13], we introduce the shifted variables so that ðu; vÞ ¼ ð0; 0Þ represents the steady state of the translated system For a given function f ðxÞ, the Taylor expansion with Lagrange remainder at the first order can be written as where c is some number between x 0 and x. Writing c ¼ x 0 þ hx, with 0 6 h 6 1, the Taylor expansion can be written as Therefore, there exist 0 6 h i 6 1; i ¼ 1; . . . ; 5, such that the following hold ð4:2Þ We now consider the function V 0 ¼ 1 , with q a strictly positive constant. By making use of expressions (4.2) and of Young's inequality AExy 6 1 2 which holds for all > 0, we can show that the derivative of V 0 along the solutions satisfies the inequality which can be obtained through an application of the Leibniz rule for the derivative of an integral with variable limits. This enables us to write dV 0 dt 6 Àu 2 f 2 0ðG Ã þ h 2 uÞ þ qf 4 ðI Ã Þ À 1 2 ð4:6Þ Expression (4.6) motivates the following definition and leads to and likewise for m 2 ðt À s 2 Þ, we can define the Lyapunov functional to be The functional V is clearly positive for all functions u; v and our derivation shows that the following inequality holds This proves the following theorem.
Theorem 1. The nonlinear system (2.1) is globally asymptotically stable if there exist strictly positive constants q; i ; i ¼ 1; . . . ; 9, such that the following inequalities are satisfied: ð4:8Þ ð4:9Þ We note that these conditions constitute a direct generalisation of those obtained in the one-delay case in [13]. To maximise the stability results, we set 4 ¼ 8  Fig. 9 a graph of the variation of f 4 ðI Ã Þ=f 0M 1 with respect to h 1 . This shows that higher values of h 1 will produce a smaller region of asymptotic stability. In accordance with recent analytical and experimental results [18], our analysis indicates that values of h 1 in the vicinity of 2 will provide a model capable of describing both asymptotically stable and oscillatory behaviours.  Following our analysis, we are in a position to suggest ranges on Hill parameters which allow an accurate modelling of the ultradian oscillations in the glucose-insulin regulatory system. These are reported in Table 2. Note that our simulations focussed on a range of ½À9; À7 for h 5 , but allows for a further expansion.

Conclusion
In this paper we introduce a new set of functions to be used in a model of the ultradian oscillations in the glucose-insulin regulation system. The four main functions used in glucose-insulin regulation models are updated by using a Hill function form which involves parameters with physiological meaning, while also removing the need for auxiliary constants. We show that the model with Hill functions is capable of accurately representing the ultradian oscillatory behaviour of the system.
An analytical study of stability of the two-delay system is performed, adopting both a local and global perspective. Existing literature does not cover analytical work related to the two-delay model in a complete sense. This paper extends the work of Li and Kuang, [12], by providing conditions for the apparition of oscillations in the local case when both time delays are nonzero. We analytically described the boundary of the stability region in the ðs 1 ; s 2 ) plane for which oscillations occur. The period of the oscillations were shown to fit the physiological ranges expected for ultradian rhythms under constant glucose infusion. Furthermore, it is very sensitive to both the physiological and the Hill parameters, in particular to glucose infusion. The analytical work is then investigated numerically using a LMI approach [20] which confirms the joint role of the delays in producing oscillations. With the goal of determining delay-dependent sufficient conditions for global stability, a Lyapunov functional approach is taken and applied to the two-delay model. By numerically investigating the stability region dependent on the delays, optimised ranges on the Hill coefficients of the model are introduced. This was done in conjunction with the need to satisfy the physiological ranges of the system. This work lays the foundations for the assessment of physiological data and could in the long term lead to an improvement in the treatment of diabetes, through for example an artificial pancreas, by providing a better understanding of the behaviour of the glucose-insulin regulation system.