1. Introduction
Cyclic rhythms are widely recognized for their significant role in regulating the function of biological and physiological systems [Reference Goldbeter and Yan16, Reference Keener and Sneyd28]. Endogenic oscillations are typically encountered in healthy individuals, while a progressive lack of control of these rhythms is often associated with system stress (for example, sleep deprivation [Reference Sweeney, Jeromson, Hamilton, Brooks and Walshe50, Reference Sweeney, Peart, Ellis and Walshe51]) and disease evolution in humans [Reference Spiga, Walker, Terry and Lightman45].
A prominent example of such endocrine oscillations in the human body is the self-regulation of blood glucose levels [Reference Grant, Wilsterman, Smarr and Kriegsfeld19]. When blood glucose levels increase, insulin is released from the pancreas. Insulin then causes blood glucose levels to decrease by stimulating body cells to absorb glucose from the blood. Conversely, when blood glucose levels fall, pancreatic
$\alpha $
-cells release glucagon stimulating hepatic glycogenolysis and neoglucogenesis. The level of blood glucose is then controlled by the rates of insulin secretion (activation by glucose) and hepatic glucose production (inhibition by insulin). Within the glucose–insulin regulatory system, both rapid oscillations of insulin (period
$\sim $
6 to 15 minutes) and ultradian oscillation of glucose and insulin (of similar period
$\sim $
100 to 150 minutes [Reference Scheen, Byrne, Plat, Leproult and Van Cauter42, Reference Sturis, Van Cauter, Blackman and Polonsky49, Reference Van Cauter, Désir, Decoster, Fery and Balasse52]) have been observed during fasting, meal ingestion, continuous enteral and intravenous nutrition [Reference O’Meara, Sturis, Van Cauter and Polonsky38].
The most important pathway to understanding the underlying mechanisms of these glucose–insulin oscillations is measuring the response to glucose infusions. A large quantity of metrics and mathematical models have been devised for that purpose. While the HBA1c metric remains an essential tool for the diagnosis, prevention and control of Type 2 diabetes [2], clinical tests involving patterns of glucose intake combined with mathematical models provide a mechanism for evaluating the efficacy of internal regulation [Reference Ajmera, Swat, Laibe, Le Novere and Chelliah1, Reference Huard and Kirkham24, Reference Makroglou, Li and Kuang34, Reference Palumbo, Ditlevsen, Bertuzzi and De Gaetano39]. The minimal model devised by Bergman and Cobelli [Reference Bergman4, Reference Bergman, Ider, Bowden and Cobelli5] provides an effective method for estimating insulin sensitivity from an intravenous or oral glucose tolerance test, although it can lead to underestimation in individuals with a large acute insulin response [Reference Ha, Muniyappa, Sherman and Quon20]. With the wider availability of continuous glucose monitors and automated insulin pumps, the ability to detect diabetic deficiencies relies on the capacity of models to reproduce more complex and realistic dynamics under various routine life conditions such as, for example, sleep deprivation [Reference Sweeney, Jeromson, Hamilton, Brooks and Walshe50]. These new capabilities potentially allow for less-invasive test strategies that do not require the detection of high-glucose levels in response to a test.
The main goal of this article is to identify the types of behaviours in a suitable mathematical model that can be expected as a response to periodic glucose uptake, specifically periodic on–off glucose infusions, which can be implemented in practice as repeated lower-dose intravenous injections (see, for example, [Reference Streeten, Gerstein, Woolfolk and Doisy46, Reference Verdy, Roussy, Tetreault and Leboeuf53]).
We focus on the capacity of the system to fall into lockstep with the frequency of a sufficiently strong glucose stimulus (so-called entrainment) which has been observed in numerous contexts at the ultradian and circadian levels in endocrinology [Reference Huard and Kirkham24, Reference Walker, Terry and Lightman54, Reference Zavala, Wedgwood, Voliotis, Tabak, Spiga, Lightman and Tsaneva-Atanasova56], but especially in models of glucose–insulin oscillations with periodic infusion [Reference Sturis, Knudsen, O’Meara, Thomsen, Mosekilde, Van Cauter and Polonsky47]. In particular, we will keep track of the peak glucose values exhibited by the model, allowing us to describe dynamical features which do not require high glucose doses.
Many efforts have been made to replicate the nonlinear response of the glycolytic system; in particular, the mathematical modelling of the delayed response of individual parts of the system by explicit time delays has proven an effective means to explain the onset of self-sustained, ultradian oscillations in the glycemic system [Reference Chen, Tsai and Wong7, Reference Cohen and Li8, Reference Li, Kuang and Mason32]. A common approach to modelling oscillatory behaviour of complex biology is to consider time delays [Reference Glass, Jin and Riedel-Kruse15]. In particular, models of endocrine regulation often incorporate explicit delays to account for the time required for the synthesis, release and action of hormones or metabolites [Reference Walker, Terry and Lightman54]. Various models of intrapancreatic rhythmic activity have been proposed recently, see [Reference Huard and Kirkham24] for a review. For example, it was shown that glucose oscillations can enhance the insulin secretory response at the
$\beta $
-cell level when tweaked at a suitable amplitude and frequency [Reference McKenna, Dhumpa, Mukhitov, Roper and Bertram35]. Negative delayed feedback has also been shown to provide a suitable explanatory mechanism for the coordinated pancreatic islet activity [Reference Bruce, Wei, Leng, Oh, Chiu, Roper and Bertram6].
Figure 1 (a) Diagrammatic overview of the glucose–insulin regulatory delayed-feedback model (2.1)–(2.2). (b)–(e) Characteristic time series of system (2.1)–(2.2) (with positive constant history) for different patterns of on–off glucose infusion. Intervals of fasting (infusion off) are indicated by a white background while intervals of glucose infusion (on) with constant rate
$G_{\text {max}}$
, period
$T_{in}$
and duration
$t_{in}$
are indicated by a light blue background. (b) and (e) Periodic oscillations; (c) damped oscillations; (d) irregular oscillations. Units are [G] mg dl
$^{-1}$
, [I] mU l
$^{-1}$
and [t] h. Infusion rates when not fasting are
$G_{\text {max}}=1.35$
mg dl
$^{-1}$
min
$^{-1}$
in panels (c)–(d) and
$G_{\text {max}}=24.3$
mg dl
$^{-1}$
min
$^{-1}$
in panel (e); period of infusion is
$T_{\text {in}}= 1$
h in panel (d) and
$T_{\text {in}}= 3$
h in panel (e); time of infusion is
$t_{\text {in}}= 30$
min in panel (d) and
$t_{\text {in}}= 5$
min in panel (e). See Section 2 for the mathematical formulation and implementation.
In this paper, we investigate a two-component, system-level mathematical model (see (2.1)–(2.2)) for blood glucose level
$G(t)$
and insulin level
$I(t)$
with two explicit time delays
$\tau _I$
and
$\tau _G$
corresponding to pancreatic insulin and hepatic glucose production pathways. The model incorporates the following physiological processes and factors that influence glucose and insulin dynamics, see Figure 1(a) for a schematic overview.
-
• Glucose uptake:
$G_{\text {in}}(t)$
represents the time-dependent glucose uptake into the blood by meal ingestion, continuous enteral or intravenous nutrition. -
• Insulin production:
$f_1$
represents the production of insulin. It is influenced by the concentration of glucose with a delay
$\tau _I$
to account for the time lag between high glucose levels triggering insulin production in the pancreas and when it becomes available for reducing glucose in the bloodstream. -
• Insulin-independent glucose utilization:
$f_2$
describes the utilization of glucose by tissues, mainly the brain, in an insulin-independent manner. It does not rely on the presence of insulin. -
• Insulin-dependent glucose utilization:
$f_3\cdot f_4$
represents the utilization of glucose by muscle tissues in an insulin-dependent manner. It reflects the capacity of tissues to use insulin for glucose uptake. -
• Glucose production by the liver:
$f_5$
represents the production of glucose by the liver. The delay
$\tau _G$
represents the time between hepatic glucose production and insulin stimulation. -
• Insulin degradation: The rate d accounts for the degradation of insulin in the body, primarily by the liver and kidneys. It combines both natural factors (for example, exercise) and artificial factors (for example, medication) that influence the rate of insulin degradation.
The nonlinear pathways
$f_1, f_2, f_3, f_4$
and
$f_5$
are represented using Hill-type functions (which are commonly used in biological modelling [Reference Gábor and Banga14]), depending on the glucose and/or insulin concentration in the blood stream at the current time point t and at two specific lag values
$t-\tau _I$
and
$t-\tau _G$
. The delays
$\tau _I$
and
$\tau _G$
are important physiological parameters encapsulating the responsitivity of the signalling and production pathways. They are assumed to be constant for the purpose of this article, although in practice, they can vary between individuals, as well as during the day and lifespan, and especially in the presence of diabetes. The model originates from the work of Sturis and collaborators who devised a model of glucose and insulin ultradian oscillations which were observed experimentally under various conditions [Reference Sturis, Polonsky, Mosekilde and Van Cauter48]. It has been extensively analyzed by various authors (including authors of this paper) in the case of constant rates of glucose infusion [Reference Huard, Bridgewater and Angelova22, Reference Huard, Easton and Angelova23, Reference Li and Kuang31, Reference Li, Kuang and Mason32]. We also remark here that the model belongs to a larger class of models incorporating delays to capture secretion processes [Reference Makroglou, Li and Kuang34, Reference Shi, Kuang, Makroglou, Mokshagundam and Li43].
We extend these earlier efforts on the analysis of the model by studying its response to periodic variations of the parameter
$G_{\mathrm {{in}}}$
, that is, periodic variations of glucose uptake, by using available software for numerical analysis of delay differential equations (see Appendix B for details on numerical implementation). In particular, we consider on–off infusion, a form of periodic infusion that is comparatively easy to implement in practice, where the rate of glucose infusion periodically switches between a positive constant value and zero. Figures 1(b)–(e) show prototypical examples for the response of system (2.1)–(2.2) for various types of glucose uptake during fasting in panel (b), glucose infusion with a (relatively high) constant rate in panel (c) and periodic on–off infusion in panels (d),(e). We first investigate the loss of ultradian oscillations under sufficiently strong constant infusion, see Figures 1(b),(c). We then aim to study the effects of different glucose infusion patterns
$G_{\mathrm {{in}}}$
on glucose homeostasis, in particular, the transiton from quasi-periodicity to entrainment accompanied by potentially large peak values of glucose, as shown in Figures 1(d),(e), a transition which—as we will show—is governed by a Torus (or Neimark–Sacker) bifurcation.
2. Glucose–insulin regulatory delayed-feedback model
We consider the system-level mathematical model
$$ \begin{align} G'(t) &= G_{\text{in}}(t) - f_2(G(t)) - f_3(G(t))f_4(I(t)) + f_5(I(t-\tau_G)),\end{align} $$
$$ \begin{align} I'(t) &= f_1(G(t-\tau_I)) - d I(t), \end{align} $$
with variables
$I(t)$
and
$G(t)$
representing the quantities of glucose (mg) and insulin (mU) in the plasma at time instant t. These are subsequently converted to concentrations in standard units (that is, mg dl
$^{-1}$
and mU l
$^{-1}$
for glucose and insulin, respectively) for all graphs in this paper by dividing by the corresponding compartmental volume (more precisely
$V_g$
and
$V_p$
, see below).
The model has been proposed and studied by Li and collaborators [Reference Li and Kuang31, Reference Li, Kuang and Mason32] as a model with explicit time delays to study the ultradian oscillations of insulin secretion, building on the original work of Sturis et al. [Reference Sturis, Van Cauter, Blackman and Polonsky49]. Specifically, system (2.1)–(2.2) explicitly depends on time delays
$\tau _I$
and
$\tau _G$
respectively representing the system’s response time to insulin production as a result to glucose uptake and the production of glucose by the liver as a result of low insulin levels. Glucose intake is modelled by parameters
$G_{\text {max}}$
,
$T_{\text {in}}$
and
$t_{\text {in}}$
, which respectively represent the maximal value of the infusion, the time between infusions and the duration of each infusion. The physiological response of the body is modelled by the nonlinearities
$$ \begin{align} \begin{aligned} f_1(G) &= \frac{R_m G^{h_1}}{G^{h_1} + (V_g k_1)^{h_1}},\quad f_2(G) = \frac{U_b G^{h_2}}{G^{h_2} + (V_g k_2)^{h_2}},\quad f_3(G) = \frac{G}{C_3 V_g},\\ \quad f_4(I) &= U_0 + \frac{ (U_m - U_0)I^{h_4}}{I^{h_4} + (1/V_i + 1/(Et_i))^{-h_4}k_4^{h_4}}, \quad f_5(I) = \frac{R_g I^{h_5}}{I^{h_5} + (V_p k_5)^{h_5}}, \end{aligned} \end{align} $$
where typical parameter values [Reference Huard, Easton and Angelova23, Reference Li, Kuang and Mason32] are provided in Table 1 with corresponding units.
Insulin degradation is modelled by a constant rate d. Throughout the paper, we fix
$d=0.06$
. The model has been considered before and has been analysed extensively including some of the authors [Reference Huard, Bridgewater and Angelova22, Reference Huard, Easton and Angelova23, Reference Li and Kuang31, Reference Li, Kuang and Mason32]. In particular, it can be shown that, for the parameter values considered and in the absence of infusion, there is a unique equilibrium solution
$(G^\ast ,I^\ast )$
[Reference Bennett and Gourley3, Reference Li and Kuang31]. We fix delay parameters used for numerical simulation to
$\tau _I=5$
and
$\tau _G=20$
if not stated otherwise. For the general theory of delay differential equations, such as existence, uniqueness and the stability of solutions, we refer the interested reader to classic textbooks on the topic [Reference Diekmann, Van Gils, Verduyn Lunel and Walther9, Reference Hale and Verduyn Lunel21].
Table 1 Values and units for parameters appearing in model functions (2.3), (from [Reference Huard, Easton and Angelova23, Reference Sturis, Van Cauter, Blackman and Polonsky49]).
For the purposes of numerical bifurcation analysis, we choose to model the (discontinuous) on–off periodic infusion by a smooth, periodic, quickly varying function between
$G_{\text {in}}(t)=0$
in mg (dl
$\cdot $
min)
$^{-1}$
(no infusion) and
$G_{\text {in}}(t)=G_{\text {max}}$
(constant glucose infusion), respectively. This continuity is also more compatible with the physiological nature of the system, wherein glucose levels can only vary continuously in response to a stimulus due to diffusion processes. We consider the specific form
$$ \begin{align} \begin{aligned} G_{\text{in}}(t) &= G_{\text{max}}\cdot s(t-\sigma_G), \\ s(t) &=h(\sin(2\pi t/T_{\text{in}}))\cdot h(\sin(2\pi(t-t_{\text{in}})/ T_{\text{in}} - \pi)), \end{aligned}\end{align} $$
where the sigmoidal function
$h(y)=(1+\exp (-k y))^{-1}$
can be considered as a smooth version of the Heaviside step function
$H(y)=0$
if
$y<0$
, and
$H(y)=1$
if
$y\geq 0$
for sufficiently large k; see Appendix B for details about the numerical implementation. The form of (2.4) was inspired by a model for auditory perception [Reference Ferrario and Rankin12]. Here,
$T_{\text {in}}$
is the time between consecutive shots of glucose/insulin with duration
$t_{\text {in}}$
, and the lag
$\sigma _G$
can be used to specify the timing of the infusion with respect to the underlying oscillation. The parameter k models the initial and terminal variations at the beginning and end of the shot application. We find
$k=100$
to be sufficiently large to generate a well-enough approximation h to the Heaviside function; in particular, we do not observe significant dynamical changes when choosing larger k. Figures 2(a) and (b) show the specific shape of the infusion patterns used for numerical computation of time series in Figure 1(d),(e).
Figure 2 Form of glucose infusion used to obtain Figures 1(d),(e). Parameters are: (a)
$G_{\text {max}}=1.35$
mg dl
$^{-1}$
min
$^{-1}$
,
$T_{\text {in}}= 1$
h and
$t_{\text {in}}= 30$
min; and (b)
$G_{\text {max}}=24.3$
mg dl
$^{-1}$
min
$^{-1}$
,
$T_{\text {in}}= 3$
h,
$t_{\text {in}}= 5$
min.
3. Results of numerical bifurcation analysis
3.1. Constant glucose infusion and ultradian oscillations
It has been shown that for a fixed constant glucose infusion
$G_{\mathrm {{in}}}$
, sufficiently large values of the response delays
$\tau _1$
and
$\tau _2$
lead to periodic oscillations in system (2.1)–(2.2) with periods in the experimentally observed range for ultradian oscillations [Reference Huard, Bridgewater and Angelova22, Reference Huard, Easton and Angelova23, Reference Li, Kuang and Mason32]. Mathematically speaking, the onset of oscillations is mediated by a supercritical Hopf bifurcation that leads to a local topological change in the solution space of system (2.1)–(2.2) from a stable equilibrium to a situation of an unstable equilibrium surrounded by a small stable limit cycle (to the right of the critical curve, that is, for larger values of the delays) close to the bifurcation point [Reference Li and Kuang31]. For details on bifurcation theory and the Hopf bifurcation, we refer the interested reader to [Reference Kuznetsov30]. To witness the bifurcation point, it is necessary to vary at least one parameter of the system. Here, we focus on the response delays
$\tau _I$
and
$\tau _G$
. Allowing these two parameters to vary simultaneously, one obtains a weakly nonlinear one-parameter curve
$\mathbf {H(\omega )}=(\tau _I(\omega ),\tau _G(\omega ))$
of Hopf bifurcation in the
$(\tau _I,\tau _G)$
-plane in terms of the Hopf frequencies
$\omega $
, (see Appendix A for a detailed derivation). The curve
$\mathbf {H(\omega )}$
corresponds to the critical curve for oscillations in system (2.1)–(2.2).
Figure 3(a) shows the curve
$\mathbf {H}$
(black) during fasting, that is, for
$G_{\text {in}(t)}=0$
, computed with the software package DDE-Biftool for Matlab [Reference Engelborghs, Lemaire, Belair and Roose10, Reference Engelborghs, Luzyanina and Roose11, Reference Sieber, Engelborghs, Luzyanina, Samaey and Roose44]. It has been numerically verified that the curve
$\mathbf {H}$
is indeed supercritical for the range of parameter values considered. Figure 3(a) can be interpreted as follows. First, for value pairs above the curve and for
$\tau _I\leq 20~\mathrm \min , \tau _G\leq 60~\mathrm \min $
, any solution of the model starting in a physiological range of glucose and insulin develops periodic oscillations (see Figure 1(b)). Second, for value pairs
$(\tau _I,\tau _G)$
below the curve
$\mathbf {H}$
, oscillations in system (2.1)–(2.2) decay and approach the equilibrium
$(G^\ast ,I^\ast ).$
This corresponds to a situation where the delay coupling is not strong enough to deviate from the equilibrium and produce detectable oscillations. In the nonfasting case, we also note that the constant administration of glucose dose that is too high to be managed in an oscillatory manner within physiological glucose and insulin ranges also leads to a loss of oscillations, compare Figure 1(c).
Figure 3 Characterization of fasting oscillations with respect to response delays. Panels show the (a) period, (b) maximum glucose value and (c) minimum glucose value as a function of response delays
$\tau _I$
(min) and
$\tau _G$
(min). Shown are the critical curve for oscillations (black, Hopf bifurcation) and iso-curves (blue) with constant period in panel (a), glucose-maxima in panel (b) and minimum of G in panel (c). The light blue rectangle indicates the physiologically amenable range of delay values for comparison. See Section 2 for the model and choice of parameters.
Figure 4 Position of the critical curve (curve of Hopf bifurcation, supercritical) in the
$(\tau _I,\tau _G)$
-plane for various values of constant glucose infusion
$G_{\text {in}}=G_{\text {max}}$
ranging from (a)
$0$
to
$0.5$
and from (b)
$0.6$
to
$1.6$
mg dl
$^{-1}$
min
$^{-1}$
(all black). The light blue rectangle shows the physiological range of delay values for comparison.
Figure 3(a) also gives an overview of the resulting period of oscillation above the critical curve
$\mathbf {H}$
shown in the form of isocurves (blue) of limit cycles with constant period. The physiological range of parameters
$(\tau _I,\tau _G)$
is highlighted by a light blue square in the background for convenience. The range of expected periods for ultradian oscillations as predicted by the model thus ranges from
$2.2$
to
$4.2$
hours during fasting. More generally, we observe the period of the limit cycle oscillation grows approximately linear with the sum of the two delay values
$\tau _I+\tau _G$
. We also observe that away from the curve
$\mathbf {H}$
, the limit cycle oscillation becomes increasingly less sinusoidal, that is, the nonlinearity of system (2.1)–(2.2) has increasingly more of an effect on the limit cycle. Figures 3(b),(c) illustrate this effect by plotting isocurves of periodic orbit with constant minimum and maximum glucose within one period of oscillation. We observe that, whereas the glucose minimum decreases approximately linearly with the sum of the delays
$\tau _I+\tau _G$
, the maximum G remains almost constant for the range of parameter values considered. Note that this predicted effect of long response delays is potentially harmful and is virtually undetectable by common testing methods.
However, we observe that, for fixed values of the delays, gradually increasing the glucose infusion leads to a loss of oscillations. This phenomenon has been observed before and can be interpreted as an insufficient insulin secretion to accommodate the infusion, forcing the system to lower glucose levels [Reference Li and Kuang31]. Figure 4 shows how the location of the curve
$\mathbf {H}$
changes for various levels of constant glucose infusion,
$G_{\text {in}}(t) = G_{\text {max}}$
. We observe two different types of change for values in the approximate ranges
$0\leq G_{\text {max}}\leq 0.55$
mg dl
$^{-1}$
min
$^{-1}$
and
$G_{\text {max}}>0.55$
mg dl
$^{-1}$
min
$^{-1}$
shown in Figures 4(a) and (b), respectively. Figure 4(a) suggests that low levels of
$G_{\text {in}}$
promote oscillations in system (2.1)–(2.2) as compared with the fasting case. This trend reverses at approximately
$G_{\text {max}}=0.55$
mg dl
$^{-1}$
min
$^{-1}$
, where the location of the curve
$\mathbf {H}$
starts moving to increasingly larger values of
$\tau _I$
and
$\tau _G$
, see Figure 4(b). Approximately at
$G_{\text {max}}=1.2$
mg dl
$^{-1}$
min
$^{-1}$
, the position of
$\mathbf {H}$
is comparable with the starting location for
$G_{\text {max}}=0$
. Further increasing
$G_{\text {max}}$
moves
$\mathbf {H}$
inside the physiological range of delay values (light blue) and finally beyond causing all oscillations to cease in the physiological parameter regime. Compare also Figures 1(b),(c) for an illustration of this transition and the loss of oscillations for
$(\tau _I,\tau _G)=(5,20)$
.
3.2. Amplitude response to on–off glucose infusion
We now investigate the effect of periodic glucose infusion on baseline fasting oscillations shown in Figure 1(b), that is, we fix
$\tau _I=5$
min,
$\tau _G=20$
min and periodically adjust the level of
$G_{\text {in}}$
between
$0$
and a positive value
$G_{\text {max}}$
to be specified. The natural frequency of ultradian oscillation in this case is
$T_0 \approx 2.2$
h. We show that the resulting glucose and insulin ranges depend sensitively on the period of the on–off infusion.
3.2.1. Long infusion time compared with period
Figures 1(d),(e) show two of the possible outcomes with different maximal infusion strength
$G_{\text {max}},$
period of infusion
$T_{\text {in}}$
and infusion duration
$t_{\text {in}}.$
Figure 1(d) shows the result of periodic infusion with
$G_{\text {in}}=1.35$
mg dl
$^{-1}$
min
$^{-1}$
for
$t_{\text {in}}=30$
min every
$T_{\text {in}}=60$
min, resulting in so-called quasi-periodic oscillations. Indeed, quasi-periodic oscillations are characterized by the presence of an oscillating envelope of the oscillation that evolves on a much slower time-scale. This is in sharp contrast with panels (b) (no infusion) and (c) (constant infusion with the same maximum rate) of Figure 1, where we have either periodic oscillations or a decay of oscillations towards the equilibrium state. Quasi-periodic oscillations can be expected to occur in oscillatory systems which are externally driven by an input with noncommensurable period, here
$T_{\text {in}}/T_0=2.2$
. In this case, the effect of infusion very much depends on its current state: when insulin is low, glucose increases quickly; when insulin is high, glucose cannot increase further and the infusion only delays the expected decrease in glucose levels.
Periodicity of the oscillations can be restored by adjusting
$G_{\text {max}}$
and
$T_{\text {in}}$
. Figure 5 summarizes the response of system (2.1)–(2.2) to periodic forcing with different values of
$G_{\text {max}}$
and
$T_{\text {in}}$
. The locus in parameter space of the quasi-periodic oscillation shown in Figure 1(d) is indicated by a green rectangle. Figure 5 shows the overall glucose maximum (in colour code) observed over a time span of
$100\cdot (T_{\text {in}}+\tau _I+\tau _G)$
minutes. The various mechanisms generating periodic rhythms can be understood from the numerically computed bifurcation curves shown in Figure 5. These correspond to curves of torus bifurcations
$\mathbf {T}$
(purple), curves
$\mathbf {F}$
(red) of fold (or saddle-node) bifurcations of periodic orbits and curves
$\mathbf {PD}$
(magenta) of period-doubling bifurcations of periodic orbits. All these curves mark the transition to periodic solutions and thus characterize the so-called entrainment of oscillations.
Figure 5 Response of model (2.1)–(2.2) to glucose infusion protocol (2.4) with maximum infusion rate
$G_{\text {max}}$
(mg (dl min)
$^{-1}$
) and length of infusion
$t_{\text {in}}=T_{\text {in}}/2$
(h). Shown is the maximum value of G (mg dl
$^{-1}$
) in colourcode (blue–white) obtained by integration for various
$T_{\text {in}}$
and
$G_{\text {max}}$
over
$100(\tau _I+\tau _G+T_{\text {in}})$
time units. The maximum data are overlaid by the curve of torus bifurcation (purple), curves of fold bifurcation of periodic orbits (red) and curves of period doubling bifurcation (magenta) bounding regions of locking to the infusion protocol. Other parameters are
$\tau _I=5$
min and
$\tau _G=20$
min. The green square indicates the parameter values leading to the quasi-periodic behaviour plotted in figure 1(d).
The curves
$\mathbf {F}$
respectively enclose deltoid-like regions—which can be viewed as resonance or locking tongues—extending from the line
$G_{\text {max}}=0$
, inside of which we observe periodic oscillations. The curves
$\mathbf {F}$
emerge pairwise from resonant points where the infusion period is a rational multiple of the natural period of the system without infusion, that is,
$pT_{\text {in}}=qT_0$
for integers
$p,q$
. Figure 5 shows the first three principal resonances of system (2.1)–(2.2), where
$p=1,2,3$
and
$q=1$
. It is expected that such resonance tongues emanate from the line
$G_{\text {max}}=0$
at every rational value of
$T_0/T_{\text {in}}$
. These higher order resonances (except
$p=4$
and
$q=1$
which is outside of the considered range of parameter values) have been omitted from the computation as they are typically very narrow and thus unlikely to be of physiological relevance.
This behaviour persists moving towards larger values of
$G_{\text {max}}$
into the regions that are bounded approximately by the curves T, where the underlying stable periodic orbit destabilizes and gives rise to a torus that corresponds to quasi-periodic oscillations. We find numerical evidence that the direction with which this torus emanates from the curve
$\mathbf {T}$
can change and gives rise to the discontinuous transition between the observed maximum values in Figure 5. The curves
$\mathbf {T}$
each emanate from either point of intersection with a curve
$\mathbf {F}$
or
$\mathbf {PD}$
. Intersections with curves
$\mathbf {PD}$
correspond to higher order locking between the ultradian oscillations and the infusion. Overall, we observe that the strength and period of the infusion have a crucial effect on the resulting amplitude of the oscillations. For instance, forcing the system periodically with
$T_0=T_{\text {in}}$
and relative amplitude
$G_{\text {max}}=1$
mg dl
$^{-1}$
min
$^{-1}$
leads to a
$40\%$
increase of the overall amplitude of the oscillation (which appears to be still in the physiological range). In contrast, stimulating the system with a gradually increasing
$G_{\text {max}}$
in the 2:1 regime first goes through a phase during which glucose amplitudes remain relatively constant before slowly increasing.
More generally, we observe that, for the assumed values of the response delays, periodic infusion with
$T_{\text {in}}=2t_{\text {in}}>T_0$
and
$G_{\text {max}}$
is sufficient for the resulting period of the resulting glucose–insulin oscillation to be set by (locked to) the period of glucose infusion.
3.2.2. Short infusion time compared with period
We note here that locking can be achieved when the same glucose dose is delivered in a shorter period of time, resulting in a more concentrated and intense infusion. To further explore this phenomenon, we conducted additional experiments using an on–off glucose infusion protocol with a fixed infusion period of
$T_{\mathrm {{in}}}=180$
min. Figure 6 showcases the results obtained from these experiments, where we varied both the infusion time
$t_{\mathrm {{in}}}$
and the average glucose dose per minute
$\bar G$
, represented by
$G_{\text {max}}\cdot t_{\mathrm {{in}}} / T_{\mathrm {{in}}}$
. We observe a locus in the parameter space that corresponds to the quasi-periodic orbit illustrated in Figure 1(e). This locus is denoted by a distinctive yellow diamond marker, which highlights the specific combination of infusion time and glucose dose that leads to the observed quasi-periodic behaviour. Additionally, we present a curve labelled as
$\mathbf {T}$
, which represents a torus bifurcation curve. This curve serves as an indicator of the critical transition point between entrainment and quasi-periodic oscillation in response to the infusion protocol.
Figure 6 Response of model (2.1)–(2.2) to glucose infusion protocol (2.4) with average infusion rate
$\bar G=G_{\text {max}} \cdot t_{\text {in}}/T_{\text {in}}$
mg dl
$^{-1}$
min
$^{-1}$
over the length of infusion
$t_{\text {in}}$
(min) with constant period
${T_{\text {in}}=180}$
(min). Shown is the maximum value of G (mg dl
$^{-1}$
) in colourcode (blue–red) obtained by integration for various
$t_{\text {in}}$
and
$ \bar G$
. The maximum data are overlaid by the torus bifurcation curve (purple). Other parameters are
$\tau _I=5$
min and
$\tau _G=20$
min.
4. Discussion
It is well documented that glucose rhythms stimulate pulsatile pancreatic insulin secretion at various timescales [Reference Satin, Butler, Ha and Sherman41, Reference Sturis, Polonsky, Mosekilde and Van Cauter48]. For example, the 1:1 entrainment mode—namely one ultradian glucose oscillation per glucose infusion cycle—was clinically shown to be present using a sinusoidal glucose infusion in individuals without diabetes [Reference O’Meara, Sturis, Van Cauter and Polonsky38, Reference Sturis, Knudsen, O’Meara, Thomsen, Mosekilde, Van Cauter and Polonsky47]. Our analysis of periodically driven ultradian oscillations highlights that a periodic on–off stimulus, closer to normal daily conditions, also possesses the ability to entrain glucose rhythms. Furthermore, the duration of each glucose input has a crucial impact on the generation of periodic rhythms, as well as on attained glycemic levels. This theoretically provides a method for delivering a fixed glucose dose while minimizing the amplitude of the resulting rhythm. This can be achieved by either altering the period of the infusion or the length of each pulse. This is most observable in Figure 6, where stretching the infusion duration leads to lower glucose amplitudes. For example, consider a scenario where glucose is infused every 180 minute over a 12-hour period. Infusing a dose with
$G_{\text {max}} = 2.4$
mg dl
$^{-1}$
min
$^{-1}$
over
$t_{in}=30$
minutes leads to a maximal glucose value around
$150$
mg dl
$^{-1}$
. In contrast, a dose with
$G_{\text {max}} = 1.2$
mg dl
$^{-1}$
min
$^{-1}$
over
$t_{in}=60$
minutes reduces the maximal glucose level to approximately
$125$
mg dl
$^{-1}$
. In both cases, the average dose per minute is
$\bar {G} = 0.4$
mg dl
$^{-1}$
min
$^{-1}$
, and a total dose of
$288$
mg dl
$^{-1}$
is infused over the 12-hour timespan.
Our study shows that the system’s response to glucose infusion patterns provides multiple pathways for the production of stable oscillatory rhythms and entrainment, which has also been shown for simpler models of glucose–insulin regulation featuring delays, for example, [Reference Li, Wang, De Gaetano, Palumbo and Panunzi33, Reference Panunzi, Palumbo and De Gaetano40, Reference Shi, Kuang, Makroglou, Mokshagundam and Li43]. However, it is clear from Figures 5–6 that only measuring maximum glucose levels is not sufficient to characterize the response of the system to periodic on–off infusion, and a more or less continuous-time measurement of, at least, glucose is required. There are several other limitations that should be mentioned here. First, let us note that while the exact location of bifurcation curves would depend on model parameters, the bifurcation types are likely to remain the same for parameter ranges representing nondiabetic individuals. Our model assumes fixed values for the delays in insulin and glucose production pathways, represented by
$\tau _I$
and
$\tau _G$
, respectively. In reality, these delays can vary between individuals and change over short and long timescales due to daily-life factors such as exercise, aging and the presence of insulin resistance. Future research could incorporate individual-specific delays to account for this variability and investigate their impact on the system’s dynamics.
It is worth noting that our model relies solely on plasma glucose and insulin measurements for prediction, which highlights the importance of accurate and reliable measurements in clinical settings. The nonlinear structure of the model allows for the description of nontrivial dynamics and enhances parameter identifiability. This aspect is crucial for developing robust and accurate models that can capture the complex dynamics of the glucose–insulin regulatory system.
Moreover, the timing of the glucose infusion does not influence the bifurcation structure (5), nor the glucose–insulin ranges of the periodic rhythms. In other words, the long-term dynamics is not dependent on the starting time of the periodic on–off glucose infusion. This is not to say that timing does not have a crucial importance. While the investigated infusion ranges ensured the positivity of glucose and insulin values, values below or above healthy physiological ranges may appear in the transient path to the limit cycle. We emphasize that we exclusively consider here a continuous glucose stimulus in this paper and that, in this context, no evidence of delay-induced uncertainty could be observed [Reference Karamched, Hripcsak, Albers and Ott25, Reference Karamched, Hripcsak, Leibel, Albers and Ott26].
In turn, additional dynamics may emerge from interactions with other physiological feedback loops or subsystems, such as the glucagon pathway or the hypothalamic-pituitary-adrenal axis, for which the alignment with glucose regulation is essential for maintaining good health [Reference Grant, Lewis and Kriegsfeld18, Reference Zavala55]. The recent incorporation of glucagon [Reference Cohen and Li8] in models of the glucose–insulin feedback system may help provide a more complete and quantitative picture of dynamical interactions occurring within the pancreas [Reference Gram Pedersen, Mosekilde, Polonsky and Luciani17, Reference Montefusco, Cortese and Pedersen36], which can be used to improve quantitative tests for the detection and measurement of insulin and glucagon resistance [Reference Morettini, Burattini, Göbl, Pacini, Ahrén and Tura37].
Another aspect to consider is the interaction between the glucose–insulin regulatory system and other physiological processes. Our model focuses solely on the glucose–insulin loop, but in reality, there are complex interactions between various metabolic pathways, hormones and organs. Integrating these interactions into a comprehensive model could provide a more complete understanding of the system’s behaviour and its response to different stimuli.
5. Conclusion
In this study, we employed a system-level mathematical model to investigate the response of the glucose–insulin regulatory system to periodic glucose infusion. By exploring different glucose infusion patterns and analyzing the resulting dynamics, we gained insights into the system’s behaviour and identified key factors (such as the infusion amplitude and period) influencing its response.
Our findings demonstrate that the glucose–insulin regulatory system exhibits a range of behaviours depending on the glucose infusion pattern. When a constant glucose infusion is applied, the system shows ultradian oscillations characterized by periodic variations in glucose and insulin levels. However, as the glucose infusion rate exceeds a certain threshold, these oscillations disappear and the system focuses on reducing glucose levels without exhibiting oscillatory behaviour. This observation suggests a physiological limit beyond which the system’s oscillatory capacity is overwhelmed.
We further investigated the effects of periodic on–off pulses, mimicking repeated intravenous glucose tolerance tests. Our analysis revealed that the period of the on–off pulses plays a crucial role in determining the system’s dynamics and glucose–insulin ranges. Different patterns of oscillations, including stable limit cycles and irregular oscillations, were observed for varying infusion periods. This highlights the importance of considering the frequency and duration of glucose stimuli in understanding the system’s response.
We identified the impact of different glucose infusion patterns on the system’s dynamics and demonstrated the importance of various types of glucose stimuli. These insights can potentially aid in the development of diagnostic and therapeutic strategies for glucose regulation and motivate new strategies in the management of metabolic disorders. Future research should aim to incorporate individual-specific delays, and consider the broader physiological context to further refine our understanding of glucose regulation and its implications for human health.
Appendix A Critical delay values for oscillatory behaviour when infusion rate is constant
The critical curve for oscillations in the (
$\tau _I,\tau _G$
)-parameter plane can be computed from the linearization of system (2.1)–(2.2) about the equilibrium solution
$(G^\ast ,I^\ast )$
and imposing the condition
$\lambda =i\omega ,~\omega>0$
(Hopf bifurcation) on solutions of the corresponding characteristic equation
$$ \begin{align} 0=\chi(\lambda):=\lambda^{2}+{\alpha}_{1} \lambda +{\alpha}_{0}+{\beta}_{1} e^{-\lambda {\tau}_{1}} + \beta_{2} e^{-\lambda {\tau}_{2}}, \end{align} $$
where
$\tau _1=\tau _I, \tau _2=\tau _I+\tau _G$
and
$\alpha _1=f_2^{\prime } (G^\ast ) + f^{\prime }_3(G^\ast )f_4(I^\ast )+d$
,
$\alpha _0=d(f_2^{\prime } (G^\ast ) + f^{\prime }_3(G^\ast )f_4(I^\ast ))$
,
$\beta _1=f^{\prime }_ 1(G^\ast )f_3(G^\ast )f^{\prime }_4 (I^\ast )$
,
$\beta _2=-f^{\prime }_1 (G^\ast )f^{\prime }_5(I^\ast )$
. A detailed derivation of (A.1) can be found in [Reference Huard, Bridgewater and Angelova22].
The equation
$0=\chi (i\omega )$
can be solved parametrically for
$\tau _1$
and
$\tau _2$
to give
$$ \begin{align}\tau_{1,2}(\omega)&=\frac{1}{\omega}\bigg\{\arctan\bigg(\frac{\alpha_1\omega}{\omega^2-\alpha_0}\bigg)\nonumber\\& \quad+ \arccos\bigg(\frac{\beta_{2,1}^2-\beta_{1,2}^2-(\omega^2-\alpha_0)^2 -\alpha_1^2\omega^2}{2\beta_{1,2}\sqrt{(\omega^2-\alpha_0)^2 +\alpha_1^2\omega^2}}\bigg)\bigg\} \end{align} $$
revealing the critical curve for oscillations
$\mathbf {H}\subset \mathbb {R}^2$
(curve of Hopf bifurcation)
$$ \begin{align} \mathbf{H}(\omega)=(\tau_I(\omega),\tau_G(\omega))=(\tau_1(\omega),\tau_2(\omega)-\tau_1(\omega)). \end{align} $$
For the considered parameter values, we have that
$\alpha _1> \alpha _0$
and
$\beta _2>\alpha _0$
, ensuring the existence of
$\mathbf {H}$
. Indeed, the curve is a sharp threshold for oscillation, as it can been shown numerically that for positive values (
$\tau _I,\tau _G$
) below
$\mathbf {H}$
, the fixed point
$(G^\ast ,I^\ast )$
is stable for any physiological range of starting values G and I. It is worth noting here that system (2.1)–(2.2) undergoes further Hopf bifurcations, respectively at
$$ \begin{align*}\tau_{I,k}(\omega)=\tau_{I}(\omega)+2\pi k/\omega, \quad \tau_{G,l}(\omega)=\tau_{G}(\omega)+2\pi l/\omega\end{align*} $$
with integers
$k,l$
; however, for the parameter values considered, we can restrict ourselves to the smallest positive such value pair to cover the physiological parameter range. The range of relevant values of
$\omega $
resulting in positive delays cannot be computed explicitly; however, straightforward calculations show that the boundaries
$\omega _{I},\omega _{G}$
satisfying
$\tau _I(\omega _{I})=0$
and
$\tau _G(\omega _{G})=0$
are given by
$$ \begin{align*} \omega_G &= \sqrt{\alpha_0-\frac{\alpha_1^2}{2} +\sqrt{\bigg(\alpha_0-\frac{\alpha_1^2}{2}\bigg)^2+(\beta_1+\beta_2)^2-\alpha_0^2}},\\ \omega_I &= \sqrt{\alpha_0+\beta_1-\frac{\alpha_1^2}{2} +\sqrt{\bigg(\alpha_0+\beta_1-\frac{\alpha_1^2}{2}\bigg)^2+\beta_2^2-\alpha_0^2}}, \end{align*} $$
with the corresponding delay values
$$ \begin{align*} \tau_I(\omega_{G}) &=\frac{1}{\omega_{G}}\arctan \bigg(\frac{\alpha_1\omega_{G}}{\omega_{G}^2-\alpha_0}\bigg) + \frac{2\pi k^\ast}{\omega_G},\\ \tau_G(\omega_{I}) &=\frac{1}{\omega_{I}}\arctan \bigg(\frac{\alpha_1\omega_{I}}{\omega_{I}^2-\alpha_0-\beta_1}\bigg)+ \frac{2\pi l^\ast}{\omega_I}, \end{align*} $$
where
$k^\ast ,l^\ast $
are the smallest integers such that
$\tau _G$
and
$\tau _I$
are positive.
We remark that the curve
$\mathbf {H}$
vaguely resembles a straight line with slope
$-1$
in the
$(\tau _I,\tau _G)$
-plane. This can be understood by exploiting the fact that parameter
$\beta _1$
is small of the order of
$10^{-3}$
. Imposing the regular perturbation ansatz
$\omega =\omega _0+\beta _1\omega _1+\mathcal {O}(\beta _1^2)$
on the imaginary part of (A.1) and comparing at zeroth and first order in
$\beta _1$
, we formally obtain
$$ \begin{align*} \omega_0 &= \sqrt{\alpha_0-\frac{\alpha_1^2}{2} +\sqrt{\bigg(\alpha_0-\frac{\alpha_1^2}{2}\bigg)^2+\beta_2^2-\alpha_0^2}},\\ \omega_1 &= \frac{\omega_0}{\alpha_1-\alpha_0+\omega_0^2}\tau_I \leq \frac{1}{2 \sqrt{\alpha_1-\alpha_0}} \tau_I. \end{align*} $$
Thus, we can approximate
$\mathbf {H}$
to first order in
$\beta _1$
$$ \begin{align*} \mathbf{H}(\omega_0 + \beta_1\omega_1(\tau_I))&\approx(\tau_I,\tau_2(\omega_0 + \beta_1\omega_1(\tau_I))-\tau_I) \end{align*} $$
by using the expression
$\tau _2(\omega )$
in (A.2). As a result,
$\mathbf {H}$
approaches the graph of the function
$\tau _I\mapsto \tau _2(\omega _0)-\tau _I$
with slope
$-1$
as
$\beta _1\to 0$
, which can be considered as a zero-order approximation of
$\mathbf {H}$
.
Appendix B Numerical bifurcation analysis of time-delay systems with periodic infusion
Numerical simulations have been obtained using pydelay [Reference Flunkert and Schoell13]. Numerical bifurcation analysis has been performed using the software package DDE-BIFTOOL for Matlab/Octave [Reference Sieber, Engelborghs, Luzyanina, Samaey and Roose44]. For a general introduction to numerical continuation methods available for delay differential equations and their application to physiological systems, see [Reference Engelborghs, Luzyanina and Roose11, Reference Krauskopf and Sieber29], respectively. Isocurves in Figure 3 have been computed using numerical continuation of periodic orbits in two parameters with the additional condition fixed period in panel (a), and fixed maximum value in panel (b) and fixed minimum value in panel (c), where in cases (b) and (c), we also relaxed the phase condition. For bifurcation analysis in the presence of periodic infusion, we append the two-dimensional ordinary differential equation
$$ \begin{align*} x^\prime (t) &= x - \omega y(t) - x(t)(x(t)^2+y(t)^2),\\ y^\prime (t) &= -\omega x(t) + y - y(t)(x(t)^2+y(t)^2), \end{align*} $$
with known stable periodic solution
$(x(t),y(t))=(\cos (\omega t),\sin (\omega t))$
to system (2.1)–(2.2). The method has been employed in several other works, (see, for example, [Reference Keane, Krauskopf and Postlethwaite27]). We achieve the specific form of infusion (2.4) by setting
$$ \begin{align*}G_{\text{in}}(t) = G_{\text{max}} h(y(t-\sigma_G))h\bigg(y\bigg(t-\sigma_G-t_{in}-\frac{\pi}{\omega}\bigg)\bigg),\end{align*} $$
where
$\omega =2\pi /T_{\text {in}}$
.
Acknowledgements
S. Ruschel was supported by UKRI Grant No. EP/Y027531/1. The authors thank Jan Sieber for helpful discussions on the implementation of the numerical continuation methods in DDE-BIFTOOL.
Open access
