SEIR模型及多染病仓室再生数的推导

SEIR模型

模型推导

在许多传染病中，易感者被感染后到有症状可以传播之前，存在一个暴露期。我们设平均的暴露期为$\frac{1}{\kappa }$，暴露类为$E$，结合易感染类$S$, 染病类$I$， 恢复类$R$和总人口规模，得到以下流程：

所以模型如下：
$$\{\begin{array}{rl}{S}^{\prime }& =-\beta SI\\ E^{\prime }& =\beta SI-\kappa E\\ I^{\prime }& =\kappa E-\alpha I\end{array}\text{\hspace{1em}(1.1)}$$
事实上，有些疾病在暴露期也存在传染性，这可以通过传染因子$\epsilon$来降低传染性假设的模拟，模型更新为：
$$\{\begin{array}{rl}{S}^{\prime }& =-\beta S(I+\epsilon E)\\ E^{\prime }& =\beta S(I+\epsilon E)-\kappa E\\ I^{\prime }& =\kappa E-\alpha I\end{array}\text{\hspace{1em}(1.2)}$$
以及初始条件
$$S(0)=S_0,E(0)=E_0,I(0)=I_0$$

再生数

下一代矩阵

为了求得再生数，我们需要引入“下一代矩阵”的概念。
假设存在$n$个疾病仓室和$m$个无病仓室，令$x\in R^n$表示各个疾病仓室的人数, $y\in R^m$表示各个无病仓室的人数。记${\mathcal{F}}_i$表示第$i$个疾病仓室中发生继发性感染的增加率，${\mathcal{V}}_i={\mathcal{V}}_i^{-}-{\mathcal{V}}_i^{+}$表示第$i$个疾病仓室的疾病进展。${\mathcal{V}}_i^{-}$表示第i个仓室的移除率，${\mathcal{V}}_i^{+}$表示其他方式转入到i仓室的变化率${}^{[1]}$。于是仓室模型可以写成以下形式：
$$\{\begin{array}{rl}{x}_i^{\prime }& ={\mathcal{F}}_i(x,y)-{\mathcal{V}}_i(x,y)\\ y_j^{\prime }& =g_j(x,y)\end{array}\text{\hspace{1em}(2.1)}$$
从上述定义和模型中，可以发现蕴含以下假设${}^{[2]}$：

1. ${\mathcal{F}}_i(0,y)=0,{\mathcal{V}}_i(0,y)=0$对所有的$y⩾0$和$i=1,2,...,n$
2. 无病系统$y^{\prime }=g(0,y)$有唯一渐进稳定平衡点，即具有形如(0,y)的初始条件的所有解当$t\to \infty$时都趋于点$(0,y_0)$。称此点为无病平衡点。
3. ${\mathcal{F}}_i(x,y)⩾0,{\mathcal{V}}_i^{-}(x,y)⩾0,{\mathcal{V}}_i^{+}(x,y)⩾0$,对所有非负$x,y$以及$i=1,2,...,n$.
4. ${\mathcal{V}}_i^{-}(x,y)=0$,当$x_i=0,i=1,2,...,n$
5. ${\sum }_{i=1}^n{\mathcal{V}}_i(x,y)⩾0$, 对所有非负$x$和$y$.

假设单个染病者进入原来没有疾病的入口，通过人员传播疾病的最初能力由上述模型(2.1)关于无病平衡点$(0,y_0)$的线性化的研究决定。从假设1中可以得到
$$\frac{\partial {\mathcal{F}}_i}{\partial y_j}(0,y_0)=\frac{\partial {\mathcal{V}}_i}{\partial y_j}(0,y_0)=0$$
$$D{\mathcal{F}}_i(x,y)=(\frac{\partial {\mathcal{F}}_i(x,y)}{\partial x_1}...\frac{\partial {\mathcal{F}}_i(x,y)}{\partial x_n}\frac{\partial {\mathcal{F}}_i(x,y)}{\partial y_1}...\frac{\partial {\mathcal{F}}_i(x,y)}{\partial y_m})(d{x}_1...d{x}_{n}d{y}_1...d{y}_m{)}^T$$
$$D{\mathcal{F}}_i(x,y)=\sum _{j=1}^n\frac{\partial {\mathcal{F}}_i(x,y)}{\partial x_j}d{x}_j+\sum _{j=1}^m\frac{\partial {\mathcal{F}}_i(x,y)}{\partial y_j}d{y}_j$$
$$D{\mathcal{F}}_i(0,y_0)=\sum _{j=1}^n\frac{\partial {\mathcal{F}}_i(0,y_0)}{\partial x_j}d{x}_j$$
$$D{\mathcal{F}}_i(0,y_0)=(\frac{\partial {\mathcal{F}}_i(0,y_0)}{\partial x_1}...\frac{\partial {\mathcal{F}}_i(0,y_0)}{\partial x_n})(d{x}_1...d{x}_n{)}^T$$
由$\frac{\partial {\mathcal{F}}_i(0,y_0)}{\partial x_j}$为常数，所以
$${\mathcal{F}}_i(0,y_0)=(\frac{\partial {\mathcal{F}}_i(0,y_0)}{\partial x_1}...\frac{\partial {\mathcal{F}}_i(0,y_0)}{\partial x_n})(x_1...x_n{)}^T$$
则方程可以写为
$$x^{\prime }=(F-V)x\text{\hspace{1em}(2.2)}$$
其中
$$F=(\frac{\partial {\mathcal{F}}_i(0,y_0)}{\partial x_j}{)}_{n\times n}和V=(\frac{\partial {\mathcal{V}}_i(0,y_0)}{\partial x_j}{)}_{n\times n}$$
由假设2，系统(2.1)的线性稳定型完全由$F-V$的线性稳定性决定。
由初始条件$x(0)=x_0,F=0$(没有继发性感染)，解得
$$x(t)=e^{-Vt}{x}_0$$
所以每个仓室内指标个案所经历的期望时间为
$${\int }_0^{\infty }t{e}^{-Vt}{x}_{0}dt=V^{-1}{x}_0$$
由指标个案产生的继发性染病的期望数可表示成患病期的期望时间与出现继发性染病率的积，$F$的$(i,j)$元素是由仓室$j$中的指标个案再仓库$i$中产生的继发性染病率，$V^{-1}$可解释成最初进入疾病仓室$j$的个体在疾病仓室$i$所经历的期望时间。所以由$F{V}^{-1}{x}_0$给出，定义$K_L=F{V}^{-1}$为系统在无病平衡点的下一代矩阵，称为有小定义域的下一代矩阵。

由下一代矩阵定义再生数

定义${\mathcal{R}}_0=\rho (F{V}^{-1})$为$K_L$的谱半径，如果${\mathcal{R}}_0<1$无病平衡点渐近稳定，${\mathcal{R}}_0>1$不稳定，则${\mathcal{R}}_0$为再生数。

下面证明${\mathcal{R}}_0$为再生数。
定义1 $Z$符号型矩阵 \qquad 如果$A=sI-B,B⩾0$，则$A$称为有$Z$符号型的。
定义2 $M$-型矩阵 $A=sI-B,B⩾0,s⩾\rho (B)$，则称$A$为$M$-矩阵。
引理1${}^{[3]}$如果$A$有$Z$符号型，则$A^{-1}⩾0$，当且仅当$A$是非奇异$M$-矩阵。
由假设知$F$非负，$V$非对角线元素非正，因此$V$有$Z$符号型。而$V$的列元素之和为正或为零，则$V$为$M$-矩阵，不妨设为非奇异的，由引理1知$V^{-1}⩾0$。因此$K_L=F{V}^{-1}$也是非负的。
引理2 \qquad 如果$F$非负，$V$是非奇异$M$-矩阵，则${\mathcal{R}}_0=\rho (F{V}^{-1})<1$，当且仅当$(F-V)$的所有特征值有负实部。
证明 \qquad 由引理1，$V^{-1}⩾0$，因此$(I-F{V}^{-1})$有$Z$符号型。又由引理1，$(I-F{V}^{-1}{)}^{-1}⩾0$当且仅当$\rho (F{V}^{-1})<1$。由等式
$$\{\begin{array}{l}(V-F{)}^{-1}=V^{-1}(I-F{V}^{-1})\\ V(V-F{)}^{-1}=I+F(V-F{)}^{-1}\end{array}$$
推得
$$\begin{array}{rl}V(V-F{)}^{-1}& =I+F(V-F{)}^{-1}\\ & =I+F{V}^{-1}(I-F{V}^{-1}{)}^{-1}\\ & =(I-F{V}^{-1})(I-F{V}^{-1}{)}^{-1}+F{V}^{-1}(I-F{V}^{-1}{)}^{-1}\\ & =(I-F{V}^{-1}{)}^{-1}\\ F(V-F{)}^{-1}& =F{V}^{-1}(I-F{V}^{-1}{)}^{-1}\end{array}$$
注意到$F⩾0,F{V}^{-1}⩾0$得到$(V-F{)}^{-1}⩾0$当且仅当$(I-F{V}^{-1}{)}^{-1}⩾0$。
而$(V-F)$有$Z$符号型，所以$(V-F{)}^{-1}⩾0$当且仅当$(V-F)$是非奇异$M$-矩阵。由于非奇异$M$-矩阵特征值实部皆正，则主对角线元素都为正(反方向可由递归法化成上三角阵证明)，非主对角线上元素都为负，则该矩阵为非奇异$M$-矩阵${}^{[4]}$。从而$(F-V)$的所有特征值有负实部。
定理 \qquad 对于模型(1)，如果${\mathcal{R}}_0<1$，则
无病平衡点局部渐近稳定，如果${\mathcal{R}}_0>1$，则
不稳定。
证明 \qquad 对系统按照上述方式求线性化的雅可比矩阵，得分块结构
$$J=\left(\begin{array}{cc}F-V& 0\\ J_{21}& J_{22}\end{array}\right)$$
若$J$的所有特征值具有负实部，则无病平衡点局部渐近稳定。显然$J$的特征值为$F-V$和$J_{22}$的特征值。由假设2知$J_{22}$的所有特征值有负实部。由引理2知$F-V$的所有特征值有负实部，当且仅当$\rho (F{V}^{-1})<1$。所以${\mathcal{R}}_0=\rho \left(F{V}^{-1}\right)<1$时，无病平衡点局部渐近稳定。

对于${\mathcal{R}}_0⩾1$的不稳定性可以由连续性建立。若${\mathcal{R}}_0⩽1$，则对$\forall \epsilon >0,((1+\epsilon )I-F{V}^{-1})$是非奇异$M$-矩阵，由引理1，$((1+\epsilon )I-F{V}^{-1}{)}^{-1}⩾0$。由引理2，$(F-(1+\epsilon )V)$的所有特征值具有负实部。因为$\epsilon >0$任意，又特征值是矩阵元素的连续函数，则$(F-V)$的所有特征值具有负实部。反之，假设$(F-V)$所有的特征值具有负实部，对任何整数$\epsilon ，(V+\epsilon I-F)$为非负$M$矩阵，由引理2，$\rho (F(V+\epsilon I{)}^{-1})<1$。由$\epsilon$的任意性，可得$\rho (F{V}^{-1})⩽1$，因此$(F-V)$至少有一个特征值具有正实部，当且仅当$\rho (F{V}^{-1})>1$，所以当${\mathcal{R}}_0>1$时，无病平衡点不稳定。

综上，${\mathcal{R}}_0=\rho (F{V}^{-1})$为再生数

SEIR模型的再生数

考虑模型(1.2)，疾病状态为$E$和$I$，得到
$$\mathcal{F}=\left(\begin{array}{c}\epsilon E\beta N+I\beta N\\ 0\end{array}\right)\mathcal{V}=\left(\begin{array}{cc}\kappa E& 0\\ -\kappa E& \alpha I\end{array}\right)$$
则
$$F=\left(\begin{array}{cc}\epsilon \beta N& \beta N\\ 0& 0\end{array}\right)V=\left(\begin{array}{cc}\kappa & 0\\ -\kappa & \alpha \end{array}\right)V^{-1}=\left(\begin{array}{cc}\frac{1}{\kappa }& 0\\ & \\ \frac{1}{\alpha }& \frac{1}{\alpha }\end{array}\right)$$
可以计算
$$K_L=F{V}^{-1}=\left(\begin{array}{cc}\frac{\epsilon \beta N}{\kappa }+\frac{\beta N}{\alpha }& \frac{\beta N}{\alpha }\\ 0& 0\end{array}\right)$$
于是${\mathcal{R}}_0=\rho (F{V}^{-1})=\frac{\epsilon \beta N}{\kappa }+\frac{\beta N}{\alpha }$

最后规模关系

对模型(1.2)的第三个式子积分，得到
$$\kappa {\int }_0^{\infty }E(s)ds=\alpha {\int }_0^{\infty }I(s)ds-I_0$$
对模型(1.2)的第一个式子从0~$\infty$积分，得
$$\begin{array}{rl}ln\frac{S_0}{S_{\infty }}& ={\int }_0^{\infty }\beta [I(s)+\epsilon E(s)]ds\\ & =\beta (\epsilon +\frac{\kappa }{\alpha }){\int }_0^{\infty }E(s)ds\\ & ={\mathcal{R}}_0(1-\frac{S_{\infty }}{N})-\frac{\epsilon \beta I_0}{\kappa }\end{array}$$
我们假设$I_0=0$，并且初始染病者在第一阶段就可以传播疾病。则最后规模关系有形式$ln\frac{S_0}{S_{\infty }}={\mathcal{R}}_0(1-\frac{S_{\infty }}{N})$，与简单SIR模型相同。

参考文献

[1]李霞. SEIR传染病模型综述[J]. 北京师范大学本科论文全文,2014-05-15.
[2]Wendi Wang,Xiao-Qiang Zhao. Threshold Dynamics for Compartmental Epidemic Models in Periodic Environments[J]. Journal of Dynamics and Differential Equations,2008,20(3).
[3]Berman A, Plemmons R J.Nonnegative Matrices in the Mathematical Sciences[M].New York:Academic press, 1979.
[4]张家驹.M矩阵的一些性质[J].数学年刊A辑(中文版),1980(01):47-50.