Analysis of a stochastic IVGTT glucose-insulin model with time delay

1.   Introduction

Diabetes mellituse has been one of the major diseases in the world due to the high percentage of diabetics in the whole population and the increasing growth rate of its onset. It is estimated that 415 million people are living with diabetes in the world, or 1 in 11 of the world adult population. In addition, 46% of people with diabetes are not diagnosed yet. The number of diabetics is expected to rise to 642 million by 2040 (https://www.diabetes.co.uk/diabetes-prevalence.html). Insulin sensitivity (IS) and glucose effectiveness (GE) are two most important physiological characteristic factors, which are not only used to assess the onset of type 2 diabetes mellituse (T2DM) in research laboratories to investigate the pathways to T2DM, but also evaluating the effectiveness of new drugs that increase the insulin sensitivity and/or glucose effectiveness in research laboratories and pharmaceutical manufacturers. To determine one's insulin sensitivity and glucose effectiveness, Glucose Clamp Test developed by DeFronzo et al. ^[1] is the gold standard for quantifying insulin sensitivity and further glucose effectiveness. However, the rigorousness for performing a clamp test and the sufferings of the subjects during and after the test enfeeble its applications. To mitigate the experiment process, intravenous glucose tolerance test (IVGTT) was developed and became a popular protocol used by research laboratories to estimate the insulin sensitivity and glucose effectiveness, in which, experimental data is fitted with a differential equation model and the resulting model parameters are used to quantify the ind ices of IS and GE. The most well known differential equation model is minimal model developed by Bergman, Cobelli and their colleagues ^[2,3], which was latter pointed out by De Gaetano and Arino ^[4] that the model is not well-post and fitting would fail in certain cases. This implies that suitable models and better parameter estimation methods are in demand (see ^{[4,5,6,7,8,9,10,11,12]}).

In real life, approximately forty two factors affect the glucose-insulin metabolic system, such as carbohydrates, medication, intense of activities, environments, and behaviors (https://diatribe.org/42factors). Many of these factors, for example, stress level and hormone cycle, are stochastic behavior. Stochastic differential equation (SDE) models are known as a powerful tool, which not only accounts for the white noise in a system of differential equations, but also is able to predict the future dynamics based on the corresponding deterministic system. There has been growing interest in utilizing SDE to model the blood sugar levels. For instance, Zhang et al. ^[13] developed a stochastic nonlinear second order differential equation with a Bayesian learning scheme to describe the blood sugar levels. Duun-Henriksen et al. applied SDE to improve the prediction and uncorrelated errors in a glucoregulatory system for type 1 diabetes mellitus ^[14]. While most studies introduced stochastic noise in modeling with the assumption that the noise process is stationary Gaussian process, interestingly, Beny$\acute{o}$ et al. ^[15] found that this stochastic term is Gaussian process but not stationary.

In the modeling of biological systems, time delay is often an important factor to be considered. The types of time delay can be divided into discrete delay, distributed delay, and internal (i.e. in the state) or external (i.e. in the input or output), see ^[16]. Discrete delay and distributed delay are found in many models, and we won't list them here. Particllarly, De la Sen ^[17] proposed a time-delay systems with non-commensurate internal point delays and discussed it's absolute stability. De la Sen ^[16] discussed the positivity properties of singular regular linear time-delay time-invariant systems subject to multiple internal and external incommensurate constant point delays. In the endocrine regulation system of blood glucose and insulin, there is a time delay for pancreatic cells to secrete insulin according to the change of blood glucose concentration. Many glucose-insulin models represent this delay in terms of discrete delays, see ^[10,18,19]. In 2017, Shi et al. ^[12] proposed a novel approach to model the insulin secretion time delay on the basis of model (1.1), which used two parameters to simulate both discrete time delay and distributed time delay.

In this paper, we propose a novel stochastic IVGTT model with a discrete time delay on basis of the work in ^[10]. We show that our model permits a unique global positive solution. Besides, we prove that although the proposed model doesn't have a positive equilibrium point, its solution perturbs around a point near the positive equilibrium point of the corresponding deterministic system without stochastic components. Moreover, we develop a maximum likelihood estimation (MLE) method to estimate the parameters involved in the model. Our numerical studies show that the proposed model with suitable parameters provided by the MLE method could improve the fitness of real data as compared to the corresponding deterministic system.

The paper is organized as follows. Section 2 presents our proposed stochastic IVGTT model with a time delay. In section 3, we prove the existence and uniqueness of the global positive solution of the proposed model and examine the asymptotic behavior of the solution. Moreover, We also discuss how the proposed model can be identified from the data. In section 4, we present a series of numerical simulations to verify our theoretical findings in the paper and reveal some intriguing dynamics of our stochastic model with respect to different noise disturbances. A brief discussion of the implications of our results is presented in the last section.

2.   A stochastic IVGTT model with a time delay

In IVGTT, after an overnight fasting so that both glucose and insulin remain to be at their vasal levels, the subject is injected with a bolus of glucose (300 mL/kg) into a subject's vein and then the blood is immediately sampled at the time mark 2', 4', 6', 8', 10', 12', 15', 20', 25', 30', 35', 40', 50', 60', 80', 100', 120', 140', 160' and 180' to measure the glucose and insulin concentrations. Among the aforementioned literature, Shi et al. ^[12] uses two parameters to simulate the insulin secretion delay and investigated impact of delay interval in the past. Li et al. ^[10] proposed a deterministic dynamic model and an approach to identify the length of time delay, which can be applied to reduce the number of parameters in fitting data.

The deterministic dynamic model studied by ^[10] is given as follows,

with initial condition $G(\theta) = \phi(\theta) > 0$ and $I(\theta) = \psi(\theta) > 0$ for $\theta\in[-\tau, 0]$, where $G(t) > 0$ and $I(t) > 0$ denote the glucose and insulin concentrations at time $t$, respectively; $b > 0$ is the rate constant of the hepatic glucose production; $S_g > 0$ is the consumption rate of the non-insulin-dependent glucose, also known as glucose effectiveness index; $S_i > 0$ is the consumption rate of the insulin-dependent glucose per unit of insulin concentration, also known as insulin sensitivity index; $\sigma > 0$ is the maximum insulin release rate; $d_i > 0$ is the constant insulin degradation rate, and $\sigma f(G(t-\tau))$ represents the insulin secretion response to glucose stimulation with time delay $\tau > 0$. According to physiology ^{[10,20,21,22]}, $f(\cdot)$ is in sigmoidal shape and takes the form as in ^[10]:

with $a > 0$ as the half$-$saturation and $r \geq 2$. Our novel stochastic model is built upon the system (2.1) proposed by Li et al. ^[10], in which, we only consider the stochastic disturbance on the glucose effectiveness index $S_g$ and insulin degradation rate $d_i$ for the sake of simplicity. More specifically, at a given time $t$, the glucose effectiveness index $S_g$ and insulin degradation rate $d_i$ are perturbed as

respectively, where $B_i(t), \; i = 1, 2,$ are two independent standard Brownian motions, and $\alpha_i (i = 1, 2)$ are positive constants with $\alpha_i^2$ representing the intensity of the white noise $B_i(t)$. As a result, we formulate following stochastic system

with the same initial condition as that in the system (2.1)

and the parameters $b, S_g, S_i, \sigma$ and $d_i$ have the same meaning as the system (2.1).

We shall investigate the influences of random disturbance on system dynamics, and use the proposed model to fit the IVGTT data sets.

3.   Theoretical properties

In this section, we study the theoretical properties of the proposed stochastic system (2.2), which show that the system (2.2) is well-post. We start with introducing some necessary notations, definition and lemmas which will be used to prove our main results.

Let $(\Omega, \mathcal{F}, \{\mathcal{F}_t\}_{t\geq0}, P)$ be a complete probability space with a filtration $\{\mathcal{F}_t\}_{t\geq0}$ satisfying the usual conditions that $\mathcal{F}_0$ contains all $P$ null sets and $\mathcal{F}_{t_{1}}\subseteq \mathcal{F}_{t_{2}}$ if $t_1\leq t_2$. Let $B_i(t), i = 1, 2$ denote two scalar Brownian motions defined on the complete probability space $\Omega$. Also let $Y(t)$ = $(G(t), I(t))^{ T}$ and $\mathbb{R}_+^2 = \mathbb{R}_+\times\mathbb{R}_+$, where $\mathbb{R}_{+}$ is the collection of all positive real numbers.

Lemma 3.1 (It$\hat{o}$ 's formula ^[23]) Let

be a $d$-dimensional Itô process, where $f\in\mathcal{L}^1(\mathbb{R}_+; \mathbb{R}^d)$, $g(t)\in\mathcal{L}^2(\mathbb{R}_+; \mathbb{R}^{d\times m})$. Let $V(X, t)\in C^{2, 1}(\mathbb{R}^d\times \mathbb{R}_+; \mathbb{R})$. Then the process $V(X, t)$ is again an Itô process, whose differential equation is given by

Definition 3.1 (Lv and Wang ^[24]) A stochastic system is said to be almost surely stochastically permanent if for any initial value $x_ 0\in \mathbb{R}^n_+$, the solution $x(t) = (x_1(t), x_2 (t), ..., x_n (t))$ satisfies

Lemma 3.2 (Lipster and Shiryayev ^[25]) Let $A(t)$ and $U(t)$ be two continuous adapted increasing process on $t\geq 0$ with $A(0) = U(0) = 0$ a.s. Let $M(t)$ be a real-valued continuous local martingale with $M(0) = 0$ a.s. Let $X_0$ be a nonnegative $\mathcal{F}_0$ -measurable random variable such that $E X_0 < \infty$. Define $X(t) = X_0+A(t)-U(t)+M(t)$ for all $t\geq0$. If $X(t)$ is nonnegative, then

Lemma 3.3 For the linear stochastic differential equation

Let $x(0) = x_0 > 0$, Eq (3.1) is almost surely stochastically permanent, i.e, $x(t)$ to Eq (3.1) satisfies

Proof. Applying Itô's formula to Eq (3.1) leads to

Integrating both sides from $0$ to $t$ gives

If we choose sufficiently small positive constant $c$ and find a suitable positive constant $C_1$, we can get that

Hence

Denote $Z_1(t) = x_0+C_1+\int_0^t\rho e^{[c(s-t)]}x(s)\mathrm{d}B(s)$. It follows from Lemma 3.2 that

On the other hand, by Itô's formula, we obtain

Setting $y(t) = x^{-1}(t)$, we can rewrite the above equations as

We deduce

Then,

We can choose suitable constant $C_2$ such that

Denote $Z_2(t) = y_0+C_2-\int_0^te^{(s-t)}\rho y(s)\mathrm{d}B(s).$ By Lemma 3.2, we have

Consquently

Namely

From (3.3) and (3.4) we know that Eq (3.1) is almost surely stochastically permanent. The lemma is completed.

3.1.   The existence and uniqueness of the global positive solution

We first demonstrate that under some mild conditions, there exists a unique global positive solution ${\bf Y}(t)$ of system (2.2) on time $t\geq0$ almost surely.

Theorem 3.1. For any given initial value (2.3), there exists a unique global positive solution ${\bf Y}(t) \in\mathbb{R}_+^2$ of system (2.2) on time $t\geq 0$ almost surely, that is, the solution will remain in $\mathbb{R}_+^2$ with probability 1.

Proof. As the system (2.2) has locally Lipschitz continuous coefficients, for any given initial value (2.3), the system (2.2) admits a unique maximal local solution ${\bf Y}(t)$ on $t\in[-\tau, \tau_e)$, where $\tau_e$ is the explosion time ^[26]. Since

and

it is easy to know that $G(t) > 0$, $I(t) > 0$ for the given initial value (2.3). i.e. the system (2.2) has a unique positive local solution ${\bf Y}(t)$ on $t\in[-\tau, \tau_e)$.

In order to verify the solution is global, we only need to prove that $\tau_e = \infty$ a.s.

From the first equation of system (2.2), we have

Let

Thus, $\Phi (t)$ is the solution of the following system

By comparison principle for stochastic differential equations ^[27], we obtain that

Lemma 3.2 indicates the system (3.5) is is almost surely stochastically permanent. Hence there exists a positive constant $\Phi_M$ satisfied $\Phi (t)\leq\Phi_M$ for all $t\geq 0$ a.s.

Similarly, from the second equation of (2.2), we can get

Then

where $\Psi_M$ is a positive constant and

On the other hand,

then

where $\varphi_m$ is a positive constant and

Then,

where $\psi_m$ is a positive constant and

From above, we have

Notice that $\varphi (t), \Phi (t), \psi(t)$ and $\Psi(t)$ are the solution of linear stochastic differential equation, and they all exist for $t\geq 0$. Hence, for the given initial value (2.3), the system (2.2) has a unique positive global solution ${\bf Y}(t)$ on $t\geq 0$. This completes the proof of Theorem 3.1.

From the above argument, we know that

and Lemma 3.2 admits that the linear stochastic differential system is almost surely stochastically permanent, we can also get the following theorem:

Theorem 3.2. System (2.2) is almost surely stochastically permanent, i.e, ${\bf Y}(t)$ to Eq (2.2) satisfies

Let $\Gamma = \{{\bf Y}(t)\in \mathbb{R}_+^2:0 < G(t), I(t)\leqslant N, t\geq 0\}$ denote the the positive invariant set of the system (2.1), where $N = \max\{\Phi_M, \Psi_M\}$, where $\Phi_M, \Psi_M$ is defined as Theorem 3.1. From Theorem 3.2, it's easy to know that $\Gamma$ is also a positive invariant set of the stochastic delayed model (2.2) almost surely, i.e., if ${\bf Y}(0) \in \Gamma$, then $P\left({\bf Y}(t) \in\Gamma\right) = 1, t > 0$.

3.2.   Analysis of asymptotic behavior

According to ^[10], there is a positive equilibrium state $E_+(G^*, I^*)$ in the deterministic system (2.1). However, it is no longer the equilibrium point of the corresponding stochastic system (2.2). In this section, we will study the asymptotic behavior of the solution of system (2.2) around the point $E_+(G^*, I^*)$.

Noting that $f'(x) = \frac{ra^rx^{r-1}}{(a^r+x^r)^2}$, it is easy to get $M'$, the maximum value of $f'(x)$, as follows

by the extremum theorem.

Theorem 3.3. If there exist $u > 0$ and $\varepsilon > 0$ such that

and

then for any ${\bf Y}(0)\in\Gamma$, the following results holds:

where $M_u = \max \{ \sigma M', uS_i \Phi_M \}$, and

and

Proof. Define

where $C\triangleq RL/2$, $D\triangleq KL/2$, $L\triangleq\sigma M'$, $R\triangleq S_g+S_iI^*$, $K\triangleq S_i\Phi_M$.

Let $V({\bf Y}(t)) = W({\bf Y}(t))+U({\bf Y}(t))$. It is straightforward to check that $V({\bf Y}(t))$ is a Liapunov function. By the mean value theorem, we get

where the value $\xi$ is between $G(t-\tau)$ and $G^*$. By the inequality $2ab\leqslant (a^2+b^2)$, we deduce

Applying the Itô's formula and the inequality $ab\leq\frac{\varepsilon}{2}a^2+\frac{1}{2\varepsilon}b^2$ for $\varepsilon > 0$, we obtain

Besides,

Hence,

Taking integral from 0 to $t$ and expectations, we get

where

Note that

where $M = \max\left\{\big\lvert M_I-I^*\big\rvert\alpha_1 \Phi_M, \big|I(0)-I^*\big|\alpha_1 G(0)\right\}$.

Therefore,

Denote

From the conditions of the theorem, it is easy to know $H_1 > 0$ and $H_2 > 0$. Taking the limit superior of both sides of (3.10) leads to

Substituting $R$, $K$ into (3.11), we can obtain the conclusion (3.7). This completes the proof of Theorem 3.3.

Remark 3.1. The stochastic system (2.2) does not have a positive equilibrium point, but the solutions of the system perturbs around a point near the positive equilibrium point of its corresponding deterministic system. From the Theorem 3.3, we deduce that the smaller the stochastic disturbance is, the closer the point is to the positive equilibrium point, and the vibration amplitude decreases as the disturbance decreases.

3.3.   The parameters estimation

Although the stochastic system (2.2) is constructed based on the deterministic system (2.1), it can be identified directly from the data without the necessity of finding (2.1) in advance. We propose to estimate the parameters in (2.2) by a maximum likelihood estimation (MLE) method. Let $\theta = (b, S_g, S_i, d_i, a, r, \sigma, \tau, \alpha_1, \alpha_2)^{\top}$ be the collection of parameters in the system (2.2). Recall that ${\bf Y}(t)$ = $(G(t), I(t))^{\top}$. Let $h({\bf Y}(t), \theta) = (b-S_gG(t)-S_iG(t)I(t), \sigma f(G(t-\tau))-d_iI(t))^{\top}$, $B(t) = (B_1(t), B_2(t))$, and $\phi(Y(t), \theta) = (\alpha_1 G(t), \alpha_2 I(t))^{\top}$. Then the system (2.2) can be represented as

As $\Delta t\rightarrow 0$, the Euler scheme produces the following discretization:

Since $B_{1}(t)$ and $B_{2}(t)$ are two independent Brownian motions, ${\bf Y}(t+\Delta t) - {\bf Y}(t)$ are a two-dimensional Gaussian independent variables with mean $h({\bf Y}(t), \theta) \Delta t$ and covariance matrix

Let $\delta_i(\theta) = {\bf Y}(t_i)- {\bf Y}(t_{i-1})- h({\bf Y}(t_{i-1}), \theta)(t_i-t_{i-1})$, where ${\bf Y}(t_0) = y_0$. We obtain the following log-likelihood function:

By maximizing the above log-likelihood function with respect to $\theta$, we can obtain a consistent estimator, denoted by $\widehat{ \theta}$, under some mild conditions ^[28].

In fact, Li et al. ^[10] showed that the delay of insulin secretion can be estimated by observing the second peak of insulin secretion. Physiological facts ensure that both glucose and insulin return to their basal level after about three hours, which allows two parameters can be expressed by other parameters. These will reduce the number of parameters in the parameter space $\theta$ and hence improves the estimation accuracy of $\widehat{\theta}$.

4.   Simulations

In this section, in order to illustrate our theoretical findings, we provide some numerical simulations of system (2.2) and the solution of the corresponding deterministic system (2.1). Simulations are performed by using Matlab Euler Maruyama method for SDE. We used experiment data listed in Tables 1 and 2 in ^[10], which originated from ^[4] and ^[8]. We consider the typical data set from subjects 6, 7 and 27.

In section 3.3, we proposed a MLE method to estimate the parameters in the system (2.2). However, since the number of data points for each subject is relatively small (around 20) and the log-likelihood function is considerably complicated, it is challenging to get a reasonable estimate of all parameters in $\theta$. Instead, we fix the values of $(b, S_g, S_i, d_i, a, r, \sigma, \tau)$ as suggested in ^[10] and maximize the log-likelihood function with respect to $(\alpha_1, \alpha_2)$ only. And we also consider the delay as approximating the time between the primary insulin release and the trough in insulin concentration determined by its secondary release. The values of $(b, S_g, S_i, d_i, a, r, \sigma, \tau)$ and the maximum likelihood estimates of $(\alpha_1, \alpha_2)$ for each subject are listed in Table 1.

To fit the data from subject 6 with the system (2.2), we choose $u = 9.2$, $\varepsilon = 1.3$ in Theorem 3.3 and consider three different sets of $(\alpha_1, \alpha_2)$: $(0.005, 0.005)$, $(0.012, 0.087)$ and $(0.07, 0.175)$, stimulating very small, adequate, and too large stochastic disturbances, respectively. By calculation, it can be seen that the conditions of Theorem 3.3 are satisfied by using these first two sets of $(\alpha_1, \alpha_2)$ but not satisfied by using the third set.

For each set of $(\alpha_1, \alpha_2)$, we conducted 1000 Monte-Carlo simulations. In each simulation, we compute the sum of the root mean square errors (RMSE) of glucose and insulin. We calculate the percentage of simulations that reduce the sum of RMSEs of glucose and insulin as compared to the system (2.1), and then compute the average reduction and the maximum reduction of those simulations. Moreover, we calculate the ratios of the average reduction and the maximum reduction to the sum of RMSEs of glucose and insulin in the system (2.1). The simulation results are reported in the following table.

The simulation results suggest that the intensity of stochastic components in the system (2.2) needs to be carefully chosen. On one hand, if the intensity is too small, the improvement of model fitness is not sufficient. On the other hand, the over-large intensity may fail to improve the model fitness. In addition, the simulation results also indicate that given a reasonable set of values of the parameters $(b, S_g, S_i, d_i, a, r, \sigma, \tau)^{\top}$, the proposed MLE method can provide satisfactory estimates of $(\alpha_1, \alpha_2)$ and the solution of the proposed system (2.2) can significantly improve the model fitness with certain probabilities.

Figures 1–3 depict the solution curves of the system (2.2) from three simulation with $(\alpha_1, \alpha_2) = (0.005, 0.005), \; (0.012, 0.087)$ and $(0.07, 0.175)$, respectively. By comparing the three figures, we note that as the values of $\alpha_1$ and $\alpha_2$ get larger, the disturbance amplitudes of the solution curves of the stochastic system (2.2) are enlarged and consequently depart further away from the positive equilibrium point of system (2.1). This phenomena is expected as $\alpha^2$'s control the intensity of disturbance, and confirms our theoretical findings in Theorem 3.3. In the process of IVGTT, there may be many disturbance of human activity that are hard to control such as psychological pressure, emotions, excessive morning exercises and so on. Our simulations also suggest that if the disturbances of human activity are enough large, there may be larger error in diagnostic results provided by a deterministic model.

To fit the data from subjects 7 and 27 with the system (2.2), we conducted 1000 Monte-Carlo simulations of the proposed system with the corresponding MLE estimates of $(\alpha_1, \alpha_2)$. Figures 4 and 5 display the solution curves from one simulation for subject 7 and one simulation for subject 27.

In each simulation, we also compare the resulting model fitness to that of the deterministic system (2.1). The summary statistics for subjects 7 and 27 are presented in Table 3. The results again show that with the $(\alpha_1, \alpha_2)$ estimated by the MLE method, the solution of the proposed (2.2) can significantly provide better data fitting than the deterministic model with certain probabilities. In fact, we have used our model to the data from many other subjects who originally appeared in ^[4,8] and obtained similar observations. We speculate that if the number of data points is sufficiently large so that all parameters in the system (2.2) can be estimated by the MLE method, the proposed stochastic model and the associated deterministic model with the MLE estimates can provide better data fitting than many existing models ^[4,8,10,12]. We will explore this estimation of all parameters in $\theta$ in future work.

5.   Discussion

Our work in this paper is an initial investigation of the stochastic dynamic modeling for the data from the protocol of intravenous glucose tolerance test. Our approach includes formulating a well-post stochastic dynamic model (2.2) according to a deterministic model (2.2) and then developing a maximum log-likelihood estimation method to improve the data fitting of the deterministic model (2.1). In general, deterministic model and stochastic model are alternate viewpoints on the same physiological metabolic phenomenon and offer complementary insights ^[29]. (More detail discussions regarding to the relations of deterministic models and stochastic models with various examples can be found in ^[29].) The model (2.2) can be viewed as a more flexible platform built upon the corresponding model (2.1) by allowing the stochastic rate of $S_g$ and $d_i$. The proposed model can possibly improve the estimations of the parameters by the MLE method studied in section 3.3. Our simulation studies numerically confirm our theoretical findings and demonstrate that the proposed model with estimated parameters can improve the fitness of clinical data.

Parameter estimation for dynamical system models has been a challenging problem in fitting data for real life problems. Being able to estimate insulin sensitivity and glucose effectiveness for an individual with great accuracy is extremely important for research in finding the pathways to T2DM and drug development. We seemly improved the fitness of IVGTT data for this protocol. Nevertheless, due to the number of each set of data points is relatively small for the number of parameters, and the log-likelihood function is noticeably complicated, it is challenging to get a reasonable estimate of all parameters in $\theta$, even though it is possible to reduce the number of parameters to be estimated as shown in ^[10]. We will continue to explore this approach for large set of parameters and/or reduced number of parameters in future work to fulfill the demand of accuracy of parameter estimation.

Acknowledgments

We would like to thank the anonymous referees for their careful reading of the original manuscript and their many valuable comments and suggestions that greatly improve the presentation of this work. This work is supported by the National Natural Science Foundation of China (No. 11701495), Scientific and Technological Key Projects of Henan Province (No. 192102310193) and Nanhu Scholars Program for Young Scholars of XYNU.

Conflict of interest

The authors declare that they have no competing interests.