根据个体状态,传染病动力学模型将总人数为 \(N\) 的人群分为:\(S(Susceptible)\)易感者,\(I(Infected)\)染病者和\(R(Recovered)\)恢复者三类个体。
- \(\beta_0\) —— 单位时间内,一位易感者接触一位染病者被传染的概率,即传染概率(infection proability)。
- \(\gamma\) —— 单位时间内,一位染病者经过恢复成为恢复者的概率。
- \(\delta\) —— 单位时间内,一位恢复者失去免疫力的概率。
- \(n\) —— 单位时间内,一位染病者接触到的人数,也称为接触率(cotract rate)。
- \(S_{(t)},I_{(t)},R_{(t)}\) —— 分别表示 \(t\) 时刻易感者,染病者和恢复者的人数。
在上面的描述中,实际上假设了在总人口数不大时,接触率与环境内的总人口成正比,即 \(n=aN\) ,此时的有效接触率为 \(\beta_0 a N\) ,令 \(\beta_0 n =\beta_0 a N=\beta N\) ,其中 \(\beta=\frac{\beta_0 n}{N}\) 表示有效接触率在总人口中 \(N\) 所占的比例,称为有效接触系数或传染系数。从而可得 \(t\) 时刻单位时间内,新增染病者人数也就是疾病的发生率为 \(\beta_0 n \frac{S_{(t)}}{N} I_{(t)} =\beta S_{(t)} I_{(t)}\) 。
当总人口数量很大时,假设接触率与总人口数成正比显然违背实际情况,因为在单位时间内,一位染病者能够接触到他人的数目是有限的。因此,此时一般假设接触率 \(n\) 为一常数 ,即有效接触率为 \(\lambda=\beta_0 n\) ,从而得到疾病的发生率为 \(\lambda \frac{S_{(t)}}{N} I_{(t)}\) ,称为标准发生率(standard incidence)。以上相关的更多解释见文献[1]。
1. SI 模型
SI模型总共包含两种状态的个体,即易感者和染病者。其模型框图为

\[\begin{cases} \frac{dS_{(t)}}{dt}= -\beta_0 n \frac{S_{(t)}}{N}I_{(t)} \\ \frac{dI_{(t)}}{dt}= \beta_0 n \frac{S_{(t)}}{N}I_{(t)} \end{cases} \qquad 或 \qquad \begin{cases} \frac{dS_{(t)}}{dt}=-\beta S_{(t)}I_{(t)}\\ \frac{dI_{(t)}}{dt}=\beta S_{(t)}I_{(t)} \end{cases}\]
\[\frac{dI_{(t)}}{dt}- \beta_0 n I_{(t)}+\frac{\beta_0 n}{N}I_{(t)}^2=0 \qquad 或 \qquad \frac{dI_{(t)}}{dt}- \beta N I_{(t)}+\beta I_{(t)}^2\]上述方程两端同时除以 \(I_{(t)}^2\) ,并令 \(y=\frac{1}{I_{(t)}}\) ,进而得以下方程
\[\frac{dy}{dt}=-\beta_0 n y+\frac{\beta_0 n}{N} \qquad 或 \qquad \frac{dy}{dt}=-\beta N y+\beta\]此微分方程的具体解法见附录 \(A\) 。结合初始条件,进而可得
\[I_{(t)}=\frac{NI_0}{I_0+(N-I_0)e^{- \beta_0 nt}} \qquad 或 \qquad I_{(t)}=\frac{NI_0}{I_0+(N-I_0)e^{- \beta Nt}}\]因为 \(\underset{t\to+\infty}{lim} e^{-t}=0\) ,因此 \(\underset{t\to+\infty}{lim} I_{(t)}=N\),进而可得 \(\underset{t\to+\infty}{lim} S_{(t)}=0\) 。此外,当 \(I_{(t)}=\frac{N}{2}\) 时,单位时间内染病者人数增加的最快,时间\[t_m=\frac{1}{\beta_0 n}ln(\frac{N}{I_0}-1) \qquad 或 \qquad t_m=\frac{1}{\beta N}ln(\frac{N}{I_0}-1)\]由上式可知,当初始染病人数不变条件下,接触率和传染概率越大,\(t_m\) 越小,传播的速度也就越快。
假设 \(N =10^7,I_0=10,n=5,\beta_0=0.1\)时,染病者人数占总人口人数的比例与时间的关系如下图所示:
当总人口数量很大时,得其微分方程为:
\[\begin{cases} \frac{dS_{(t)}}{dt}= -\lambda \frac{S_{(t)}}{N}I_{(t)} \\ \frac{dI_{(t)}}{dt}= \lambda \frac{S_{(t)}}{N}I_{(t)} \end{cases}\]类似上述求解,可得
\[I_{(t)}=\frac{NI_0}{I_0+(N-I_0)e^{-\lambda t}}\]和
\[t_m=\frac{1}{\lambda}ln(\frac{N}{I_0}-1)\]
2. SIS 模型
SIS模型总共包含两种状态的个体,即易感者和染病者。其模型框图为

\[\begin{cases} \frac{dS_{(t)}}{dt}= -\beta_0 n \frac{S_{(t)}}{N}I_{(t)} + \gamma I_{(t)}\\ \frac{dI_{(t)}}{dt}= \beta_0 n \frac{S_{(t)}}{N}I_{(t)} - \gamma I_{(t)}\end{cases} \qquad 或 \qquad \begin{cases} \frac{dS_{(t)}}{dt}=-\beta S_{(t)}I_{(t)} +\gamma I_{(t)}\\ \frac{dI_{(t)}}{dt}=\beta S_{(t)}I_{(t)}-\gamma I_{(t)} \end{cases}\] 在初始时刻,即 \(t=0\) 时,有\(S_{(0)}=S_0, I_{(0)}=I_0\) ,且对于\(t\geqslant 0\),\(S_{(t)}+I_{(t)}=N\) 均成立。因此,\(S_{(t)}=N-I_{(t)}\) 代入第二个微分方程,于是
\[\frac{dI_{(t)}}{dt}+(\gamma- \beta_0 n) I_{(t)}+\frac{\beta_0 n}{N}I_{(t)}^2=0 \qquad 或 \qquad \frac{dI_{(t)}}{dt}+(\gamma- \beta N) I_{(t)}+\beta I_{(t)}^2=0\]上述方程两端同时除以 \(I_{(t)}^2\) ,并令 \(y=\frac{1}{I_{(t)}}\) ,进而得以下方程
\[\frac{dy}{dt}=(\gamma-\beta_0 n) y+\frac{\beta_0 n}{N} \qquad 或 \qquad \frac{dy}{dt}=(\gamma-\beta N) y+\beta\]此微分方程的具体解法见附录 \(A\) 。结合初始条件,进而可得
\[I_{(t)}=\frac{(\gamma-\beta_0 n)N }{\beta_0 n}/\left((\frac{(\gamma-\beta_0 n)N}{I_0\beta_0 n}+1)e^{(\gamma- \beta_0 n) t}-1\right)\]或
\[I_{(t)}=\frac{(\gamma-\beta N)}{\beta}/\left((\frac{(\gamma-\beta N)}{I_0\beta}+1)e^{(\gamma- \beta N) t}-1\right)\] 假设 \(N =10^7,I_0=10,n=5,\beta_0=0.1,\gamma=0.05\)时,染病者人数占总人口数的比例与时间的关系如下图所示:
当总人口数量很大时,得其微分方程为:
\[\begin{cases} \frac{dS_{(t)}}{dt}= -\lambda \frac{S_{(t)}}{N}I_{(t)}+\gamma I_{(t)} \\ \frac{dI_{(t)}}{dt}= \lambda \frac{S_{(t)}}{N}I_{(t)}-\gamma I_{(t)} \end{cases}\]类似上述求解,可得
\[I_{(t)}=\frac{(\gamma-\lambda) N}{\lambda}/\left((\frac{(\gamma-\lambda )N}{\lambda I_0}+1)e^{(\gamma-\lambda)t}-1 \right)\]
3. SIR 模型
SIR模型总共包含三种状态的个体,即易感者、染病者和恢复者。其模型框图为

\[\begin{cases} \frac{dS_{(t)}}{dt}= -\beta_0 n \frac{S_{(t)}}{N}I_{(t)} \\ \frac{dI_{(t)}}{dt}= \beta_0 n \frac{S_{(t)}}{N}I_{(t)} - \gamma I_{(t)}\\\frac{dR_{(t)}}{dt}=\gamma I_{(t)}\end{cases} \qquad 或 \qquad \begin{cases} \frac{dS_{(t)}}{dt}= -\beta S_{(t)}I_{(t)} \\ \frac{dI_{(t)}}{dt}= \beta S_{(t)}I_{(t)} - \gamma I_{(t)}\\\frac{dR_{(t)}}{dt}=\gamma I_{(t)}\end{cases}\] 在初始时刻,即 \(t=0\) 时,有\(S_{(0)}=S_0, I_{(0)}=I_0,R_{(0)}=R_0\) ,且对于\(t\geqslant 0\),\(S_{(t)}+I_{(t)}+R_{(t)}=N\) 均成立。
上述微分方程没有解析解,下面通过其动力系统分析以得到它解信息。上述方程两端相加得
\[\frac{dS_{(t)}}{dt}+\frac{dI_{(t)}}{dt}+\frac{dR_{(t)}}{dt}=0\]两侧同时积分得\(S_{(t)}+I_{(t)}+R_{(t)}=C(常数)\)。
由于微分方程的前两个式子不含 \(R\) ,因此,只需讨论前两个方程
\(\underset{t\to+\infty}{lim} S_{(t)}=S_\infty\)。
由上述微分方程的第二个式子除以第一个式子得
\[\frac{dI_{(t)}}{dS_{(t)}}=\frac{\sigma}{S_{(t)}}-1,\sigma=\frac{\gamma N}{\beta_0 n} \qquad 或\qquad \frac{dI_{(t)}}{dS_{(t)}}=\frac{\rho}{S_{(t)}}-1,\rho=\frac{\gamma}{\beta}\]当 \(S_{(t)}=\sigma(S_{(t)}=\rho)\) 时,\(I_{(t)}\) 达到极大值。从而得到该系统的轨迹线图如下

假设 \(N =10^7,I_0=10,n=5,\beta_0=0.2,\gamma=0.5\)时,染病者人数、易感者人数和恢复者人数点总人口数的比例与时间的关系如下图所示:

\[\begin{cases} \frac{dS_{(t)}}{dt}= -\lambda \frac{S_{(t)}}{N}I_{(t)} \\ \frac{dI_{(t)}}{dt}= \lambda \frac{S_{(t)}}{N}I_{(t)} - \gamma I_{(t)}\\\frac{dR_{(t)}}{dt}=\gamma I_{(t)}\end{cases}\]类似上述求解可得
\[R_0=\frac{\lambda}{\gamma N}S_0=\frac{S_0}{\tau}(\tau=\frac{\gamma N}{\lambda})\]4. SIRS 模型
SIRS模型总共包含三种状态的个体,即易感者、染病者和恢复者。其模型框图为

\[\begin{cases} \frac{dS_{(t)}}{dt}= \delta R_{(t)}-\beta_0 n \frac{S_{(t)}}{N}I_{(t)} \\ \frac{dI_{(t)}}{dt}= \beta_0 n \frac{S_{(t)}}{N}I_{(t)} - \gamma I_{(t)}\\\frac{dR_{(t)}}{dt}=\gamma I_{(t)}-\delta R_{(t)} \end{cases} \qquad 或 \qquad \begin{cases} \frac{dS_{(t)}}{dt}=\delta R_{(t)} -\beta S_{(t)}I_{(t)} \\ \frac{dI_{(t)}}{dt}= \beta S_{(t)}I_{(t)} - \gamma I_{(t)}\\\frac{dR_{(t)}}{dt}=\gamma I_{(t)}-\delta R_{(t)}\end{cases}\] 在初始时刻,即 \(t=0\) 时,有\(S_{(0)}=S_0, I_{(0)}=I_0,R_{(0)}=R_0\) ,且对于\(t\geqslant 0\),\(S_{(t)}+I_{(t)}+R_{(t)}=N\) 均成立。
上述微分方程没有解析解,下面通过其动力系统分析以得到它解信息。令 \(R_{(t)}=N-S_{(t)}-I_{(t)}\)
\[\begin{cases} \frac{dS_{(t)}}{dt}= \delta(N-S_{(t)}-I_{(t)})-\beta_0 n \frac{S_{(t)}}{N}I_{(t)} \\ \frac{dI_{(t)}}{dt}= \beta_0 n \frac{S_{(t)}}{N}I_{(t)} - \gamma I_{(t)}\\\frac{dR_{(t)}}{dt}=\gamma I_{(t)}-\delta (N-S_{(t)}-I_{(t)}) \end{cases} \qquad 或 \qquad \begin{cases} \frac{dS_{(t)}}{dt}=\delta (N-S_{(t)}-I_{(t)}) -\beta S_{(t)}I_{(t)} \\ \frac{dI_{(t)}}{dt}= \beta S_{(t)}I_{(t)} - \gamma I_{(t)}\\\frac{dR_{(t)}}{dt}=\gamma I_{(t)}-\delta (N-S_{(t)}-I_{(t)})\end{cases}\]此时 \(S_{(t)}、I_{(t)}、R_{(t)}\) 无解析解,只可求得其数值解。
假设 \(N =10^7,I_0=10,n=5,\beta_0=0.2,\gamma=0.2,\delta=0.3\)时,染病者人数、易感者人数和恢复者人数点总人口数的比例与时间的关系如下图所示:

\[\begin{cases} \frac{dS_{(t)}}{dt}= \delta R_{(t)}-\lambda \frac{S_{(t)}}{N}I_{(t)} \\ \frac{dI_{(t)}}{dt}= \lambda \frac{S_{(t)}}{N}I_{(t)} - \gamma I_{(t)}\\\frac{dR_{(t)}}{dt}=\gamma I_{(t)}-\delta R_{(t)} \end{cases}\]令 \(R_{(t)}=N-S_{(t)}-I_{(t)}\) ,则可得
\[\begin{cases} \frac{dS_{(t)}}{dt}= \delta (N-S_{(t)}-I_{(t)})-\lambda \frac{S_{(t)}}{N}I_{(t)} \\ \frac{dI_{(t)}}{dt}= \lambda \frac{S_{(t)}}{N}I_{(t)} - \gamma I_{(t)}\\\frac{dR_{(t)}}{dt}=\gamma I_{(t)}-\delta (N-S_{(t)}-I_{(t)}) \end{cases}\]此时 \(S_{(t)}、I_{(t)}、R_{(t)}\) 无解析解,只可求得其数值解。
附录A一阶线性微分方程 (常数变易法)
\[\frac{dy}{dx}=P(x)y+Q(x) \qquad (a_1)\]其中\(P(x),Q(x)\)在考虑的区间上是\(x\)的连续函数。若\(Q(x)=0\),式\((a_1)\)变为
\[\frac{dy}{dx}=P(x)y \qquad (a_2)\]
称式\((a_2)\)为一阶齐次线性微分方程。若\(Q(x)\neq 0\),则称\((a_1)\)为一阶非齐次线性微分方程。式\((a_2)\)为变量可分离方程,其通解为
\[y=Ce^{\int P(x)dx} \qquad (a_3)\]此处\(C\)为任意常数。
下面讨论一阶非齐次线性微分方程通解的求法。以上不难看出,式\((a_2)\)为式\((a_1)\)的特殊情形,对式\((a_2)\)形式求解,可得到其通解为
\[y=C(x)e^{\int P(x)dx} \qquad (a_4)\]此为式\((a_3)\)中将常数\(C\)变为\(x\)的待定函数\(C(x)\)所对应的形式。
对式\((a_4)\)两边同时对\(x\)求导,得
\[\frac{dy}{dx}=\frac{dC(x)}{dx}e^{\int P(x)dx}+C(x)P(x)e^{\int P(x)dx} \qquad (a_5)\]将式\((a_4)和(a_5)\)代入式\((a_1)\),得到
\[\frac{dC(x)}{dx}=Q(x)e^{-\int P(x)dx} \qquad (a_6)\]两边对\(x\)积分得
\[C(x)=\int Q(x)e^{-\int P(x)dx}dx+ C \qquad (a_7)\] 将式\((a_7)\)代入式\((a_4)\),得到式\((a_1)\)的通解为
\[y=e^{\int P(x)dx}\left(\int Q(x)e^{-\int P(x)dx}dx+ C\right)\]此外,式\((a_1)\)也可以通过积分因子法求解其通解,但相对于常数变易法略烦琐,常数变易法和积分因子法的详细过程见参考文献[2]。
参考文献
[1] 马知恩,周义仓,王稳地,靳祯. 传染病动力学的数学建模与研究. 科学出版社,2004.
[2] 王高雄,周之铭,朱思铭,王寿松. 常微分方程(第四版). 高等教育出版社,2020.
