Amplitude and Frequency Variation in Nonlinear Glucose Dynamics with Multiple Delays via Periodic Perturbation

Characterising the glycemic response to a glucose stimulus is an essential tool for detecting deficiencies in humans such as diabetes. In the presence of a constant glucose infusion in healthy individuals, it is known that this control leads to slow oscillations as a result of feedback mechanisms at the organ and tissue level. In this paper, we provide a novel quantitative description of the dependence of this oscillatory response on the physiological functions. This is achieved through the study of a model of the ultradian oscillations in glucose-insulin regulation which takes the form of a nonlinear system of equations with two discrete delays. While studying the behaviour of solutions in such systems can be mathematically challenging due to their nonlinear structure and non-local nature, a particular attention is given to the periodic solutions of the model. These arise from a Hopf bifurcation which is induced by an external glucose stimulus and the joint contributions of delays in pancreatic insulin release and hepatic glycogenesis. The effect of each physiological subsystem on the amplitude and period of the oscillations is exhibited by performing a perturbative analysis of its periodic solutions. It is shown that assuming the commensurateness of delays enables the Hopf bifurcation curve to be characterised by studying roots of linear combinations of Chebyshev polynomials. The resulting expressions provide an invaluable tool for studying the interplay between physiological functions and delays in producing an oscillatory regime, as well as relevant information for glycemic control strategies.


Introduction
Homoeostasis refers to the body's ability to maintain certain variables within a narrow range and is a result of the negative feedback loops which occur within the body (Cannon 1932;Thomas et al. 1995). These loops keep the parameters involved at, or close to, a healthy value, or steady state (Thomas et al. 1995). The steady state may be stable, leading to constant levels over time, or unstable, leading to sustained oscillations (Thomas et al. 1995). Examples of homoeostatic processes include: the regulation of temperature; blood pressure; glucose concentration; percentage of water within the cells; and calcium levels. While many biological systems can be modelled with ordinary differential equations (ODEs), delays can so often prove crucial in realistically replicating core aspects of these systems (Thomas et al. 1995). Indeed, a lot of work has gone into determining the role of these delayed recurrent loops in physiology (Bocharov and Rihan 2000), especially in hormonal systems (Walker et al. 2012(Walker et al. , 2010. For example, in the case of modelling the glucose-insulin regulatory system, an ODE system with subsystems replicating delay can replicate the ultradian rhythms that occur naturally in the system without the need for postulating an internal pulsatile insulin pacemaker (Levy 2001;Sturis et al. 1991).
This paper focuses on the ultradian rhythms that occur within the glucose-insulin system. As discussed in O' Meara et al. (1993), insulin resistance has been observed to gradually dampen the oscillations, culminating in a loss of synchronisation between glucose and insulin. While insulin resistance may be present in patients without diabetes, it is generally associated with patients with type 2 diabetes (T2DM), and so we take the accurate tuning of the ultradian oscillations to be a sign of healthy regulation (Huard et al. 2017). We then propose the following question: What is the effect of diabetic parameters on the amplitude and period of the ultradian rhythms?
To answer this question, we introduce a perturbative scheme for the periodic solutions of a two-compartment delay differential equation (DDE) model of the ultradian rhythms in the glucose-insulin regulatory system based on the Poincaré-Lindstedt (P-L) method. The model is a polynomial expansion of the system presented in Huard et al. (2017), which was originally created in Sturis et al. (1991). A large number of authors developed and studied this model to characterise its local and global stability properties and characterise its periodic solutions (Bennett and Gourley 2004;Engelborghs et al. 2001;Giang et al. 2008;Huard et al. 2015;Kissler et al. 2014;Li and Kuang 2007;Li et al. 2006;Wang et al. 2009). As far as we are aware, this is the first time that the P-L method has been applied to a model of glucose and insulin regulation. It has been applied to other biological systems such as in Verdugo and Rand (2008), where it was used to predict the amplitude and frequency for a two-compartment DDE model of gene expression with one delay. Similarly, in Brandt et al. (2006), the frequency of a two-compartment DDE model for describing two couple Hopfield neurons was calculated. To the best of our knowledge, in all cases where this technique has been applied to dynamical systems with n components and multiple delays τ i 's, the characteristic equation defining the eigenmodes λ's is an exponential polynomial of the form c n λ n + c n−1 λ n−1 + · · · + c 0 + e −λ i τ i = 0, c i ∈ R, τ i ∈ R + , λ ∈ C, (1) in which the polynomial part is typically Hurwitz stable. Thus, for these systems, delays interact in a linear manner to give rise to a supercritical Hopf bifurcation at a point in the space of delays which can be obtained, along with the characteristic frequency, by using imaginary root crossing methods (Cooke and Van Den Driessche 1986;Kuang 1993).
The system studied in this paper, which is presented in Sect. 2, contains two delays, τ 1 and τ 2 , and poses the additional challenge that the characteristic equation involves two exponentials. Although crossing curves can be described in this situation using geometric arguments (Gu et al. 2005), we show that characteristic frequencies can be obtained by studying the zeroes of linear combinations of Chebyshev polynomials of the first type. This is achieved by assuming the commensurateness of delays, that is, τ 2 = κτ 1 , where κ ∈ Z + is a coupling parameter. This assumption does not lead to additional restrictions and allows us to define a group of lines in the space of delays along which the perturbative scheme can be defined. Furthermore, it enables the study of solutions for physiologically relevant ranges of the model parameters.
In summary, the paper aims to provide a way of quantifying the contribution of model parameters to the amplitude and period of the ultradian rhythms in the glucoseinsulin regulatory system. An analysis of the effect of each physiological feature on the limit cycles, with a particular focus on insulin resistance, is then described. The work is divided as follows. After presenting the model in Sect. 2, local stability properties and conditions for the boundedness of trajectories are detailed in Sect. 3. In Sect. 4, the P-L method is applied in order to obtain approximate formulas for the amplitude and period of the oscillations. Section 5 is dedicated to the study of the influence of model parameters on the oscillations. Finally, physiological implications are explored, along with concluding outlooks.

Model
The model under consideration follows the framework given in Fig. 1 and is a twocompartment first-order polynomial DDE system given bẏ where G(t) is the plasma glucose concentration (in mg) and I (t) is the plasma insulin concentration (in mU ). By noting that, within the plasma, the volume of glucose space is around 100dl while that of insulin distribution is approximately 3l (Sturis et al. 1991), G(t) and I (t) are converted to concentrations (mg/dl and mU /l, respectively), for the purpose of all figures. All parameters are assumed to be nonzero and have units as defined in Table 1. The coefficients a 1 and a 2 describe the insulin independent and dependent glucose utilisations respectively. The parameter a 0 = G in + C includes two contributions, namely G in which corresponds to a constant glucose infusion, and the constant C which partially represents the hepatic glucose production. Indeed, following the work of Huard et al. (2015), Li et al. (2006) and Sturis et al. (1991), we note that hepatic glucose production can be represented by a rational function of the type where τ 2 corresponds to the time taken for variations in insulin levels to have an observable effect on hepatic glucose production, in minutes. In the insulin balance equation, b 1 is the insulin production capability of the individual, and b 2 is the insulin degradation rate. While the insulin secretion term appears to be unbounded, we note that it is usually modelled using a sigmoidal function (Huard et al. 2015;Li et al. 2006;Sturis et al. 1991), where τ 1 is the time taken for insulin to be secreted and transported to the interstitial space in response to elevated glucose levels. Hence, the secretion function is a local approximation of a bounded function. Thus, although the polynomial form of the system contrasts with previous models in which the hepatic and pancreatic secretions are represented by bounded sigmoidal functions (Bennett and Gourley 2004;Huard et al. 2015;Li et al. 2006), it is appropriate for studying dynamics in the neighbourhood of the limit cycle. To further aid in the study of the model dynamics, we shall later take the delay τ 2 to be commensurate with τ 1 . It is worth noting that the convergence speed to this limit cycle may be altered through this approximation. Boundedness and positivity of trajectories can also be affected through the polynomial approximation, as investigated in Sect. 3.2.1.

Local Stability Analysis
The presence of Hopf bifurcations in model (2) strongly depends on the values of delays. We divide our analysis into two parts. Firstly, we investigate conditions for the one-delay model (with a 3 = 0) to undergo a Hopf bifurcation. Secondly, the bifurcation curve in the two-delay model is investigated by assuming commensurateness of delays.

Constant Hepatic Glucose Production
Here, we look to find conditions on parameter values in the polynomial model (2) with a 3 = 0 ensuring the presence of oscillations in an appropriate range. Equivalently, this amounts to investigating the second-order approximation of (2) when p > 2. In this case, only the constant term from hepatic production remains, namelẏ Although τ 2 does not appear in the resulting system, the minimal model (3) captures the role of the delay τ 1 in producing an oscillatory regime. Incidentally, it was observed numerically that the contribution of τ 1 to the amplitude and period of the oscillations is more prominent than that of τ 2 (Li et al. 2006). Furthermore, the positivity and bound-edness of trajectories for system (3) is easily demonstrated following the arguments formulated in Bennett and Gourley (2004) and Shi et al. (2017).

Value of Parameters
Given an oscillatory solution, the inverse problem of choosing parameters in model (3) can be addressed in the following way. We note that the system's steady state (Ḡ, I ) satisfies the equations By Descartes' rule of signs, (4) has exactly one positive root forḠ, and so one can always find a positive (Ḡ,Ī ). Here we assume that the target basal levels can be identified with the steady state of the system.
1. It has been shown that the insulin clearance is proportional to the plasma insulin concentration (Koschorreck and Gilles 2008;Topp et al. 2000). Numerical fitting procedures have rendered values for b 2 in the range (0.03, 0.3) (Chen et al. 2010). As in Huard et al. (2017), Huard et al. (2015) and Li et al. (2006), we shall choose an initial value of b 2 = 0.06 for most numerical computations with the reduced model. The value of b 1 is then obtained as b 1 = b 2ĪḠ −n . 2. The system must be in an oscillatory state. Therefore, the delay τ ∈ R + must be larger than the critical value τ 0 such that the characteristic equation possesses a set of conjugate purely imaginary root λ = ±iω 0 . This requirement implies that the critical pair (ω 0 , τ 0 ) satisfies the system This in turns implies that ω 0 satisfies a quartic equation Requiring that Eq. (7) possesses a positive root for ω 0 gives explicit conditions on the coefficients for the existence of a Hopf bifurcation. Since the middle term in (7) is always positive, one must have n > a 1 a 2Ī + 1 > 1, in order to have a bifurcation. Consequently, the periodic perturbation scheme presented in Sect. 4.1 shall focus on the case n = 2, which implies that a 1 < a 2Ī . The value of τ 0 is then obtained from (6) as where K ∈ Z + is the smallest integer such that (8) defines a positive value. Larger values of K give successive τ 0 values for which stability switches may occur in the linear system. However, for the nonlinear system (3), it is numerically observed that oscillations are present whenever τ ≥ τ 0 . 3. The constant a 0 = G in +C is then obtained from the steady-state equation, a 0 = a 1Ḡ + a 2ḠĪ . The parameters a 1 and a 2 must be chosen such that oscillations are present for a physiologically relevant value of the critical delay τ 0 . For different values of b 2 , one can numerically compute the range of achievable values for τ 0 using Eq. (8) (see Figure 9 in Bridgewater et al. 2019).
This approach provides a model able to replicate the nonlinear oscillations within an appropriate physiological range (Fig. 2).

Extension to Two Delays
We now consider the problem of studying periodic solutions in the two-delay system (2). In the first instance, we investigate the question of boundedness and positivity of trajectories of the system. In the second, the characterisation of periodic solutions in the two-dimensional space of delays is achieved by assuming the commensurateness of delays, i.e. τ 2 = κτ 1 , with κ ∈ Z + .

Positivity and Boundedness
Here we look for conditions for the positivity and boundedness of solutions in the polynomial initial value probleṁ where φ, ϕ are continuous strictly positive functions on [− max (τ 1 , τ 2 ), 0], p ∈ Z + and n is a strictly positive even integer, which is the typical case in insulin secretion modelling studies (e.g. Topp et al. 2000;Yang et al. 2016). As mentioned previously, positivity and boundedness of solutions for model (9) when a 3 = 0 can be established using the same techniques as in Bennett and Gourley (2004) and Shi et al. (2017). Here, the additional non-positive and possibly nonlinear term when a 3 = 0 can lead to unbounded trajectories and to the development of singularities in finite time. Nonetheless, the positivity of I (t) can be easily established.
Proposition 1 For all solutions of (9) and ∀t > 0, we have that I (t) > 0. Furthermore, ∀t > τ 2 , the following inequality can be established Proof Indeed, from the second equation in (9), we have which implies the positivity of I (t) since n is assumed to be a positive even integer. Equation (10) also leads to Due to the biological nature of the problem, unbounded solutions are not in the scope for applications, and therefore, we shall focus on parameter values not leading to such solutions. The following holds true.

Proposition 2 If G(t) is strictly positive ∀t > 0, then G(t) and I (t) are bounded from above.
Proof Let G(t) be positive for all times t > 0. Then, since I (t) > 0, we have thaṫ The converse can be established in the following case.
Proof From the first equation in system (12), it can be easily seen that Hence G(t) is strictly positive for all t if a 0 > a 3 I p + and G(0) > 0, and is therefore bounded by Proposition 2.
It is noteworthy to highlight that depending on the values taken for n and p, some solutions may become unbounded in finite time. A subclass of such blow-up solutions is represented by trajectories possessing poles for which the location depends on initial conditions. Such singularities, called movable poles, can be detected in the delay-free system using a local analysis of the Painlevé type (see e.g. Conte 2012; Goriely andHyde 1998, 2000). As shown in Appendix 1, the values n = p = 2 may lead to this kind of singular solution. Consequently, the periodic perturbation scheme presented in Sect. 4.2 shall focus on the case n = 2, p = 1 to avoid unbounded solutions.

Periodic Solutions in System with Commensurate Delays
We now consider the problem of characterising periodic solutions in system (2) where delays are assumed to be commensurate, i.e. τ 2 = κτ 1 , with κ ∈ Z + . A straightforward generalisation can be made for the case when the delays are rationally related. This assumption allows to perform a perturbative analysis of the periodic solutions along the line τ 2 = κτ 1 , given that the point (τ 1 , κτ 1 ) remains sufficiently close to the curve of Hopf bifurcations in the delay space (τ 1 , τ 2 ). Geometrically, this approach provides a discrete set of critical points, corresponding to the intersection between lines and the threshold curve (Fig. 3). We show that these crossing points can be described by studying the zeros of linear combinations of Chebyshev polynomials of the first kind.
In its most general form, model (2) with commensurate delays becomeṡ Any equilibrium (Ḡ,Ī ) of (12) obeys the equations which always possess a unique strictly positive root since a i , b i > 0 and n, p ∈ Z + . The linearisation about this steady state (Ḡ,Ī ) reads as and similarly for y. The characteristic equation of (14) is a quasi-polynomial of the form where (15) is Hurwitz stable for τ = 0 and we now look for conditions on κ and τ ensuring that it undergoes a Hopf bifurcation. Hence, setting λ = iω 0 , ω 0 > 0 leads to the system Eliminating polynomial occurrences of ω 0 , one is led to the following trigonometric equation where T n and U n−1 are Chebyshev polynomials of the first and second kinds, respectively. Equation (18) implies that z satisfies the following polynomial equation Further using the following well-known identities, for n ≥ m, we obtain that (19) can be rewritten as a linear combination of Chebyshev polynomials of the first type For a given κ ∈ Z + , any real root of (20) with |z| < 1 which gives a nonzero solution for ω 0 through equations (16) and (17), which can be rewritten in terms of z as leads to a non-constant periodic solution.
Example 1 Let us consider, for illustration, the special case It is easily seen that the root z = −1/2 is the unique value leading to a null frequency (ω 0 = 0), since U 1 (z) is linear. The presence of an additional root in the interval (−1, 1) can be assessed using, for instance, Sturm sequences (see, e.g. Bochnak et al. 2013). For the polynomial 2z 3 + z 2 − z − 1, the Sturm chain H (z) is given by which, when evaluated at the end points, gives #sign changes(H (−1)) − #sign changes(H (1)) = 1.

Hopf Bifurcation Formulae
We now turn to the main objective of this work, namely the analysis of the variation of amplitude and frequency in the nonlinear model with respect to the model parameters. This is achieved through the perturbation of periodic solutions in the linear model in a neighbourhood of the critical manifold. This technique is usually referred to as the P-L expansion and was extended to differential equations with explicit delay in Casal and Freedman (1980) and applied to a limited number of delayed models to highlight the contribution of parameters to the amplitude of oscillations (Brandt et al. 2006;Casal and Freedman 1980;Rand and Verdugo 2007;Verdugo and Rand 2008). We first consider the case when hepatic glucose production is assumed to be constant (see model (25)), in which case the small bifurcation parameter > 0 represents the distance from the critical τ 0 . Secondly, we extend our considerations to the two delay case by assuming that the second delay is a constant multiple of the first one, thus defining a fan of expansion lines in the space of delays (see model (39)). Throughout this section, it is assumed that the frequency ω 0 of the periodic solution of the linearised system does not satisfy an equation of order lower than the characteristic equation.

Constant Hepatic Glucose Production
Recall from Sect. 3.1 that the simplified model with a constant hepatic glucose production term, with n = 2, is given bẏ Defining X (t) = X = G(t)−Ḡ, Y (t) = Y = I (t)−Ī , as deviations from the positive equilibrium, substituting in (25) and eliminating Y allows us to write the following nonlinear second-order DDE for X We now introduce the bifurcation parameter , which is defined as the distance from the critical delay τ 0 as = √ τ − τ 0 . The variables are scaled as where s corresponds to a new time variable ensuring that u(s) has a period of 2π .
Here, is also assumed to have an -expansion where ω 0 is the frequency associated to the critical value τ 0 . Finally, we expand the delayed term u τ along with u and its derivatives Substituting in (26) and collecting terms (up to and including O( 2 )) gives with u iτ = u i (s − τ 0 ω 0 ). Without loss of generality, the seed solution is chosen as u 0 (s) = A 0 cos(s) where, following (27), A 0 is related to the amplitude of X (denoted byĀ) byĀ = A 0 . Choosing and substituting in (32) and comparing coefficients of the cos(s) and sin(s) terms shows that: (1) A 1 , B 1 are arbitrary (2) ω 1 = 0. From comparison of the cos(2s) and sin(2s) coefficients, it can be shown that where F, G, H and K are some functions of a 1 , a 2 , b 1 , b 2 ,Ḡ,Ī , ω 0 . Due to their length, they are not reproduced here. By substituting u 2 (s) = A 2 sin(s) + B 2 cos(s) + C 2 sin(2s) + D 2 cos(2s) + E 2 sin(3s) + F 2 cos(3s) + G 2 along with lower-order terms into (33), and comparing the coefficients of the cos(s) and sin(s) terms, the dominant term for the amplitude,Ā, of the limit cycle is expressible asĀ while the dominant term of the amplitude of the insulin oscillations,B, and the secondorder term for frequency correction, ω 2 , are given bȳ Here, we introduced the following definitions with b ρ = b 2 2 + ρ and b 4ρ = b 2 2 + 4ρ. The period of the limit cycle, T , is given by Simulations making use of expressions (36), (37), and (38) are presented in Sect. 5.1.

Linear Hepatic Glucose Production
We now turn to study the effect of a non-constant hepatic glucose production, and hence a second delay, on the limit cycles (see model (12)). As described earlier, the P-L method has been used in a limited number of studies to investigate the effect of model parameters on the amplitude of the resulting oscillatory solutions. For example, in Brandt et al. (2006) a coupled first-order DDE model describing a two-neuron system with delay was explored. While the system had two separate delays, these were combined giving a characteristic equation of the form with A 1 , A 2 constant and A 3 a function of the model parameters. In contrast, the characteristic equation of model (9), given by (15), contains a second exponential term which leads to additional challenges in finding points of bifurcation, as discussed in Sect. 3.3. For conciseness and in order to avoid the blow-up of solutions in finite time, we assume that p = 1 although the technique can be applied to higher orders under some restrictions (see Appendix 1). The resulting equations are thus given bẏ where the dimensionless parameter κ is used to represent the commensurateness of the time delays. In line with the calculation in the previous section, we introduce X = G(t) −Ḡ, Y = I (t) −Ī . Let (ω 0 , τ 0 ) be a critical pair obtained using the algorithm described in Sect. 3.3. Then by setting = √ τ − τ 0 , we can scale the variables X and Y as X (t) = u(s), Y (t) = v(s), where s is the scaled time variable as defined in (27). The expansions for v(s) and v(s − κ τ ) are given by Substituting in (39) and collecting terms (up to and including O( 3 )) gives with m = 0, 1, 2, u mτ = u m (s − τ 0 ω 0 ), v mτ = v m (s − κτ 0 ω 0 ) and where the inhomogeneous terms g m and h m are related to the solutions of previous orders. Here we have g 0 = 0, h 0 = 0, and By imposing the initial conditions u 0 (0) = C 0 , v 0 (0) = C 0 R 1 on the solutions of (43) and (44), with m = 0 we find that , Q = 2(ω 0 (a 2Ḡ + a 3 cos(κτ 0 ω 0 )) + a 3 b 2 sin(κτ 0 ω 0 )).
We note that C 0 is related to the amplitude of X (t) (denoted byC) byC = C 0 and that the amplitude of Y (t) (denoted byD) is given bȳ As we necessitate that the solutions u m and v m are periodic in s with period 2π , we must find conditions such that the inhomogeneities do not contain secular terms. Hence, we proceed as in Brandt et al. (2006) and expand u m , v m , g m and h m as Fourier series to get Substituting (50) and (51) in (43) and (44), it can be seen that the coefficients, α (m) j,1 , β (m) j,1 (with j = 1, 2), in the inhomogeneities g m and h m must satisfy the following conditions 2,1 cos(τ 0 ω 0 ) = 0. (53) The system for k = 1 leads to equations of the form 1 (a 1 , a 2 , a 3 , b 1 , b 2 where Z 1 , Z 2 are functions of a 1 , a 2 , a 3 , b 1 , b 2 , κ,Ḡ,Ī , ω 0 , τ 0 . If C 0 = 0, then we obtain the trivial solution. Additionally, it can be seen that Z 1 and Z 2 do not vanish modulo the characteristic curve. Therefore, following our assumption that ω 0 does not satisfy any polynomial equation of lower order, this implies that ω 1 = 0. For further details of the derivation of these conditions, the reader is referred to Appendix 1. Substituting these conditions in (43) and (44) leads to the following expressions where h 1 , h 2 , h 3 and h 4 are functions of a 1 , a 2 , a 3 , b 1 , b 2 , κ,Ḡ,Ī , ω 0 , τ 0 . Expressions for α (2) 1,1 , β 1,1 , α 1,2 and β 1,2 can then be obtained. Finally, using conditions (52), (53) and solving the resulting system of equations, it can be shown that the dominant term for the amplitude,C, the second-order term for frequency correction, ω 2 , and the period, T , can be expressed as where W 1 , W 2 , V 1 and V 2 are functions of a 1 , a 2 , a 3 , b 1 , b 2 ,Ḡ,Ī , κ, ω 0 and τ 0 . Given their length, the expressions are not reproduced here but are used in simulations in Sect. 5.2.

Parameter Analysis
The closed-form expressions for the limit cycles presented in this paper allow for the effect of changes in each of the model parameters on the amplitude and period to be more easily studied. In this section, we first study the effect of changes in a 2 (insulin resistance), b 1 (insulin secretion) and b 2 (insulin degradation) on the important characteristics of the waveforms for the constant hepatic production model solutions. Secondly, we investigate how the amplitude and period of the solutions of model (39) vary with respect to the delay coupling parameter κ.

Constant Hepatic Glucose Production
In this section, the closed-form expressions for the amplitude and period of the solutions to model (3), as given by (36), (37) and (38), will be analysed using two different Parameter Sets, which are given in Table 2. The values used in Parameter Set 1 are based off the values used in Huard et al. (2015), Li et al. (2006) and Sturis et al. (1991) and represent a patient under a constant glucose infusion, while Parameter Set 2 looks to replicate the values that may be observed in a patient under a larger constant glucose infusion.

On the Influence of Insulin Sensitivity
To begin, we analysed the relationship between the insulin sensitivity parameter a 2 and the closed-form expressions for the amplitude and period of model (3). It can be seen Fig. 4 Amplitudes of the oscillations,Ā andB, as a function of a 2 using Parameter Set 1 (left) and Parameter Set 2 (right) Fig. 5 Period of the oscillations, as given by (38), as a function of a 2 using Parameter Set 1 (left) and Parameter Set 2 (right) from Fig. 4 that the amplitude of X (t), as given by (36), varies between 4 and 15 for Parameter Set 1, and 1 and 38 mg/dl for Parameter Set 2. Additionally, the amplitude of Y (t) from (37) is observed to be approximately between either 1 and 5, or 1 and 29 mU /l. These values are within a physiologically acceptable range. Furthermore, in Fig. 5 we note that the values of the period, as defined by (38), vary between 74.5 and 76 minutes for Parameter Set 1, and 76 and 78 min for Parameter Set 2, and hence are also within an acceptable range. It can be seen in Fig. 6 that increasing a 2 has little effect on the oscillations, while decreasing a 2 has a more profound effect. For example, for Parameter Set 1 a 20% increase in a 2 from 0.0017 increases the amplitude by 15%. However, a 20% decrease in a 2 reduces the amplitude by almost 80%. Similarly for Parameter Set 2, a 20% increase in a 2 from 0.0004 increases the amplitude by less than 5% while a 20% decrease in a 2 reduces the amplitude by approximately 30%.

Insulin Secretion Capacity b 1 and Insulin Degradation b 2 versusĀ andB
Next, we looked at the relationship between the insulin secretion capacity b 1 and the closed-form expression of the amplitude of X (t) defined by (36). As shown in Fig. 7a, the amplitude variation with respect to b 1 is between 0 and 19, regardless of the value of a 2 used. However, a 2 does have an effect on the decline of amplitude observed with an increased b 1 . Indeed as a 2 increases, the observed value of b 1 such that the amplitude begins to decrease decreases. This effect is also observed in Fig. 7b.  Fig. 6 Percentage change of the closed-form expressions of the amplitude (blue) and period (black) of Y (t), which are given by (37) and (38) respectively, versus the percentage change in a 2 for Parameter Set 1 (left) and Parameter Set 2 (right). The initial value used for a 2 was: 0.0017 (left); and 0.0004 (right) Figure 7c, d show the effect of b 2 on the amplitude of X (t) and Y (t) respectively for Parameter Set 1. It is observed from both that the oscillations are lost when b 2 drops below ≈ 0.059 regardless of the value of a 2 . Indeed, when looking at the amplitude as a function of b 2 , a 2 has very little effect onĀ and only a small effect onB for Parameter Set 1. However, as seen in Fig. 7e, f, this is not true for Parameter Set 2, where a 2 has a more profound effect on both the amplitude of X (t) and Y (t). Irrespective of this, it is observed for both Parameter Sets that an increase in b 2 leads to an increase in the amplitude of the oscillations. While the size of the oscillations of Y (t) in Fig. 7d, f is plausible for all values of b 2 , the size of the oscillations in Fig. 7c, e becomes too large with b 2 . Indeed, eventually the value ofḠ −Ā drops below 70 mg/dl in both cases. Physiologically, this would mean the onset of hypoglycaemia. The ranges 0.06 < b 2 < 0.062 for Parameter Set 1, and 0.064 < b 2 < 0.075 for Parameter Set 2 ensure that glucose levels are kept within a realistic physiological range.

Non-constant Hepatic Glucose Production
We now move on to investigate the effect of the two delays on the closed-form expressions for the amplitude and period of model (39). In particular, we focus on the relationship between the amplitudes of X (t) and Y (t), given by (58) and (49) respectively, and the commensurate delay parameter κ. Using Parameter Set 1 with a 2 = 0.0017, it can be seen in Fig. 8 that when a 3 is small, κ has a negligible effect on the amplitude. Additionally, when a 3 < O(10 −2 ), changes in a 3 also have a negligible effect on the amplitudes. However, when a 3 = 0.1 it is observed that the amplitude of both X (t) and Y (t) varies with κ in an oscillatory manner.
To further investigate how the closed-form expressions for the amplitude and period vary with κ, we now look into the variation using the parameter values given in Table 3. From Fig. 9, we can see that κ has a large influence on both the amplitude and period for Parameter Set 3. Increasing κ from 1 to 7, the amplitude increases by ≈ 700% and the period by ≈ 500%. However, it must be noted that for κ > 1 the values of the amplitude and period are not in a physiological range. This is most likely due to the size of 2 . Indeed, when κ = 7 we note that 2 = 9.49476 compared with 0.2806 when κ = 1. Therefore, instead of setting τ = 16 we shall instead choose τ such that Fig. 7 Amplitude of X (t), as given by (36), and Y (t), as given by (37), versus b 1 and b 2 for Parameter Sets 1 and 2 2 = 0.16. The results can be seen in Fig. 10. Here we observe that the amplitude remains in a physiological range for 1 ≤ κ ≤ 7 and that the period of the oscillations is within an acceptable physiological range for insulin levels. We also note that there is very little variation in the amplitude when κ increases. This implies that more accurate values of the two delays can be chosen without losing the physiological accuracy of the oscillations.  (58) and (49) respectively, as a function of κ for Parameter Set 1 with a 2 = 0.0017

Discussion
In this contribution, we have analysed two polynomial DDE models of the glucoseinsulin regulatory system, one with a constant hepatic glucose production and one with a linear hepatic production containing a commensurate delay, to investigate the effect of various diabetic parameters on the ultradian oscillations. Indeed, by performing a P-L perturbation method we have been able to obtain analytical expressions that are based on the model parameters for the linearised amplitude and period of the two models. From a mathematical point of view, the accuracy of these expressions is of a high degree. Indeed, from Fig. 11 we can see that the closed-form expression for the amplitude of X (t) given by (36) is almost an exact match for the amplitude obtained through numerical simulations. Furthermore, Fig. 12 shows that the first-order solutions for G(t) and I (t) of (3) obtained using the P-L technique are a good approximation to the solutions calculated using a classical Runge-Kutta method. Increasing the accuracy by taking into account the second and third-order terms computed in Sect. 4 leads to an approximate solution that is indistinguishable from the one obtained numerically.
From a physiological perspective, it is important to note while the values of the amplitude of X (t) in Fig. 4 are within a physiologically acceptable range, the decrease in amplitude in the presence of mild insulin resistance is most notably seen experimentally in insulin levels, while the amplitude in glucose oscillations has been observed to remain almost constant (O'Meara et al. 1993). Therefore, the expressions derived in Fig. 9 Amplitude of X (t) (left) and Y (t) (right), as defined by (58) and (49) respectively, versus κ for Parameter Set 3

Fig. 10
Amplitude of X (t) (left) and Y (t) (right), as given by (58) and (49) respectively, versus κ for fixed . All other parameters are as defined in Table (3) Fig. 11 Comparison of two methods for calculating the amplitude of X (t) from model (3) using Parameter Set 1. The blue line represents the P-L approximation given by (36), and the black dots represent the numerical approximations obtained using a Runge-Kutta method this paper could theoretically be used to obtain estimates for insulin sensitivity through the matching of clinical insulin data to the two models. Finally, we note that strategies aiming to restore glucose-insulin oscillations could make use of the fact that the amplitude and frequency of oscillations show very little A.Ġ ∼ −a 2 G I B.Ġ ∼ −a 3 I p , p > 1.
Looking for dominant terms of the form one is led to the following solution for each truncation A. G(t) ∼ 2 1/n −1 The following conclusions can be drawn.
Truncation A. Since all parameters are assumed to be strictly positive, we see that when n is even, equations (61), (62) do not lead to an expression with real coefficients. Therefore, no open set of initial conditions leads to a movable pole (Goriely andHyde 1998, 2000). Truncation B. In order to have a pole when p > 1, the quantities z 1 = 1+ p np−1 and z 2 = 1+n np−1 need to be positive integers, which can only happen when n = p = 2. Indeed, since n is assumed to be even and p > 1, the quantity np − 1 is larger than 1. Hence, for z 1 and z 2 to be positive integers, one needs to have the inequalities which can be rewritten as thus leading to n = p = 2. Under these conditions, it is easily seen that the coefficients defined by equations (63), (64) are real. Further looking for additional powers in the expansion allowing for the presence of arbitrary constants (socalled resonances), we obtain the following values for r , where r 0 = −1 corresponds to the arbitrariness of t 0 . The coefficient r 1 being a positive integer by construction, it provides the order at which an additional arbitrary constant may appear, which is necessary to satisfy the initial value problem.
Thus, it is shown that Truncation B (with n = p = 2) may possess finite-time blowup solutions. Of course, these considerations do not exclude the potential presence of other types of singularities for other values of n, p.