糖尿病(DM) 是一种代谢紊乱引起的疾病, 其产生的主要原因为胰岛素绝对或相对分泌不足和(或) 胰岛素利用障碍. 糖尿病主要分为三类, 包括1型糖尿病(T1DM), 2型糖尿病(T2DM) 和妊娠糖尿病. 目前糖尿病的诊断方法主要为空腹血糖测定和糖耐量试验^[1-2].

糖尿病治疗的目标是控制患者血糖浓度维持在正常范围内, 常用的治疗方式是服用降糖药和注射胰岛素, 其中连续地皮下胰岛素注射治疗可以通过安置胰岛素泵(开环方法) 来实现, 但是传统的治疗方式容易使患者发生低血糖或者因血糖控制不佳而导致各种并发症, 并危及患者的生命. 为了解决这个问题, 学者们开始考虑使用现代控制技术制造出能够代替胰腺的"人工胰腺", 通过预测血糖的变化来给患者供给适量的胰岛素, 最终达到控制血糖的目的. 为此, 学者们提出了不同的数学模型来模拟血糖-胰岛素相互间的作用, 希望为糖尿病提供更有效、更经济的疗法. 数学模型已成为了解葡萄糖-胰岛素调节系统的重要工具, 并给糖尿病患者提供更合理的胰岛素使用方法. 基于Bolie^[3]的开创性工作, 许多著名的胰岛素-葡萄糖模型相继被研究^[4-9]. 例如, Ackerman等^[4-5]建立了一种口服葡萄糖耐量测试模型, 该模型可测量患者利用特定量葡萄糖的能力. Li等模拟了葡萄糖-胰岛素调节系统^[8]并研究了具有两个时滞的胰岛素调节系统^[9]. Topp等^[10]通过考虑包含$ \beta $细胞质量, 胰岛素和血糖相互作用的血糖- 胰岛素调节模型, 给出了预测长时间高血糖症存在的三种途径. Hernandez等^[11]考虑了胰岛素利用障碍, 进一步扩展了Topp模型. Gallenberger等^[12]提出了一个与$ \beta $细胞相关的周期葡萄糖-胰岛素浓度动力学模型. 其他有关糖尿病的预防和血糖监测的研究可参考文献[13-19]. 在本文中, 基于Li, Wang和Huang等^{[8-9, 20-21]}提出的模型, 考虑以下模型

(1.1) $ \begin{equation} \left\{ \begin{array}{llc} { } {\rm d}G(t) = \Big(G_{in}-S_gG(t)-S_iG(t)\frac{I(t)}{n+I(t)}+\frac{S_k}{r^m+I^m(t)}\Big){\rm d}t, \\ { } {\rm d}I(t) = \Big(\frac{\sigma G^2(t)}{\alpha^2+G^2(t)}-d_iI(t)\Big){\rm d}t, \end{array} \right. \end{equation} $

其中$ G_{in} $是葡萄糖平均输入速率, $ \frac{\sigma G^2(t)}{\alpha^2+G^2(t)} $是受葡萄糖浓度升高刺激而分泌的胰岛素, $ S_gG(t) $是不依赖胰岛素的葡萄糖摄取, $ S_iG(t)\frac{I(t)}{n+I(t)} $代表依赖胰岛素的葡萄糖摄取, $ \frac{S_k}{r^m+I^m(t)} $为肝糖原水解产生的葡萄糖, $ d_iI(t) $表明胰岛素以恒定速率$ d_i $降解.

众所周知, 饮食、运动、年龄和体重等会影响胰岛素-葡萄糖动力学. 事实上, 血糖的摄取不仅取决于某些激素(例如皮质醇或生长激素) 的浓度变化, 还取决于食物摄取的频率、身体活动的突然变化或情绪压力. 代谢器官和组织还受到各种内部和外部的影响, 这些影响会随时间而变化(例如: 血流量、能量需求、激素水平、组织本身的细胞代谢). Liu等^[22]指出葡萄糖耐受性, 胰岛素对葡萄糖激发的反应, 胰岛素敏感性和$ \beta $细胞的形态均会受到随机波动的影响. 同时, 与确定性系统相比, 随机模型能更好地描述实际情况. 因此, 我们利用Yuan等^[23]使用的技术和方法, 考虑如下随机模型

(1.2) $ \begin{equation} \left\{ \begin{array}{llc} { } {\rm d}G(t) = \Big(G_{in}-S_gG(t)-S_iG(t)\frac{I(t)}{n+I(t)}+\frac{S_k}{r^m+I^m(t)}\Big){\rm d}t+\alpha_1G(t){\rm d}B_1(t), \\ { } {\rm d}I(t) = \Big(\frac{\sigma G^2(t)}{\alpha^2+G^2(t)}-d_iI(t) \Big){\rm d}t+\alpha_2I(t){\rm d}B_2(t), \end{array}\right. \end{equation} $

其中$ B_i(t)(i = 1, 2) $为定义在完备概率空间$ (\Omega, {\cal F}, \{{\cal F}_t\}_{t\geq0}, P) $中的相互独立布朗运动, 非负常数$ \alpha_i(i = 1, 2) $为噪声强度.

在本文中, 除非另有说明, $ (\Omega, {\cal F}, \{{\cal F}_t\}_{t\geq0}, P) $为完备的概率空间, 相应的滤子$ \{{\cal F}_t\}_{t\geq0} $满足通常条件(即$ {\cal F}_t $是右连续的且$ {\cal F}_0 $包含所有的$ {\Bbb P} $零集).

本文安排如下: 在第2节, 得到系统正平衡点的存在性并证明其是全局渐进稳定的. 第3节研究随机系统的全局正解的存在唯一性、随机持久性和平稳分布的存在性. 第4节通过数值模拟来验证理论分析的结果. 在第5节, 我们进行简要的讨论并总结主要的结果.

确定性模型(1.1) 的正平衡点$ E(G^*, I^*) $满足

(2.1) $ \begin{equation} G_{in}-S_gG^*-S_iG^*\frac{I^*}{n+I^*} +\frac{S_k}{r^m+(I^*)^m} = 0, \frac{\sigma (G^*)^2}{\alpha^2+(G^*)^2}-d_iI^* = 0. \end{equation} $

由方程(2.1) 可得

(2.2) $ \begin{equation} G_{in}-S_gG^*-S_iG^*\frac{\frac{\sigma (G^*)^2}{d_i(\alpha^2+(G^*)^2)}}{n+ \frac{\sigma (G^*)^2}{d_i(\alpha^2+(G^*)^2)}} +\frac{S_k}{r^m+(\frac{\sigma (G^*)^2}{d_i(\alpha^2+(G^*)^2)})^m} = 0. \end{equation} $

令

$ f(x) = G_{in}-S_gx-S_ix\frac{\frac{\sigma x^2}{d_i(\alpha^2+x^2)}}{n+ \frac{\sigma x^2}{d_i(\alpha^2+x^2)}} +\frac{S_k}{r^m+(\frac{\sigma x^2}{d_i(\alpha^2+x^2)})^m}, $

易得, $ f(x) $为单调减函数. 又$ f(0) = G_{in}+\frac{S_k}{r^m}>0, \lim\limits_{x\rightarrow +\infty}f(x) = -\infty, $因此, $ f(x) = 0 $有唯一的正根, 即系统(1.1) 存在唯一的正平衡点$ E(G^*, I^*) $. 注意到系统(1.1) 在$ E(G^*, I^*) $处的Jacobian矩阵为

$ J(E) = \left( \begin{array}{cc} { } -S_g-\frac{S_iI^*}{n+I^*} &{ } -\frac{S_iG^*}{n+I^*}(1-\frac{I^*}{n+I^*}) { } -\frac{mS_k(I^*)^m}{(r^m+(I^*)^m)^2I^*} \\ { } \frac{2\sigma G^*}{(G^*)^2+\alpha^2}(1-\frac{(G^*)^2}{(G^*)^2+\alpha^2}) & -d_i \end{array} \right), $

其特征方程为$ \lambda^2+p_1\lambda+p_2 = 0, $其中

$ p_1 = d_i+S_g+\frac{S_iI^*}{n+I^*}>0, $

$ \begin{eqnarray*} p_2& = &d_i(S_g+\frac{S_iI^*}{n+I^*})+\frac{2\sigma G^*}{(G^*)^2+\alpha^2}\frac{S_iG^*}{n+I^*}(1-\frac{(G^*)^2}{(G^*)^2+\alpha^2})(1-\frac{I^*}{n+I^*})\\ &&+\frac{2\sigma G^*}{(G^*)^2+\alpha^2}(1-\frac{(G^*)^2}{(G^*)^2+\alpha^2})\frac{mS_k(I^*)^m}{(r^m+(I^*)^m)^2I^*}>0. \end{eqnarray*} $

由Routh-Hurwitz判据可得, 系统(1.1) 的正平衡点$ E(G^*, I^*) $为局部渐进稳定的.

定理2.1  系统(1.1) 存在唯一的正平衡点, 且该平衡点是局部渐进稳定的.

下面我们将证明其全局稳定性.

定义

$ B(G, I) = \frac{1}{I}, {\quad} f_1 = G_{in}-S_gG-S_iG\frac{I}{n+I}+\frac{S_k}{r^m+I^m}, $

$ f_2 = \frac{\sigma G^2}{\alpha^2+G^2}-d_iI. $

计算可得

$ \frac{\partial (Bf_1)}{\partial G} +\frac{\partial (Bf_2)}{\partial I} = \frac{-S_g-\frac{S_iI}{n+I}}{I}-\frac{\sigma G^2}{I^2(G^2+\alpha^2)}<0. $

由Dulac判别法可知, 系统(1.1) 无闭轨线, 因此正平衡点$ E(G^*, I^*) $全局渐进稳定.

定理2.2  系统(1.1) 存在唯一的正平衡点, 且该平衡点是全局渐进稳定的.

本节我们将研究随机模型(1.2) 的动力学. 首先我们证明对于任意的正初始值, 系统(1.2) 存在唯一的全局正解.

定理3.1  对于任意给定的初始条件$ \big(G(0), I(0)\big)\in {\Bbb R}_+^2 $, 系统(1.2) 存在唯一的正解$ \big(G(t), I(t)\big) $, 满足对任意的$ t\geq0 $, 有$ G(t)>0 $且$ I(t)>0 $是几乎确定的(a.s.).

证  显然, 系统(1.2) 系数满足局部Lipschitz条件, 因此对任意的$ \big(G(0), I(0)\big)\in {\Bbb R}^2_+ $, 系统在$ t\in [0, \tau_e) $时存在唯一正解$ \big(G(t), I(t)\big) $, 其中$ \tau_e $为爆破时间. 为了证明这个解是全局的, 只需证明$ \tau_{e} = +\infty $, a.s..

令$ k_{0}>0 $足够大, 使得$ (G(t), I(t))\in[\frac{1}{k_0}, k_{0}]\times[\frac{1}{k_0}, k_{0}] $, 对于任意的正数$ k\geq k_{0} $, 定义一个停时序列$ \tau_k = \inf\left\{t\in [0, \tau_e):G(t)\notin (\frac{1}{k}, k);I(t)\notin (\frac{1}{k}, k)\right\}. $

定义$ \inf\Phi = +\infty $ ($ \Phi $代表一个空集). 显然当$ k\rightarrow \infty $时, $ \tau_{k} $是单调递增的, 且$ \tau_{k}<\tau_{e} $, 于是有$ \tau_{\infty} = \lim\limits_{k\rightarrow +\infty}\tau_{k} $, 其中$ \tau_{\infty}\leq\tau_{e} $ a.s.. 因此, 只需证明$ \tau_{\infty}\rightarrow \infty $ a.s..

假设$ \tau_{\infty}\nrightarrow\infty $, 则存在常数$ T_1\geq0 $, $ \epsilon\in(0, 1) $和一个正整数$ k_{1}\geq k_{0} $, 使得

(3.1) $ \begin{equation} P\{\tau_{k} \leq T_1\}\geq\varepsilon , \forall k\geq k_{1}. \end{equation} $

定义如下$ V(G, I) \in C^2({\Bbb R}_+^2, {\Bbb R}_+) $:

$ V(G, I) = (G-\ln G-1)+(I-\ln I-1). $

由It$ \hat{\rm o} $公式可得

$ {\rm d}V = LV{\rm d}t+\alpha_1(G-1)dB_1+\alpha_2(I-1){\rm d}B_2, $

其中

$ \begin{eqnarray*} LV& = &(1-\frac{1}{G})(G_{in}-S_gG-\frac{S_iGI}{n+I}+\frac{S_k}{r^m+I^m}) +(1-\frac{1}{I})(\frac{\sigma G^2}{\alpha^2+G^2}-d_iI)+\frac{1}{2}(\alpha_1^2+\alpha_2^2)\\ &\leq&(G_{in}-S_gG+\frac{S_k}{r^m+I^m}+S_g+\frac{S_i I}{G+I}) +(\frac{\sigma G^2}{\alpha^2+G^2}-d_iI)+d_i +\frac{1}{2}(\alpha_1^2+\alpha_2^2)\\ &\leq&G_{in}+\frac{S_k}{r^m}+S_g+S_i+\sigma+d_i+\frac{1}{2}(\alpha_1^2+\alpha_2^2)\\ :& = &M, \end{eqnarray*} $

其中$ M $为正常数. 因此

(3.2) $ \begin{equation} {\rm d}V = LV{\rm d}t+\alpha_1(G-1)dB_1+\alpha_2(I-1){\rm d}B_2\leq M{\rm d}t+\alpha_1(G-1){\rm d}B_1+\alpha_2(I-1){\rm d}B_2. \end{equation} $

对式(3.2) 两边同时从$ 0 $到$ \tau_k\wedge T $积分可得

(3.3) $ \begin{equation} V[G(\tau_k\wedge T), I(\tau_k\wedge T)]\leq V[G(0), I(0)]+\int_0^{\tau_k\wedge T} M {\rm d}s+M_1+M_2, \end{equation} $

其中$ M_1 = \int_0^{\tau_k\wedge T}{\alpha_1(G-1)}{\rm d}B_1, M_2 = \int_0^{\tau_k\wedge T}{\alpha_2(I-1)}{\rm d}B_2. $

对式(3.3) 两边同时取期望可得

(3.4) $ \begin{equation} E[V(G(\tau_k\wedge T), I(\tau_k\wedge T))]\leq V[G(0), I(0)]+M(\tau_k\wedge T). \end{equation} $

令$ \Omega_k = \{\tau_k\leq T\} $($ k\geq k_1 $), 则由式(3.1), $ P(\Omega_k)\geq\varepsilon $. 注意到, 对任意的$ \omega\in \Omega_k $, 至少有一个$ G(\tau_k, \omega) $或$ I(\tau_k, \omega) $等于$ k $或$ \frac{1}{k} $, 因此$ V(G, I) $不小于$ k-\ln k-1 $或$ \frac{1}{k}+\ln k-1 $, 即

$ V\big(G(\tau_k, \omega), I(\tau_k, \omega)\big)\geq (k-\ln k-1) \wedge (\frac{1}{k}+\ln k-1). $

由式(3.4) 可知

$ \begin{eqnarray*} V\big(G(0), I(0)\big)+MT\geq &E\Big[1_{\Omega_k}(\omega)V\Big(G(\tau_k, \omega), I(\tau_k, \omega)\Big)\Big]\\ \geq &\varepsilon [(k-\ln k-1) \wedge (\frac{1}{k}+\ln k-1)], \end{eqnarray*} $

其中$ 1_{\Omega_k} $是$ \Omega_k $的示性函数. 令$ k\rightarrow \infty $可得$ \infty>V\big(G(0), I(0)\big)+MT = \infty. $因此有$ \tau_\infty = \infty $, 系统(1.2) 存在唯一的全局正解.

注3.1  定理3.1表明, 对于任意的正初始值, 系统(1.2) 依概率1存在唯一正解. 由系统(1.2) 可知

$ \begin{eqnarray*} {\rm d}G(t)& = &\Big(G_{in}-S_gG(t)-S_iG(t)\frac{I(t)}{n+I(t)}+\frac{S_k}{r^m+I^m(t)}\Big){\rm d}t+\alpha_1G(t){\rm d}B_1(t)\\ &\geq&\Big(G_{in}-(S_g+S_i)G(t)\Big){\rm d}t+\alpha_1G(t){\rm d}B_1(t). \end{eqnarray*} $

我们考虑如下方程

$ \left\{ \begin{array}{lll} {\rm d}X(t) = \Big(G_{in}-(S_g+S_i)X(t)\Big){\rm d}t+\alpha_1X(t){\rm d}B_1(t), \\ X(0) = G(0). \end{array} \right. $

由文献[24]可知$ \lim\limits_{t\rightarrow \infty}\frac{1}{t}\int_{0}^t X(s){\rm d}s = \frac{G_{in}}{S_g+S_i}, $即$ { } \liminf\limits_{t\rightarrow \infty} X(t)>0. $

通过随机微分方程比较定理可得$ { } \liminf\limits_{t\rightarrow \infty} G(t)>0. $此外, 由定理3.1可得$ t\in(0, \infty) $时, $ G(t)>0 $. 因此, 对任意的$ t\in(0, \infty) $, 存在常数$ \ell>0 $使得

(3.5) $ \begin{equation} G(t)>\ell. \end{equation} $

注3.2  通常而言, 相较于唯一的全局正解, 系统的持久性更有实际意义, 因为它预示着糖尿病病症动力学的长期水平. 因此, 我们将在下一部分探讨随机模型(1.2) 的持久性.

本节将讨论随机系统(1.2) 的持久性. 首先, 我们给出随机持久的定义.

定义3.1^[25]  对于任意给定的初始条件$ \big(G(0), I(0)\big)\in {\Bbb R}_+^2 $, 系统(1.2) 存在唯一的正解$ \big(G(t), I(t)\big) $, 如果对任意$ 0<\epsilon <1 $, 存在正常数$ \delta = \delta(\epsilon) $和$ \bar{\delta} = \bar{\delta}(\epsilon) $使得下式成立

(3.6) $ \begin{equation} \begin{array} {ll} { } \liminf\limits_{t\rightarrow \infty} P\{G(t)\geq\delta\}\geq1-\epsilon, \quad \quad \liminf\limits_{t\rightarrow \infty} P\{I(t)\geq\delta\}\geq1-\epsilon, \\ { } \liminf\limits_{t\rightarrow \infty} P\{G(t)\leq\bar{\delta}\}\geq1-\epsilon, \quad \quad \liminf\limits_{t\rightarrow \infty} P\{I(t)\leq\bar{\delta}\}\geq1-\epsilon, \end{array} \end{equation} $

则称随机系统(1.2) 为随机持久的.

定理3.2  随机系统(1.2) 是随机持久的.

证  设$ (G(t), I(t)) $是系统(1.2) 具有任意正初始值$ (G(0), I(0)) \in {\Bbb R}^2_+ $的解, 我们将定理的证明分为两步.

首先我们证明$ G(t) $为随机持久的.由定理3.1可知, 对任意$ t\geq 0 $, 系统(1.2) 的解$ (G(t), I(t)) \in {\Bbb R}^2_+ $ a.s.. 定义如下$ V(G(t), t) \in C^2({\Bbb R}_+^2, {\Bbb R}_+) $:

$ V = e^{\lambda t} G(t), $

其中$ 0<\lambda<S_g $. 由It$ \hat{\rm o} $公式可得

$ \begin{eqnarray*} {\rm d}V& = &e^{\lambda t}{\rm d}G+\lambda e^{\lambda t}G{\rm d}t\\ & = &e^{\lambda t}\Big(G_{in}-S_gG-S_iG\frac{I}{n+I}+\frac{S_k}{r^m+I^m}+\lambda G\Big){\rm d}t +\alpha_1e^{\lambda t}G {\rm d}B_1\\ & <&e^{\lambda t}\Big(G_{in}+\frac{S_k}{r^m}\Big){\rm d}t+\alpha_1e^{\lambda t}G {\rm d}B_1. \end{eqnarray*} $

对两边从$ 0 $到$ t $积分得

$ \begin{eqnarray*} e^{\lambda t}G(t)-G(0)<\frac{1}{\lambda}\Big(G_{in}+\frac{S_k}{r^m}\Big)(e^{\lambda t}-1)+\int_0^t \alpha_1e^{\lambda s}G {\rm d}B_1, \end{eqnarray*} $

两边同时取期望可得

$ e^{\lambda t}E\big(G(t)\big)<G(0)+\frac{1}{\lambda}\Big(G_{in}+\frac{S_k}{r^m}\Big)(e^{\lambda t}-1), $

即$ E\big(G(t)\big)<G(0)e^{-\lambda t}+\frac{1}{\lambda}\big(G_{in}+\frac{S_k}{r^m}\big)(1-e^{-\lambda t}), $因此

$ \limsup\limits_{t\rightarrow \infty}E\big(G(t)\big) <\frac{1}{\lambda}\Big(G_{in}+\frac{S_k}{r^m}\Big). $

对任意$ \epsilon_1>0 $, 取$ \bar{\delta}_1 = \frac{\frac{1}{\lambda}\big(G_{in}+\frac{S_k}{r^m}\big)}{\epsilon_1} $, 由Chebyshev不等式可得

$ \limsup\limits_{t\rightarrow \infty} P\{G(t)>\bar{\delta}_1\}\leq \limsup\limits_{t\rightarrow \infty} \frac{E\big(G(t)\big)}{\bar{\delta}_1}< \epsilon_1. $

于是有

(3.7) $ \begin{equation} \liminf\limits_{t\rightarrow \infty} P\{G(t)\leq\bar{\delta}_1\}\geq1-\epsilon_1. \end{equation} $

即$ G(t) $是随机有界的. 此外，由It$ \hat{\rm o} $公式易得

$ \begin{eqnarray*} {\rm d}G^{-1}& = &\Big\{-G^{-2} \Big(G_{in}-S_gG-S_iG\frac{I}{n+I}+\frac{S_k}{r^m+I^m}\Big) +\frac{\alpha_1^2}{G}\Big\}{\rm d}t-\alpha_1G^{-1}{\rm d}B_1\\ &\leq&\Big(-\frac{G_{in}}{G^2}+\frac{S_g+S_i+\alpha_1^2}{G}\Big){\rm d}t -\alpha_1G^{-1}{\rm d}B_1. \end{eqnarray*} $

令$ W = \frac{e^t}{G(t)} $, 则对$ W $应用It$ \hat{\rm o} $公式可得

$ {\rm d}W = \frac{e^t}{G(t)}{\rm d}t+e^t{\rm d}G^{-1} \leq e^t\Big(-\frac{G_{in}}{G^2}+\frac{1+S_g+S_i+\alpha_1^2}{G}\Big){\rm d}t -\alpha_1e^tG^{-1}{\rm d}B_1. $

显然, 存在常数$ M' = \frac{(1+S_g+S_i+\alpha_1^2)^2}{4G_{in}}>0 $使得对任意$ G>0 $, $ -\frac{G_{in}}{G^2}+\frac{1+S_g+S_i+\alpha_1^2}{G}< M' $. 因此

$ {\rm d}W<M'e^t{\rm d}t-\alpha_1e^tG^{-1}{\rm d}B_1. $

对两边从$ 0 $到$ t $积分并取期望可得

$ \limsup\limits_{t\rightarrow \infty} E\big(W(t)\big)<W(0)+M'(e^t-1), $

即

$ \limsup\limits_{t\rightarrow \infty} E\big(G^{-1}(t)\big)<M'. $

对任意$ \epsilon_1>0 $, 取$ \delta_1 = \frac{\epsilon_1}{M'} $, 由Chebyshev不等式可得

$ \limsup\limits_{t\rightarrow \infty} P\{G^{-1}(t)>\delta^{-1}_1\}\leq \limsup\limits_{t\rightarrow \infty} \delta_1 E\big(G^{-1}(t)\big)< \delta_1 M' = \epsilon_1, $

即

$ \limsup\limits_{t\rightarrow \infty} P\{G<\delta_1\}<\epsilon_1, $

因此

(3.8) $ \begin{equation} \liminf\limits_{t\rightarrow \infty} P\{G\geq\delta_1\}\geq1-\epsilon_1. \end{equation} $

于是, $ G(t) $为随机持久的.

下面, 我们证明$ I(t) $是随机持久的. 类似地, 定义如下$ V_1(I(t), t) \in C^2({\Bbb R}_+^2, {\Bbb R}_+) $:

$ V_1 = e^{\lambda_1 t} I(t), $

其中$ 0<\lambda_1<d_i $. 由It$ \hat{\rm o} $公式可得

$ {\rm d}V_1 = e^{\lambda_1 t}{\rm d}I+\lambda e^{\lambda_1 t}I{\rm d}t <\sigma e^{\lambda_1 t} {\rm d}t+\alpha_2e^{\lambda_1 t}I {\rm d}B_2. $

与不等式(3.8) 证明类似, 对任意$ \epsilon_2>0 $, 取$ \bar{\delta}_2 = \frac{\sigma}{\lambda_1\epsilon_2} $, 则有

(3.9) $ \begin{equation} \liminf\limits_{t\rightarrow \infty} P\{I(t)\leq\bar{\delta}_2\}\geq 1-\epsilon_2. \end{equation} $

最后我们证明, 对任意$ 0<\epsilon_2<1 $, 存在正常数$ \delta_2 $使得

$ \liminf\limits_{t\rightarrow \infty} P\{I(t)\geq\delta_2\}\geq 1-\epsilon_2. $

定义$ W_1 = \frac{e^t}{I(t)} $. 由It$ \hat{\rm o} $公式及不等式(3.5) 可得

$ \begin{eqnarray*} {\rm d}W_1& = &\frac{e^t}{I(t)}{\rm d}t+e^t{\rm d}I^{-1} \leq e^t\Big(\frac{-\sigma \ell^2}{(\alpha^2+\ell^2)}\frac{1}{I^2}+(d_i+\alpha_2^2+1)\frac{1}{I}\Big){\rm d}t-\frac{\alpha_2e^t}{I}{\rm d}B_2\\ &<&\left(\frac{(d_i+\alpha_2^2+1)^2(\alpha^2+\ell^2)}{4\sigma \ell^2}\right)e^t{\rm d}t-\frac{\alpha_2e^t}{I}{\rm d}B_2. \end{eqnarray*} $

同样地, 对任意$ \epsilon_2>0 $, 取$ \delta_2 = \frac{4\sigma \ell^2\epsilon_2}{(d_i+\alpha_2^2+1)^2(\alpha^2+\ell^2)} $, 则有

(3.10) $ \begin{eqnarray} \liminf\limits_{t\rightarrow \infty} P\{G\geq\delta_2\}\geq 1-\epsilon_2. \end{eqnarray} $

综上所述, 由式(3.7)–(3.10) 可得随机系统(1.2) 是随机持久的.

为了证明系统(1.2) 存在唯一的遍历平稳分布, 我们先给出一些已有的结论.

设$ X(t) $是$ {\Bbb R}^d $中的自治$ Markov $过程, 可表示为如下随机微分方程

$ {\rm d}X(t) = f(X(t)){\rm d}t+\sum\limits_{r = 1}^k g_r(X(t)){\rm d}B_r(t), $

其扩散矩阵为

$ A(X) = a_{ij}(X)\mbox, {\qquad} a_{ij}(X) = \sum\limits_{r = 1}^k g_r^i(X)g_r^j(X). $

引理3.1^[26]  对Markov过程$ X(t) $, 若存在具有正则边界的有界区域$ U \in {\Bbb R}^d $具有如下性质:

$ \rm (P.1) $ 对任意的$ x \in U $, 扩散矩阵$ A(X) $为严格正定的;

$ \rm (P.2) $ 存在非负$ C^2 $ -函数$ V $, 使得在$ {\Bbb R}^d \backslash U $上$ LV $为负数;

则Markov过程$ X(t) $存在唯一的遍历平稳分布$ \mu(\cdot) $.

定理3.3  对任意初始值$ \big(G(0), I(0)\big)\in {\Bbb R}^2_+ $, 系统(1.2) 存在唯一的遍历平稳分布$ \mu(\cdot) $.

证  为了证明定理3.3, 我们需要验证引理3.1中的条件$ \rm (P.1) $、$ \rm (P.2) $. 系统(1.2) 的扩散矩阵为

$ A = \left( \begin{array}{cc} \alpha_1^2G^2{\quad} & 0 \\ 0 {\quad} & \alpha_2^2I^2 \end{array} \right). $

显然满足条件$ \rm (P.1) $. 下面我们验证条件$ \rm (P.2) $. 由定理3.1, 对任意初始值$ \big(G(0), I(0)\big)\in {\Bbb R}^2_+ $, 系统(1.2) 存在唯一的全局正解$ \big(G(t), I(t)\big) $. 定义函数$ V(G, I) \in C^2({\Bbb R}^2_+, {\Bbb R}_+) $:

$ V(G, I) = (G-\ln G-1)+(I-\ln I-1). $

由It$ \hat{\rm o} $公式可得

$ \begin{eqnarray*} LV& = &(1-\frac{1}{G})(G_{in}-S_gG-\frac{S_iGI}{n+I}+\frac{S_k}{r^m+I^m})+(1-\frac{1}{I})(\frac{\sigma G^2}{\alpha^2+G^2}-d_iI)+\frac{1}{2}(\alpha_1^2+\alpha_2^2)\\ &\leq&(G_{in}-S_gG+\frac{S_k}{r^m+I^m}-\frac{G_{in}}{G}+S_g+\frac{S_i I}{n+I})+(\frac{\sigma G^2}{\alpha^2+G^2}-d_iI-\frac{\sigma G^2}{I(\alpha^2+G^2)})\\ &&+d_i+\frac{1}{2}(\alpha_1^2+\alpha_2^2)\\ &\leq&G_{in}-S_gG+\frac{S_k}{r^m}-\frac{G_{in}}{G} +S_g+S_i+\sigma-d_iI-\frac{\sigma G^2}{I(\alpha^2+G^2)} +d_i+\frac{1}{2}(\alpha_1^2+\alpha_2^2)\\ :& = &M_1'-S_gG-\frac{G_{in}}{G}-d_iI-\frac{\sigma G^2}{I(\alpha^2+G^2)}, \end{eqnarray*} $

其中$ M_1' = {G_{in}+\frac{S_k}{r^m}+S_g+S_i+\sigma+d_i+\frac{1}{2}(\alpha_1^2+\alpha_2^2)} $.

考虑如下开集$ U_{\epsilon} = \{(G, I)\in {\Bbb R}^2_+|\epsilon<G<\frac{1}{\epsilon}, \epsilon^3<I<\frac{1}{\epsilon^3}\} $, 其中$ 0<\epsilon<1 $为一个充分小的常数且满足

$ \epsilon<\min\bigg\{\frac{G_{in}}{M_1'}, \frac{S_g}{M_1'}, \frac{\sigma}{\alpha^2M_1'}, \sqrt[3]{\frac{d_i}{M_1'}}\bigg\}. $

将$ U^{C}_{\epsilon} $分为如下四个区域:

$ U_1 = \{(G, I)\in {\Bbb R}^{2}_+|0<G\leq\epsilon\}, \quad U_2 = \{(G, I)\in {\Bbb R}^{2}_+|0<I\leq\epsilon^3, G>\epsilon\}, $

$ U_3 = \{(G, I)\in {\Bbb R}^{2}_+|G\geq\frac{1}{\epsilon}\}, \quad U_4 = \{(G, I)\in {\Bbb R}^{2}_+|I\geq\frac{1}{\epsilon^3}\}. $

于是有

情形1  在区域$ U_1 $上, 有

$ LV\leq M_1'-S_gG-\frac{G_{in}}{G}-d_iI-\frac{\sigma G^2}{I(\alpha^2+G^2)}\leq M_1'-\frac{G_{in}}{G} \leq M_1'-\frac{G_{in}}{\epsilon} <0. $

情形2  在区域$ U_2 $上, 有

$ \begin{eqnarray*} LV&\leq&M_1'-S_gG-\frac{G_{in}}{G}-d_iI-\frac{\sigma G^2}{I(\alpha^2+G^2)} \leq M_1'-\frac{\sigma G^2}{I(\alpha^2+G^2)}\\ &\leq & M_1'-\frac{\sigma \epsilon^2}{I(\alpha^2+\epsilon^2)} \leq M_1'-\frac{\sigma \epsilon^2}{\epsilon^3(\alpha^2+\epsilon^2)} <0. \end{eqnarray*} $

$ \rm 情形 3 $  在区域$ U_3 $上, 有

$ LV\leq M_1'-S_gG-\frac{G_{in}}{G}-d_iI-\frac{\sigma G^2}{I(\alpha^2+G^2)} \leq M_1'-S_gG \leq M_1'-\frac{S_g}{\epsilon} <0. $

$ \rm 情形 4. $  在区域$ U_4 $上, 有

$ LV\leq M_1'-S_gG-\frac{G_{in}}{G}-d_iI-\frac{\sigma G^2}{I(\alpha^2+G^2)} \leq M_1'-d_iI \leq M_1'-\frac{d_i}{\epsilon^3} <0. $

综上所述, 对任意的$ (G, I)\in U^{C}_{\epsilon} $, 有$ LV(G, I)<0 $, 即引理3.1中条件$ (P_2) $成立, 因此系统(1.2) 存在唯一遍历平稳分布$ \mu(\cdot) $.

为了验证上述理论结果, 我们将给出对确定性模型(1.1) 以及随机模型(1.2) 数值模拟的结果. 对于随机系统(1.2), 我们采用Milstein高阶方法^[27], 其对应的离散方程为

(4.1) $ \begin{equation} \left\{ \begin{array}{rl} G(j+1) = &{ } G(j)+\Big(G_{in}-S_gG(j)-\frac{S_iG(j)I(j)}{n+I(j)} +\frac{S_k}{r^m+I^m(j)}\Big)\Delta t+\alpha_1G(j)\sqrt{\Delta t}\xi_j\\ &{ } +\frac{\alpha_1^2}{2}G(j)(\xi_j^2-1)\Delta t, \\ I(j+1) = &{ } I(j)+\Big(\frac{\sigma G^2(j)}{\alpha^2+G^2(j)}+d_iI(j) \Big)\Delta t +\alpha_2I(j)\sqrt{\Delta t}\zeta_j+\frac{\alpha_2^2}{2}I(j)(\zeta_j^2-1)\Delta t, \end{array} \right. \end{equation} $

其中$ \xi_j $和$ \zeta_j $, $ j = 1, 2, \cdots , n $, 为相互独立的高斯随机变量. 本节的主要目的是通过探究确定性(无噪声) 和对应的随机(有噪声) 葡萄糖-胰岛素系统的全局动力学, 进一步研究环境波动如何影响胰岛素对血糖的产生和利用.

对于确定性系统(1.1) 和随机系统(1.2), 取初始值$ (G(0), I(0)) = (100, 40) $并选取适当的参数: $ G_{in} = 2.16 $, $ S_k = 2 $, $ S_g = 5\times 10^{-6} $, $ S_{i} = 0.1 $, $ n = 80 $, $ m = 4 $, $ r = 80 $, $ \sigma = 6.27 $, $ \alpha = 105 $, $ d_{i} = 0.08 $ (详见表 1).

首先, 由定理2.1和定理2.2可知, 系统存在唯一的正平衡点, 且该平衡点是全局渐进稳定的(如图 1). 下面, 改变$ \alpha_1, \alpha_2 $的值来探究噪声强度对系统(1.2) 动力学的影响.

其次, 取定噪声强度$ \alpha_1 = 0.01 $, $ \alpha_2 = 0.01 $. 数值模拟结果(见图 2) 表明, 在一段时间之后, 胰岛素浓度和血糖浓度将在确定性模型的平衡点$ (80.83, 29.17) $附近波动. 图 2还分别给出血糖浓度和胰岛素浓度的密度函数图, 从图中可以明显看出, 它们在$ 80.83 $和$ 29.17 $附近呈正态分布.

提高环境噪声强度至$ \alpha_1 = 0.05 $, $ \alpha_2 = 0.05 $. 由定理3.3, 系统(1.2) 仍存在遍历平稳分布(见图 3). 值得注意的是, 此时血糖浓度和胰岛素浓度在平衡点附近波动的强度更大, 其密度函数图也能很好地反映这一现象. 对比图 2和图 3可知, 随着噪声强度的减少, 随机模型的波动将逐渐减少并趋向于确定性模型的动力学.

本文建立并研究了确定性葡萄糖-胰岛素模型及其相应的随机模型的动力学性质. 我们通过理论分析得到模型(1.1) 存在唯一全局渐进稳定的正平衡点$ E(G^*, I^*) $. 为了研究随机模型(1.2) 的长期动力学行为, 我们首先证明该模型全局正解的存在唯一性. 随后, 我们研究了随机模型(1.2) 的持久性并证明模型存在唯一的遍历平稳分布. 这些结果表明, 血糖浓度和胰岛素浓度是动态变化的, 随机模型的解不会趋于一个稳定的值, 而是会在其确定性模型的解附近波动. 因此尽管无法准确的预测血糖浓度及胰岛素浓度的准确值, 我们仍可以通过平稳分布估计出其密度函数. 最后, 数值模拟的结果也清晰地表明, 随着噪声强度的降低, 随机模型的波动将逐渐减小并趋于其确定性模型的动力学.

本文所做的工作将会为临床医师对糖尿病患者的诊疗及胰岛素治疗糖尿病策略调整提供帮助. 此外, 最近的一些文献[28-35]表明白噪声无法描述种群在自然界中遭受突然的灾难性冲击或外来物种入侵的现象, 而电报噪声(即颜色噪声) 可以描述两个或多个环境状态之间的切换. 因此, 我们将在未来的工作中进一步考虑更现实但复杂的模型, 如带有时滞和电报噪声的葡萄糖-胰岛素模型.