微分方程建模实例(2)

传染病动力学是用数学模型研究某种传染病再某一地区是否蔓延下去,成为当地的”地方病”,或最终该病将被消除。下面以Kermack和Mckendrick提出的阈值模型为例说明传染病动力学数学模型的建模过程。
STEP01 模型假设
(1)被研究人群是封闭的,总人数为n。s(t),i(t)和r(t)分别表示t时刻人群中易感染者、感染者(病人)和免疫者的人数。起始条件为s0个易感染者,i0个感染者,n-s0-i0个免疫者。
(2)易感人数的变化率与当时的易感人数和感染人数之积成正比,比例系数为λ。
(3)免疫者人数的变化率与当时的感染者人数成正比,比例系数为μ。
(4)三类人总的变化率代数和为零。
STEP02 模型建立
$$
\begin{cases}
\frac{ds}{dt} = -\lambda si \\
\frac{di}{dt} = \lambda si - \mu i \\
\frac{dr}{dt} = \mu i \\
s(t)+i(t)+r(t)=n
\end{cases}
$$
以上模型被称为Kermack-Mckendrick模型。
STEP03 模型分析和求解
对于Kermack-Mckendrick模型方程组,无法求出s(t)、i(t)和r(t)的解析解,于是转到平面s-i上来讨论解的性质。
将方程组中前两个方程消去dt,可得
$$
\begin{cases}
\frac{di}{ds} = \frac{1}{\sigma s}-1 \\
i\vert_{s=s_{0}} = i_{0}
\end{cases}
$$
其中σ=λ/μ,是一个传染期内每个患者有效接触得平均人数,称为接触数。用分离变量法可求出它的解为
$$
i=(s_{0}+i_{0})-s+ \frac{1}{\sigma }ln\frac{s}{s_{0}}
$$
当s0小于等于1/σ时,传染病不会蔓延。患者人数一直在减少并逐渐消失。当s0大于1/σ时,患者人数会增加,传染病开始蔓延,健康者的人数在减少。当s(t)减少至1/σ时,患者在人群中的比例达到最大值,然后患者数逐渐减少到零。由此可知,1/σ是一个阈值,要想控制传染病的流行,应控制s0使之小于此阈值。

传染病控制两个途径
一是提高卫生和医疗水平。卫生水平越高,传染性接触率λ就越小;医疗水平越高,恢复系数μ就越大。这样,阈值1/σ就越大。
二是降低s0。由s0+i0+r0=n可知,要想减少s0可采取提高r0的方法,而这需要通过预防接种和群体免疫等措施来实现。

STEP04 参数估计

参数σ的值可以由实际数据估计来得到,记s,i分别是传染病流行结束后的健康者人数和患者人数。
当流行结束后,患者都将转化为免疫者,所以i=0,可得
$$i_{\infty}=0=s_{0}+i_{0}+r_{0}+\frac{1}{\sigma }ln\frac{s_{\infty}}{s_{0}}$$
解出σ得
$$\sigma = \frac{lns_{0}-lns_{\infty }}{s_{0}+i_{0}-s_{\infty }}$$
例:以1950年上海市某全托幼儿所发生得一起水痘流行过程为例,应用K-M模型进行模拟。该所儿童总人数n为196人,既往患过水痘而此次未感染者40人,查不出水痘患病史而本次流行期间感染水痘者96人,既往无明确水痘史,本次又未感染得幸免者60人。
解:以初始值s0=155(196-40-1),s=60代入公式得σ=0.0099

1
2
3
4
5
6
7
8
9
10
import numpy as np
from scipy.integrate import odeint
s0=155.0
i0=1.0
s_inf=60.0
sigma=(np.log(s0)-np.log(s_inf))/(s0+i0-s_inf)
print("sigma=",sigma)
S=np.array([155, 153, 139, 101])
I=(s0+i0)-S+1/sigma*np.log(S/s0)
print("所求的解为:\n",I)
0%