一类HCV模型的随机稳定性和分岔分析

Stochastic Dynamics Analysis of HCV Virus Model

摘要: 在HCV模型中加入高斯白噪声，建立带有随机激励的非线性微分方程组，应用随机降阶和随机平均法相关定理将原方程化为二维的伊藤微分方程。然后，本文应用最大Lyapunov指数和奇异边界理论分析了该随机系统平凡解的局部和全局稳定性，得到了相应的稳定性条件；利用Hamilton理论中的随机平均法，对HCV系统的动态分岔和唯象分岔行为做了研究。最后，通过数值模拟得到HCV流行病在噪声的影响下传染速率会有明显的变化。

Abstract: Gauss white noise was added into the nonlinear HCV model, nonlinear differential equations were established, and the original equation was reduced to a two-dimensional Ito ̂ differential equation by applying the theorem of stochastic central manifold and the correlation theorem of stochastic average method. Then, the maximum Lyapunov index and singular boundary theory are applied to analyze the local stochastic stability and global stochastic stability of the stochastic system respectively, and the stability conditions are obtained. The dynamic bifurcation and phenomenological bifurcation behavior of systems are studied by using stochastic average method for quasi-integrable Hamiltonian systems. Finally, the numerical simulation shows that the infection rate of HCV epidemic will change obviously under the influence of noise.

1. 引言

丙型肝炎病毒(HCV)是引起肝脏疾病的主要因素，会导致肝炎疾病的发生。估计这种慢性疾病的全球流行率已达到7100万。约60%~80%被HCV感染的患者会进展为慢性阶段，而其余患者在急性感染后可自行恢复。约15%至30%的慢性病发展为肝硬化和肝细胞癌。HCV入侵宿主后，会进入血小板，它的一些RNA会残留在血小板。此外，它可以引起血小板的自身免疫反应。文献 [1] [2] 中的模型能够解释HCV感染过程中T淋巴细胞(CTL)和免疫应答的联系。在文献 [3] 中，通过对复杂的HCV模型研究得到了树突状细胞(dc)和CTL的作用。另外，在病毒感染过程中通过对非细胞溶解治疗，HCV的自发清除也被纳入。最近，在文献 [4] 中研究了一个在无限空间传播延迟的非局部反应扩散HCV模型，以用来探究高迁移率基团盒1 (HMGB1)在阻断病毒离子方面的作用。文献 [5] 研究了乙型肝炎病毒(HBV)和乙型肝炎表面抗原(HBsAg)之间抗体对其的阻断作用。文献 [6] 研究了病毒和细胞传播和细胞介导免疫反应对于一般病毒动力学模型的影响。

文献 [7] 提出了一类传染病数学模型，其主要研究了病毒到细胞和细胞到细胞传播以及抗体应答在感染肝细胞非细胞溶解治疗中的影响。分析了该模型，讨论了确定性环境下细胞间传播和治愈率的影响。然而，由于病毒的随机行为和免疫系统的复杂性，确定性模型并不能提供对病毒动力学的充分理解。由于温度、湿度等环境参数的随机波动，在文献 [8] [9] 中提出并分析了几种基于随机方法的传染病流行模型。为了在模型分析中引入随机性，通常使用不同的方法，如参数摄动和连续时间马尔可夫链，以及加入相关的噪声。如今对于流行病的随机分析 [10] 的研究仍然很少，尤其对于丙型肝炎病毒的研究更少。

本文工作在于对HCV动力学的随机建模。因此，参考先前发表的HCV动力学确定性模型 [7]，在此基础上加入白噪声对系统进行随机激励分析。首先，将带有激励的HCV模型利用降阶和平均法理论进行降阶；然后，对其平凡解的稳定性进行研究。最后，模拟分岔图验证理论的有效性。

2. 模型介绍

常见的流行病模型有SIR模型和HIV模型，除此之外，HCV模型的动力学行为也十分丰富。基于确定性的HCV模型：

$\{\begin{array}{l}\frac{\text{d}T}{\text{d}t}=\lambda -{\beta }_{1}TV-{\beta }_{2}TI-d_{1}T+aI\\ \frac{\text{d}I}{\text{d}t}={\beta }_{1}TV+{\beta }_{2}TI-d_{2}T-aI\\ \frac{\text{d}V}{\text{d}t}=kI-d_{3}V-\rho VZ\\ \frac{\text{d}Z}{\text{d}t}=cVZ-d_{4}V\end{array}$ (1)

方程中 $T,I,V$ 和Z分别代表未感染肝细胞、感染肝细胞、病毒和b细胞的密度。这里假设 $\lambda$ 是未受感染的肝细胞固定的繁殖速率， $d_1$ 是未受感染的肝细胞的自然死亡速率。假设 ${\beta }_1$ 是未感染肝细胞被病毒感染的速率， ${\beta }_2$ 是被感染肝细胞感染的速率。感染的肝细胞通过自然死亡以 $d_2$ 的速率被清除，并通过非细胞溶解过程以a的速率转化为未感染的肝细胞。在血液中游离的病毒在感染的肝细胞中产生的速率和衰变的速率分别为k和 $d_3$。病毒进入宿主细胞后，抗体以c的速率刺激自身产生 细胞，使病毒以p的速率中和。最后， $d_4$ 是b细胞的自然死亡率。以上所有参数均为正数。易求得系统的非负平衡点为

$E_1\left(0,\frac{\lambda }{d_1},0,0\right)$， $E_2\left(I^{\ast },T^{\ast },V^{\ast },Z^{\ast }\right)$，其中：

$I^{\ast }=\frac{k\lambda {\beta }_1+\lambda {\beta }_2{d}_3-d_1{d}_2{d}_3-a{d}_1{d}_3}{d_2\left(k{\beta }_1+{\beta }_2{d}_3\right)}$, $T^{\ast }=\frac{d_3\left(a+d_2\right)}{k{\beta }_1+{\beta }_2{d}_3}$, $V^{\ast }=\frac{k\left(k\lambda {\beta }_1+\lambda {\beta }_2{d}_3-d_1{d}_2{d}_3-a{d}_1{d}_3\right)}{d_2{d}_3\left(k{\beta }_1+{\beta }_2{d}_3\right)}$

确定性系统(1)在受到细胞环境、医疗条件以及个人身体素质等随机因素的影响下，则其中的参数 $d_1,d_2,d_3,d_4,k,p,c$ 会变得不稳定。考虑到随机激励给系统带来的影响和参数的实际意义以及运算的复杂程度，我们只考虑参数 $d_1,d_2,k,c$ 受白噪声的影响，对其加入随机项： $d_1\to d_1-\delta \xi \left(t\right)$, $d_2\to d_2-\delta \xi \left(t\right)$, $k\to k-\delta \xi \left(t\right)$, $c\to c-\delta \xi \left(t\right)$，则得到一个新的随机非线性微分方程：

$\{\begin{array}{l}\stackrel{\dot{}}{T}=\lambda -{\beta }_{1}TV-{\beta }_{2}TI-d_{1}T+aI+\delta \left(T-T^{\ast }\right)\xi \left(t\right)\\ \stackrel{\dot{}}{I}={\beta }_{1}TV+{\beta }_{2}TI-d_{2}I-aI+\delta \left(I-I^{\ast }\right)\xi \left(t\right)\\ \stackrel{\dot{}}{V}=kI-d_{3}V-\rho VZ+\delta \left(V-V^{\ast }\right)\xi \left(t\right)\\ \stackrel{\dot{}}{Z}=cVZ-d_{4}V+\delta \left(Z-Z^{\ast }\right)\xi \left(t\right)\end{array}$ (2)

其中 $\delta$ 代表白噪声强度， $\xi \left(t\right)$ 表示零均值白噪声，即 $E\left(\xi \left(t\right)\right)=0$, $E\left(\xi \left(t\right)\xi \left(t+\tau \right)\right)=\delta \left(\tau \right)$，令 $X=T\left(t\right)-T^{\ast }$, $Y=I\left(t\right)-I^{\ast }$, $U=V\left(t\right)-V^{\ast }$, $W=Z\left(t\right)-Z^{\ast }$，代入系统(2)可得

$\{\begin{array}{l}\stackrel{\dot{}}{X}=AX-BY-{\beta }_1{T}^{\ast }U-{\beta }_{1}XU-{\beta }_{2}XY+\delta X\xi \left(t\right)\\ \stackrel{\dot{}}{Y}=CX+DY+{\beta }_1{T}^{\ast }U+{\beta }_{1}XU+{\beta }_{2}XY+\delta Y\xi \left(t\right)\\ \stackrel{\dot{}}{U}=EY-FU-\rho UY+\delta U\xi \left(t\right)\\ \stackrel{\dot{}}{W}=HW+c{Z}^{\ast }U+cWU+\delta W\xi \left(t\right)\end{array}$ (3)

其中：

$A=-\left({\beta }_1{V}^{\ast }+{\beta }_2{I}^{\ast }+d_1\right)$, $B={\beta }_2{T}^{\ast }-a$, $C={\beta }_1{V}^{\ast }+{\beta }_2{I}^{\ast }$, $D={\beta }_2{T}^{\ast }-d_2-a$, $E=k-\rho V^{\ast }$, $F=d_3+\rho I^{\ast }$, $H=d_4-c{V}^{\ast }$

3. 随机HCV模型的初步处理

由随机稳定性的等价定义 [11] 知，系统(3)在平衡点 $E_0\left(0,0,0,0\right)$ 处的稳定性等价于系统(2)在平衡点 $E_2\left(I^{\ast },T^{\ast },V^{\ast },Z^{\ast }\right)$ 处的稳定性，则系统(3)的特征方程可表示如下：

$f=\left(\lambda +H\right)\left({\lambda }^3+a{\lambda }^2+b\lambda +c\right)=0$ (4)

其中：

$a=F-D-A$, $b=AD-E{\beta }_1{T}^{\ast }-AF-BC-DF$, $c=AE{\beta }_1{T}^{\ast }+CE{\beta }_1{T}^{\ast }+ADF-BCF$

若令 $b_0=-\frac{c}{a},b=b_0+\epsilon \frac{2\left[{\left(a+2\epsilon \right)}^2-\frac{c}{a}\right]}{a+2\epsilon }$，其中 $\epsilon$ 为足够小的参数，则通过计算可得特征方程(4)的特征值为： ${\lambda }_1=-H$, ${\lambda }_2=-a-2\epsilon$, ${\lambda }_{3,4}=\epsilon \pm \sqrt{\frac{c}{a+2\epsilon }-\frac{{\epsilon }^2}{4}}i$，当 $\epsilon =0$ 时， $b=-\frac{c}{a}$，则特征值为： ${\lambda }_1=-H$, ${\lambda }_2=-a$, ${\lambda }_{3,4}=\pm \sqrt{b}i$ 是系统(3)的退化临界焦点，那么当 $\epsilon$ 足够小时，系统(3)具有不变流形：

$W_{loc}^c\left(O\right)=\left\{\left(x_1,x_2,x_3,x_4\right)|x_1=h\left(x_2,x_3\right),x_4=g\left(x_2,x_3\right),\left|x_2\right|+\left|x_3\right|\le 1\right\}$,

则(3)系统可近似代替为。

$\{\begin{array}{l}{\stackrel{\dot{}}{x}}_2=-\sqrt{b}{x}_3+\varphi \left(x_2,x_3\right)+\delta \left(l_{12}{x}_2+l_{13}{x}_3\right)\xi \left(t\right)\\ {\stackrel{\dot{}}{x}}_3=\sqrt{b}{x}_2+\phi \left(x_2,x_3\right)+\delta \left(m_{12}{x}_2+m_{13}{x}_3\right)\xi \left(t\right)\end{array}$ (5)

其中：

$\begin{array}{c}\varphi \left(x_2,x_3\right)=\left(l_2{h}_1+l_8{g}_1\right)x_2^3+\left(l_2{h}_2+l_3{h}_1+l_8{g}_2+l_9{g}_1\right)x_2^2{x}_3\\ +\left(l_2{h}_4+l_3{h}_2+l_8{g}_4+l_9{g}_2\right)x_2{x}_3^2\\ +l_4{x}_2^2+l_5{x}_2{x}_3+l_6{x}_3^2+\left(l_3{h}_4+l_9{g}_4\right)x_3^3\end{array}$

$\begin{array}{c}\phi \left(x_2,x_3\right)=\left(m_2{h}_1+m_8{g}_1\right)x_2^3+\left(m_2{h}_2+m_3{h}_1+m_8{g}_2+m_9{g}_1\right)x_2^2{x}_3\\ +\left(m_2{h}_4+m_3{h}_2+m_8{g}_4+m_9{g}_2\right)x_2{x}_3^2\\ +m_4{x}_2^2+m_5{x}_2{x}_3+m_6{x}_3^2+\left(m_3{h}_4+m_9{g}_4\right)x_3^3\end{array}$

在文献 [12] 中，罗和郭将极坐标变换 $x=r\cos \theta$ 和 $y=r\sin \theta$ 与伊藤公式相结合，并用随机平均法将系统(3)改写为伊藤随机微分方程：

$\{\begin{array}{l}\text{d}r=\left[\left({\mu }_1+\frac{1}{16}{\mu }_2\right)r+\frac{1}{8}{\mu }_3{r}^3\right]\text{d}t+{\left(\frac{1}{8}{\mu }_4{r}^2\right)}^{\frac{1}{2}}\text{d}{W}_r\left(t\right)\\ \text{d}\theta =\left[{\mu }_5+\frac{1}{8}{\mu }_6{r}^2\right]\text{d}t+{\left(\frac{1}{8}{\mu }_2\right)}^{\frac{1}{2}}\text{d}{W}_{\theta }\left(t\right)\end{array}$ (6)

$W_r\left(t\right)$ 和 $W_{\theta }\left(t\right)$ 是独立的标准Wiener过程，并有以下标记：

${\mu }_1=0$，

${\mu }_2={\delta }^2{l}_{12}^2+{\delta }^2{m}_{12}^2+{\delta }^2{l}_{13}^2+3{\delta }^2{m}_{13}^2$，

${\mu }_3=3\left(l_2{h}_1+l_8{g}_1\right)+\left(l_2{h}_4+l_3{h}_2+l_8{g}_4+l_9{g}_2\right)+\left(m_2{h}_2+m_3{h}_2+m_8{g}_2+m_9{g}_1\right)+\left(m_3{h}_4+m_9{g}_4\right)$，

${\mu }_4=3{\delta }^2{l}_{12}^2+{\delta }^2{m}_{12}^2+{\delta }^2{l}_{13}^2+{\delta }^2{m}_{13}^2$，

${\mu }_5=4\sqrt{b}+{\delta }^2\left(l_{12}{l}_{13}-m_{12}{m}_{13}\right)$，

${\mu }_6=3\left(m_2{h}_1+m_8{g}_1\right)-\left(m_2{h}_4+m_3{h}_2+m_8{g}_4+m_9{g}_2\right)-\left(l_2{h}_2+l_3{h}_2+l_8{g}_2+l_9{g}_1\right)-\left(l_4{h}_3+m_9{g}_4\right)$

因此由以上可得，其振幅为：

$\text{d}r=\left[\frac{1}{16}{\mu }_{2}r+\frac{1}{8}{\mu }_3{r}^3\right]\text{d}t+{\left(\frac{1}{8}{\mu }_4{r}^2\right)}^{\frac{1}{2}}\text{d}W\left(t\right)$ (7)

4. 随机稳定性分析

4.1. 局部随机稳定性

讨论线性伊藤随机微分方程的稳定性，就等于方程(7)中 ${\mu }_3=0$ 时的稳定性，则有：

$\text{d}r=\frac{1}{16}{\mu }_{2}r\text{d}t+{\left(\frac{1}{8}{\mu }_4{r}^2\right)}^{\frac{1}{2}}\text{d}W\left(t\right)$ (8)

根据解线性伊藤微分方程精确解析解的方法，可得

$r\left(t\right)=r\left(0\right)\exp ({\displaystyle {\int }_0^t\left[m\left(0\right)-\frac{{\delta }^2\left(0\right)}{2}\right]\text{d}s}+{\displaystyle {\int }_0^t\delta \left(0\right)\text{d}W(s))}$

其中 $m\left(0\right)=\frac{{\mu }_2}{16}$， $\delta \left(0\right)=\frac{{\mu }_4}{8}$。

根据拟不可积Hamilton系统的原理，线性伊藤随机微分方程Lyapunov指数为 $\Vert r\left(t\right)\Vert ={\left(r\left(t,W\right)\right)}^{\frac{1}{2}}$，则：

$\lambda =\underset{t\to +\infty }{\lim }\frac{1}{t}\ln \left(\Vert Z\left(t,z_0\right)\Vert \right)=\underset{t\to +\infty }{\lim }\frac{1}{t}\ln {\left(r\left(t\right)\right)}^{\frac{1}{2}}=\frac{m\left(0\right)-\frac{{\delta }^2\left(0\right)}{2}}{2}=\frac{{\mu }_2-{\mu }_4}{32}$

由乘积遍历性定理 [11] 知，系统的平凡解以概率1渐进稳定的充份必要条件是：最大Lyapunov指数 $\lambda <0$，因此可以得到系统保持局部随机稳定的条件是： ${\mu }_2-{\mu }_4<0$。

因为乘积遍历性定理的最大Lyapunov指数只适用于判断系统平凡解的局部稳定性，对全局稳定性不能判别，所以我们下面讨论一下系统的全局随机稳定性。

4.2. 全局随机稳定性

在随机激励的影响下，漂移系数 $m\left(r\right)$ 和扩散系数 ${\delta }^2\left(r\right)$ 会变得不稳定，经常会出现奇异边界，而得知这种边界上的性态就等价于得知整个扩散过程的性质。所以根据奇异边界理论的边界类别，对系统的平凡解的全局稳定性做出判别。

分为两种情况：

情形1：当 ${\mu }_3=0$ 时，系统此时变为如下形式：

$\text{d}r=\frac{1}{16}{\mu }_{2}r\text{d}t+{\left(\frac{1}{8}{\mu }_4{r}^2\right)}^{\frac{1}{2}}\text{d}W\left(t\right)$ (9)

当 $r=0$ 时，其为系统的第一类奇异边界；当 $r=+\infty$ 时， $m_r=+\infty$, $r=+\infty$ 为系统的第二类奇异边界。

根据奇异边界理论，可分别计算出在 $r=0$ 边界处的扩散指数、漂移指数和特征标值，即

, , ,

若，则是排斥自然边界；若，则是吸引自然边界；若，则是严格自然边界。

同理，可分别计算出在边界处的扩散指数、漂移指数和特征标值，即

, ,

若，则是吸引自然边界；若，则是排斥自然边界；若，则是严格自然边界。

因此得出系统全局稳定的条件：当，是吸引自然边界，是排斥自然边界，则表

明线性的伊藤随机微分方程的平凡解在概率为1的意义下是稳定的，这也表明系统(2)在含有随机激励时在平衡点处是概率稳定的。

情形2：当时，系统此时变为如下形式：

(10)

当时，其为系统的第一类奇异边界；当时，,为系统的第二类奇异边界。

根据奇异边界理论，可分别计算出在边界处的扩散指数、漂移指数和特征标值，即

, ,

若，则是排斥自然边界；若，则是吸引自然边界；若，则是严格自然边界。

同理，可分别计算出在边界处的扩散指数、漂移指数和特征标值，即

, , ,

若，即则是吸引自然边界；若，即则是排斥自然边界；若，即则是严格自然边界。

当是吸引自然边界，是排斥自然边界，系统所有解曲线都从右边界进入到系统的内部被

左边界吸引，当左边界r趋于0处时，所有平凡解是稳定的，因此得出系统全局稳定的条件：，。

5. 随机Hopf分岔

随机分岔是随机激励系统在系统的某个或某些参数发生微小变化时它的性态也会发生变化的现象，这一部分将应用拟不可积Hamilton系统随机平均法分析系统的随机分岔动力学行为。随机分岔分为两类：动态分岔(D-分岔)与维象分岔(P-分岔)。

5.1. 动态分岔

(1) 当时，系统为：

(11)

当时系统是一个确定性的线性系统，不会发生分岔现象。因此下面讨论的情形，令,，由方程生成的连续动态系统为：

(12)

它是方程(11)以x为初值的唯一强解。当,时，0是的一个不动点，对此不动点，是stratonovich随机微分方程意义下的随机参激，设有界，对所有满足椭圆性条件，保证了最多只有一个平稳概率密度，由振幅的伊藤微分方程得到(11)式对应的FPK方程：

(13)

令得到与方程(13)相应的平稳概率密度：

(14)

其中c为归一化常数，故上述动态系统有不动点和非平凡平稳运动两种平稳状态，为前者的不变测度的密度，(14)式为后者的不变测度的密度，为便于研究动态分岔，需分别计算这两个不变测度的Lyapunov指数。

考虑线性的伊藤微分方程，得到方程(11)的解：

(15)

动态系统关于不变测度的Lyapunov指数可定义如下：

(16)

将(15)式代入方程(16)，这里,得不动点Lyapunov指数如下：

(17)

以(14)式为密度的不变测度，将(15)式代入方程(16)，假定有界，并且可积，得到Lyapunov指数：

(18)

令，于是可得，当时不动点的不变测度是稳定的，当时非平凡状态不变测度

是稳定的，所以是系统(2)一个D分岔点。

化简系统(14)可得以下方程：

(19)

其中c为归一化常数，令，显然当时，。是一个函数,不可积，当时，。是在该区间上取得最大值的一个点。时系统(2)发生D分岔，即是系统发生D分岔的临界条件。当时不存在一个点使得取到极大值，则系统(2)不会发生 分岔。

(2) 当时，令，且，则系统(7)变为：

(20)

当时，经分析可得系统有确定性叉形分岔。，时，只有唯一的不变测度其密度为；，时，有三个不变测度分别为：与两个非平凡平稳状态测度，其密度为：

(21)

(22)

式中，则有(17)式可得关于的Lyapunov指数为：

(23)

根据遍历性理论，由方程(18)，(21)和(23)可计算得出关于两个非平凡平稳状态测度的Lyapunov指数：

(24)

由此可知，在处发生D分岔。

由于(21)式有极值，于是可知在处系统(2)会发生P分岔。

5.2. 维象分岔

这部分考虑线性伊藤微分方程的稳态概率密度，从不变测度极值来分析系统的分岔行为。令

，则系统(7)变为：

(25)

则根据伊藤方程的振幅得到方程(25)的FPK方程如下：

(26)

由初值条件，，其中是扩散过程的转移概率密度，的不变测度是稳态概率密度，得到其退化系统的解如下：

(27)

通过计算可得：

(28)

由Namachivays [11] 理论知，非线性随机动力系统稳态性主要是通过不变测度的极值显示的，能显示出系统最根本的稳态性行为。当噪声时，的极值几乎可以表示确定系统的稳态行为。若在处可取得极大值，则轨线会在处的邻域内停留较长时间，即在概率意义下是稳定的点。

6. 数值模拟

由以上得知，原系统的响应过程为，所以响应过程的联合概率密度函数为。选取参数，根据实际情况给参数赋予一定的值，使得，然后选取为分岔参数。设，下图(图1~3)就是模拟的平稳概率密度函数图和联

合概率密度函数图。当时，图1中的平稳概率密度函数图像是一个单调递减的函数，联合概率密度函数图像是单峰尖状的；当时，图2中平稳概率密度函数图像从单峰值函数变为单调递减函数，从时的联合概率密度函数图可以看出系统未发生分岔；当时，在图3中可以看到，在时发生分岔，联合概率密度函数图像从单峰形状变为火山口形状。

(a)的平稳概率密度图 (b)的联合概率密度图

Figure 1. The stationary probability density graph and the corresponding joint probability density graph when

图1. 当时的平稳概率密度图和相对应的联合概率密度图

(a)的平稳概率密度图 (b)的联合概率密度图

Figure 2. The stationary probability density graph and the corresponding joint probability density graph when

图2. 当时的平稳概率密度图和相对应的联合概率密度图

(a)的平稳概率密度图 (b)的联合概率密度图

Figure 3. The stationary probability density graph and the corresponding joint probability density graph when

图3. 当时的平稳概率密度图和相对应的联合概率密度图

7. 结论

本文主要运用随机激励的Hamilton理论和随机平均法分析了HCV模型的动力学行为。结果表明，HCV系统在受到白噪声影响时会变得不稳定，并且发生随机分岔行为，根据参数的实际意义选取 作为分岔参数，当时，图3中的图3(a)图从单调函数变为单峰函数，图3(b)图像由单峰状变化为火山口状，系统发生了分岔。说明此刻可能不利于病毒的控制，需采取一定的策略，以减缓病毒的感染速率，防止病情的恶化。

基金项目

国家自然科学基金(批准号：61863022)。

NOTES

^*第一作者。