传染率为βSI/(cN + I)的SIRS模型及数值模拟

A SIRS Model with Infection Rate of βSI/(cN + I) and Its Numerical Simulation

摘要: 本文应用经典的SIRS传染病动力学模型，提出一种传染率为的新的SIRS模型，采用数值模拟的方法，研究相应模型解的稳定性问题。通过大量数值试验，分析了重要参数对模型解的稳定性影响。

Abstract: In this paper, a new SIRS model with the infection rate of is proposed by using the classical SIRS epidemic dynamics model. The stability of the solution of the corresponding model is studied by numerical simulation. The influence of important parameters on the stability of the model solution is analyzed through a large number of numerical tests.

1. 引言

在传染病的数学模型中，最早Kermack和Mckendrick在1927年提出了一个著名的SIR“仓室”模型。所谓SIR仓室模型就是针对某类传染病将该地区的人群分成三个仓室 $S\left(t\right),I\left(t\right),R\left(t\right)$。 $S\left(t\right)$ 表示t时刻未染病但有可能被该类传染病传染的人数， $I\left(t\right)$ 表示t时刻已被感染成病人而且具有传染率的人数， $R\left(t\right)$ 表示t时刻已从染病者类移出的人数。而且他们假设此环境的总人口保持一个常数，假设一个病人一旦与易感者接触就必然具有一定的传染力，而且一个病人能感染的易感者数目与此环境内的易感者总数成正比。很显然这种假设太简单，后来Capasso和Serio在对1973年巴里(意大利东南部港市)发生的霍

乱这一地方病传播的研究后引入一种饱和的传染率 $g\left(I\right)S$，其中， $g\left(I\right)=\frac{kI}{1+\alpha I}$ 这种传染率的重要性在于染病者和易感者的有效接触数量将随着人群的拥挤作用或易感者采取保护措施而出现饱和、甚至下降。

文献 [1] 考虑具有饱和发生率 $\frac{\beta S{I}^p}{1+\alpha I^q}$ 的传染病模型，文献 [2] 采用了当 $p=q=1$ 的情况，研究了人类免

疫缺陷病毒1型和丙型肝炎病毒，文献 [3] 讨论了当 $p=q=1$ 时的SIRS型传染病模型，并分析了这类模型的分支情况，文献 [4] - [14] 同样也研究了饱和发生率的传染病模型。

基于此，本文将研究具有代表意义的具有常数输入且疾病发生率为非线性的动力系统模型，以期为进一步理解种群动力学和流行病动力学的机理提供支撑，并提出一种新传染率下的SIRS模型，对模型应用数值模拟方法，研究解的稳定性，以及重要参数对稳定性的影响。

2. 传染率为 $\frac{\beta SI}{cN+I}$ 的新的SIRS模型

2.1. 基本模型与假设

将所研究的种群分为易感者、染病者和恢复者三类，以 $S=S\left(t\right)$ 、 $I=I\left(t\right)$ 和 $R=R\left(t\right)$ 分别表示t时刻易感者、染病者和恢复者的数量，建立模型的基本假设为：

1) $bN$ 表示对种群的常数输入率，且均为易感染者；

2) $f(N)$ 表示种群的自然死亡率，与种群数量N有关，满足 $f^{\prime }\left(N\right)>0$， $0<f\left(0\right)<b<f\left(+\infty \right)$ ；当种群中无传染病存在时，种群(N)的变化率符合方程 $N^{\prime }=\left[b-f\left(N\right)\right]N$ ；

3) $\alpha$ 表示染病者的因病死亡率；

4) $\gamma$ 表示染病者的恢复率；

5) $\epsilon$ 表示恢复者的免疫失去率，即恢复者中一部分人将失去免疫力，从而又成为易感者；

6) 通过有效接触传染病传播的传染率为 $\frac{\beta SI}{cN+I}$。

用框图表示上述传染病假设的传播规律如下：

根据上述仓式框图建立模型：

$\{\begin{array}{l}\frac{\text{d}S}{\text{d}t}=bN-f\left(N\right)S-\frac{\beta SI}{cN+I}+\epsilon R\triangleq P\left(S,I,R\right)\\ \frac{\text{d}I}{\text{d}t}=\frac{\beta SI}{cN+I}-\left(f\left(N\right)+\alpha +\gamma \right)I\triangleq Q\left(S,I,R\right)\\ \frac{\text{d}R}{\text{d}t}=\gamma I-\left(f\left(N\right)+\epsilon \right)R\triangleq W\left(S,I,R\right)\end{array}$ (2-1)

2.2. 模型的分析

由于 $N\left(t\right)=S\left(t\right)+I\left(t\right)+R\left(t\right)$，所以有，

$N^{\prime }=\left(b-f\left(N\right)\right)N-\alpha I.$

做变量替换：

$x=\frac{S}{N},y=\frac{I}{N},z=\frac{R}{N}$,(2-2)

则 $x\left(t\right),y\left(t\right),z\left(t\right)$ 分别表示t时刻易感者、感染者和治愈者在种群中多占的比例。直接计算可得：

$\{\begin{array}{l}{x}^{\prime }=b-bx-\left(\frac{\beta }{c+y}-\alpha \right)xy+\epsilon z\\ y^{\prime }=y\left[\frac{\beta }{c+y}+\alpha y-\left(b+\alpha +\gamma \right)\right]\\ z^{\prime }=\gamma y-\left(b+\epsilon \right)z+\alpha yz\end{array}$ (2-3)

和

$N^{\prime }=\left[b-\alpha y-f\left(N\right)\right]N$ (2-4)

根据模型的实际意义，设模型(2-3)和(2-4)的初始条件为

$x\left(0\right)>0,y\left(0\right)\ge 0,z\left(0\right)\ge 0,N\left(0\right)>0$

由 $N=S+I+R$ 有 $x+y+z=1$。将 $x=1-y-z$ 代入系统(2-3)的第二个方程可得：

$\{\begin{array}{l}{y}^{\prime }=y\left[\frac{\beta \left(1-y-z\right)}{c+y}+\alpha y-\left(b+\alpha +\gamma \right)\right]\triangleq yF\left(y,z\right)\\ z^{\prime }=\gamma y-\left(b+\epsilon \right)z+\alpha yz\triangleq G\left(y,z\right)\end{array}$ (2-5)

显然， $y=0$ 是模型(2-5)的解。有解的存在惟一性知：当初值 $y\left(0\right)>0$ 时，对于 $t>0$ 都有 $y\left(t\right)>0$，进而由模型(2-3)的第一个和第三个方程有 $z\left(t\right)>0$， $x\left(t\right)>0$，所以模型(2-5)有正不变集

$\Omega =\left\{\left(y,z\right):y\ge 0,z\ge 0,y+z\le 1\right\}$

2.3. 模型求解的Matlab实现

记控制传染病的阈值 $R_0=\frac{\beta }{c\left(b+\alpha +\gamma \right)}$，当 $R_0<1$ 该模型(2-5)有无病平衡点 $P_0\left(y_0,z_0\right)=\left(0,0\right)$，且在该 $P_0$ 点是全局渐近稳定的，模型的解应是趋向于无病平衡点的。

给定新模型的初值和参数如下：

初值： $t=0$ 时， $S\left(t\right)=80\text{/}万人，I\left(t\right)=10\text{/}万人，R\left(t\right)=10\text{/}万人$

参数： $b=0.01193$， $f\left(N\right)=0.00614+0.0001\times N$， $\epsilon =0.97$， $\alpha =0.00905$， $\gamma =0.07$， $c=0.02$。

有 $R_0$ 与 $\beta$ 的关系不难求出当 $\beta <c\left(b+\alpha +\gamma \right)=0.0018$ 时，模型的解在P₀点是全局渐近稳定的；当 $\beta >\text{0}\text{.0018}$ 时，该模型的解趋向于另一个点P₁ (此处不做具体研究)，但是此时在P₀点不是全局渐进稳定的。

下面对当 $\beta =\text{0}\text{.0018}$ 时运行本文中编写的4阶R-K [15] 程序，并用所得到的数据画图，见图1和图2：

Figure 1. Time-varying image of susceptible S

图1. 易感者S随时间的变化图像

Figure 2. Time-varying image of susceptible patient I and healer R

图2. 染病者I和治愈者R随时间的变化图像

从上面的图像中不难看出，随时间变化，染病者和恢复者占总人口数的比例趋向于零，而易感染者所占比例趋向于1，即此时，模型的解趋向于无病平衡点。

增大 $\beta$ 的值使 $R_0>1$，不妨取 $\beta =0.02$，运行程序，同样用得到的数据作图，见图3和图4：

Figure 3. Time-varying image of susceptible S

图3. 易感者S随时间的变化图像

Figure 4. Time-varying image of susceptible patient I and healer R

图4. 染病者I和治愈者R随时间的变化图像

从图像中以看出，此时模型的解并不趋向于点 $P_0\left(S_0,I_0,R_0\right)=\left(\frac{K}{d},0,0\right)$，即模型的解在无病平衡点是不稳定的；但是从图像中可以看出模型的解趋向于另一个点 $P_1$，该点可能就是地方病平衡点。

2.4. 结论

从上述两次数值模拟中可以看到，如果取 $\beta \le 0.0018$，则感染者和治愈者所占比例均趋向于零，而易感染者比例趋向于1，故传染病没有蔓延；但是当 $\beta >0.0018$，感染者和治愈者所占比例均不趋向于零，故传染病可能会蔓延。

所以，在本模型给定上述参数不变的情况下，控制 $\beta <0.0018$ 即可有效的控制传染病，避免传染病的蔓延。

3. 总结

本文主要研究的是在传染病模型中其他参数不变的前提下通过控制传染病的阈值 $R_0<1$ 确定一个重要参数 $\beta$ 的范围，从而验证无病平衡点的全局渐近稳定性。

另外考虑增大和减小 $\beta$ 的值，观察传染病模型解的趋向变化。

新的SIRS传染病模型是在已有的经典传染病模型基础上的改进，同样由 $R_0<1$ 给出 $\beta$ 的一个范围，使新的传染病模型的解在无病平衡点是全局稳定的，并且从图像中观察传染病模型解的趋向；另外给出一个使 $R_0>1$ 的 $\beta$ 值，同样作图，通过对比两次数值实验，可以看出两个 $\beta$ 值使模型的解呈现两种趋向，从而比较传染病是否蔓延的两种状态。