image
KongHung

Opportunity favors only the prepared mind!

免责声明:网站内容仅供个人学习记录,禁做商业用途,转载请注明出处。

版权所有 © 2017-2020 NEUSNCP个人学习笔记 辽ICP备17017855号-2

详解传染病动力学模型


        近期查看一下传染病动力学方面的部分资料,包网络资料和专业书籍,比如知呼上写的传染病动力学数学模型。因为去年年初首先从武汉爆发的新型冠状病毒(COVID-19),网络上出现了大量有关传染病动力学模型的网络资料。查看之后发现大部分均是直接给出相应的数学模型,对相关参数的含义以及涉及到常微分方程的求解等内容,均阐述得不特征细致,让初次接触者容易混淆相关参数,对模型含义理解不清楚,不深刻。为此,我经以整理了常见的传染病SI模型、SIS模型、SIR模型和SIRS模型的相关知识,并附上采用变易系数法求解常微分方程的具体过程,以帮助初学者快速了解相关知识。
        根据个体状态,传染病动力学模型将总人数为 \(N\) 的人群分为:\(S(Susceptible)\)易感者,\(I(Infected)\)染病者和\(R(Recovered)\)恢复者三类个体。
  • \(\beta_0\) —— 单位时间内,一位易感者接触一位染病者被传染的概率,即传染概率(infection proability)。
  • \(\gamma\) —— 单位时间内,一位染病者经过恢复成为恢复者的概率。
  • \(\delta\)  —— 单位时间内,一位恢复者失去免疫力的概率。
  • \(n\) —— 单位时间内,一位染病者接触到的人数,也称为接触率(cotract rate)。
  • \(S_{(t)},I_{(t)},R_{(t)}\) —— 分别表示 \(t\) 时刻易感者,染病者和恢复者的人数。
        一位易感者与一位染病者接触,就会以概率 \(\beta_0\) 被传染从而变为染病者。我们把带有传染概率 \(\beta_0\) 的接触称为有效接触,此时的接触率称为有效接触率(adequate contact rate),即 \(\lambda=\beta_0 n\) 为单位时间内一位染病者能够传染的人数,表示一个染病者传染他人的能力,反映了一个染病者的活动能力和环境条件等因素。在与染病者的接触者中,若接触对象为易感者,则接触者可能被传染;若拉解对象为非易感者,则接触者不会发生传染。易知,易感者在总人口所占的比例为 \(\frac{S_{(t)}}{N}\) ,因此,每一位染病者平均对易感者的有效接触率为 \(\beta_0 n \frac{S_{(t)}}{N}\) ,也是每一位染病者平均对易感者的传染率,简称为传染率。从而 \(t\) 时刻单位时间内,所有染病者能够传染易感者的人数(即单位时间内易感者减少人数,也是染病增加人数)为 \(\beta_0 n \frac{S_{(t)}}{N} I_{(t)}\) ,此人数称为疾病的发生率(incidence)。
        在上面的描述中,实际上假设了在总人口数不大时,接触率与环境内的总人口成正比,即 \(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}\]
        在初始时刻,即 \(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}- \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\)时,染病者人数占总人口人数的比例与时间的关系如下图所示:
        从上图可知,在第27或28天时,增加的染病者人数达到最快,约在第40天时,易感者基本上已全部被传染。
        当总人口数量很大时,得其微分方程为:

\[\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\)时,染病者人数占总人口数的比例与时间的关系如下图所示:
        从上图可知,大约在第31天时,增加的染病者人数达到最快,约在第45天时,染病者占总人口数的比例基本上趋于稳定,约为占90%。与SI模型相比,染病者人数增加速度最快的时间和染病者占总人口数的比例稳定时间均有所延后。在染病者有以恢复率为 \(\gamma\) 恢复为易感者的条件下,易感者不会全被染病者传染。
        当总人口数量很大时,得其微分方程为:
\[\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\) ,因此,只需讨论前两个方程
\[\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)}\\\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)}\\\end{cases}\]
由 \(\frac{dS_{(t)}}{dt}<0\) ,则 \(S{(t)}\) 是严格单调递减的且有下界 \(0\) ,同时对于 \(t \geqslant 0\) ,都有 \(0 \leqslant S_{(t)} \leqslant N\) 成立。于是存在极限
\(\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)}\) 达到极大值。从而得到该系统的轨迹线图如下
当初始时刻易感者人数 \(S_{(0)}=S_0>\sigma (S_{(0)}=S_0>\rho)\) 时,随着时间增长,染病者人数 \(I_{(t)}\) 将先增加达到最大值  \(I_{(\sigma)} (I_{(\rho)})\) ,然后再逐渐减少而最终消亡。该现象表明,只要 \(S_0>\sigma (S_0>\rho)\) ,即 \(S_0\frac{\beta_0 n}{\gamma N}>1(S_0\frac{\beta}{\gamma}>1)\),疾病就会流行。令\[R_0=\frac{\beta_0 n}{\gamma N}S_0=\frac{S_0}{\sigma}(R_0=\beta \frac{1}{\gamma}S_0=\frac{S_0}{\rho})\]此时,当  \(R_0>1\) 时,疾病流行,当 \(R_0<1\) 时,疾病不会流行,染病者人数 \(I_{(t)}\) 将单下降而趋向于 \(0\) 。因此,\(R_0=1\) 是区分疾病流行与否的阈值,要防止疾病流行,必须减小 \(R_0\) ,使其小于 \(1\) 。\(\frac{1}{\gamma}\) 表示平均移出时间,也就是平均染病期。由恢复率 \(\frac{1}{\gamma}\) 的定义知,若染病者人数为 \(I_{(t)}\) ,则单位时间内恢复者的数目为 \(\gamma I_{(t)}\) ,所以经过时间 \(\frac{1}{\gamma}\) ,染病者全部恢复。
        
假设 \(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.









 
参考资料: #
最近更新: 2021年3月10日 21:40:34
浏览:2542