数学建模——传染病模型
百度文库- 让每个人平等地提升自我
传染病模型
摘要
当今社会，人们开始意识到通过定量地研究传染病的传播规律，建立传染病的传播模型，可以为预测和控制传染病提供可靠、足够的信息。本文利用微分方程稳定性理论对传统传染病动力学建模方式进行综述，且针对甲流，SARS等新生传染病模型进行建模和分析。
不同类型的传染病的传播过程有其各自不同的特点，我们不是从医学的角度一一分析各种传染病的传播，而是从一般的传播机理分析建立各种模型，如简单模型，SI模型，SIS模型，SIR模型等。本文中，我们应用传染病动力学模型来描述疾病发展变化的过程和传播规律，运用联立微分方程组体现疫情发展过程中各类人的内在因果联系，并在此基础上建立方程求解算法。然后，通过借助Matlab程序拟合出与实际较为符合的曲线并进行了疫情预测，评估各种控制措施的效果，从而不断完善文中的模型。
本文由简到难、全面地评价了该模型的合理性与实用性，而后对模型和数据也做了较为扼要的分析，进一步改进了模型的不妥之处。同时，在对问题进行较为全面评价的基础上又引入更为全面合理的假设，运用双线性函数模型对卫生部的措施进行了评价并给出建议，做好模型的完善与优化工作。
关键词：传染病模型,简单模型，SI，SIS,SIR，微分方程，Matlab。
一、问题重述
有一种传染病（如SARS、甲型H1N1）正在流行，现在希望建立适当的数学模型，利用已经掌握的一些数据资料对该传染病进行有效地研究，以期对其传播蔓延进行必要的控制，减少人民生命财产的损失。考虑如下的几个问题，建立适当的数学模型，并进行一定的比较分析和评价展望。
1、不考虑环境的限制，设单位时间内感染人数的增长率是常数，建立模型求t 时刻的感染人数。
2、假设单位时间内感染人数的增长率是感染人数的线性函数，最大感染时的增长率为零。建立模型求t时刻的感染人数。
3、假设总人口可分为传染病患者和易感染者，易感染者因与患病者接触而得病，而患病者会因治愈而减少且对该传染病具有很强的免疫功能，建立模型分析t 时刻患病者与易感染者的关系，并对传染情况（如流行趋势，是否最终消灭）进行预测。
二、问题分析
1、这是一个涉及传染病传播情况的实际问题，其中涉及传染病感染人数随时间的变化情况及一些初始资料，可通过建立相应的微分方程模型加以解决。
2、问题表述中已给出了各子问题的一些相应的假设。
3、在实际中，感染人数是离散变量，不具有连续可微性，不利于建立微分方程模型。但由于短时间内改变的是少数人口，这种变化与整体人口相比是微小的。
因此，为了利用数学工具建立微分方程模型，我们还需要一个基本假设：感染人数是时间的连续可微函数。
三、模型假设
模型二和模型三的假设条件：
假设一：在疾病传播期内所考察地区的总人数N不变，即不考虑生死，也不考虑迁移。人群分为易感染者（Susceptible）和已感染者（Infective）两类（取两个词的第一个字母，称之为SI模型），以下简称健康者和病人。时刻t这两类人在总人数中所占比例分别记作s(t)和i(t)。
假设二：每个病人每天有效接触的平均人数是常数，称为日接触率。当病人与健康者接触时，使健康者受感染变为病人。
假设三：模型三在假设一和假设二的基础上进行考虑，然后设病人每天治愈的比例为μ，称为日治愈率。病人治愈后成为仍可被感染的健康者，显然1/μ是这种传染病的平均传染期。
模型四的假设条件：
假设四：总人数N不变。人群分为健康者、病人和病愈免疫的移出者（Removed）三类，称SIR模型。三类人在总数N中占的比例分别记作s(t),i(t)和r(t)。
假设五：病人的日接触率为λ，日治愈率为μ（与SI模型相同），传染期接触为σ=λ/μ。
四、符号说明
t ·······························某一具体时刻
x(t)·····························病人人数λ·······························每天每个病人有效接触的人数N································总人数
s(t)·····························健康者总人数
i(t)·····························病人总人数
······························初始时刻病人的比例
i
····························病人的最大值
t
m
`μ····························日治愈率
1/μ···························平均传染率σ·····························接触率
r(t)···························移出者
s
·····························初始时刻健康者的比例
五、模型的建立与求解
模型1
在这个最简单的模型中，设时刻t 的病人人数x(t)是连续、可微函数，并且每天每个病人有效接触（足以使人致病的接触）的人数为常数λ，考察t 到
病人人数的增加，就有
t t x t x t t x ?=-?+)()()(λ
程有个病人，即得微分方时有再设00x t =
)1()0(,d d 0x x x t
x
==λ
方程（1）的解为
)2()(0t
e x t x λ=
结果表明，随着t 的增加，病人人数x(t)无限增长，这显然是不符合实际的。
建模失败的原因在于：在病人有效接触的人群中，有健康人也有病人，而其中只有健康人才可以被传染为病人，所以在改进的模型中必须区别这两种人。
模型2（SI 模型）
的增加率，即
就是病人数个健康者被感染，于是有，所以每天共为病人数为个健康者变为病人，因天可使根据假设，每个病人每Ni Nsi t i t Ns t Ni t s λλλ)()()()(
)3(
d d Nsi t
i
N
λ=
又因为
)4(1
)()(=+t i t s
，则病人的比例为再记初始时刻0)0(i t =
)5()0(,)1(d d 0
i i i i t
i
=-=λ
方程（5）是Logistic 模型。它的解为
)6(11110t e i λ-???
? ??-+
所示。
和图的图形如图和
21~d d ~)(i t
i
t t i
，这个时刻为到达最大值时
可知，第一，当式及图由m
t i t i i ???
??=d d d d 2/11)6(),5(
)7(11ln 01?
??
?
??-=-i t m λ
这时病人增加的最快，可以认为是医院的门诊量最大的一天，预示着传染病高潮的到来，是医疗卫生部门关注的时刻。
况。
，这显然不符合实际情将被传染，全变为病人即所有人终时到来。第二，当可以推迟传染病高潮的健设施、提高卫生水平以改善保越小卫生水平越高。所，表示该地区的卫生水平成反比，因为日接触率与,1→∞→i t t m λλλ 其原因是模型中没有考虑到病人可以治愈，人群中的健康者只能变成病人，病人不会再变成健康者。
模型3（SIS 模型）
有些传染病如伤风、痢疾等愈合后免疫力很低，可以假定无免疫性，于是病人被治愈后变成健康者，健康者还可以被感染再变成病人，所以这个模型成为SIS 模型。
考虑到这一模型的假设条件，于是有
[]Ni(t)-t (t)i(t)i(t)-t)i(t N μλ?=?+Ns （8）
可得微分方程
i i(0) i -i)-(1==μλi dt
di
0 （9） 定义
μλσ= （10）
其中σ是整个传染期内每个病人有效接触的平均人数，称为接触数。 得到
[])1/-(1--di σλi i dt
= （11）
模型4（SIR 模型）
大多数传染者如天花 流感 肝炎 麻疹等治愈后均有很强的免疫力，所以病愈的人既非健康者（易感染者），也非病人（已感染者），因此他们将被移除传染系统，我们称之为移除者，记为R 类。
SIR 模型是指易感染者被传染后变为感染住，感病者可以被治愈，并会产生免疫力，变为移除者。人员流动图为：S-I-R 。
1.模型构成： 在假设1中显然有：
s(t) + i(t) + r(t) = 1 （12）
对于病愈免疫的移出者的数量应为
r
t
d N
Ni d μ= （13）
不妨设初始时刻的易感染者、染病者、恢复者的比例分别为0s （0s ＞0），0i （0i ＞0），0r =0，则SIR 基础模型用微分方程组表示如下：
di
dt ds
dt dr
dt si i si i λμλμ?=-???=-???=?? （14）
s(t) ， i(t)的求解极度困难，在此我们先做数值计算来预估计s(t) ， i(t)的一般变化规律。
2.数值计算
在方程（3）中设λ=1，μ=，i （0）= ，s （0）=，用MATLAB 软件编程：
function y=ill(t,x)
a=1;b=;
y=[a*x(1)*x(2)-b*x(1);-a*x(1)*x(2)];
ts=0:50;
x0=[,];
[t,x]=ode45('ill',ts,x0);
plot(t,x(:,1),t,x(:,2))
pause
plot(x(:,2),x(:,1))
输出的简明计算结果列入表1。i(t) , s(t)的图形以下两个图形，i~s图形称为相轨线,初值i(0)=,s(0)=相当于图2中的P0点,随着t的增,(s,i)沿轨线自右向左运动.由表1、图1、图2可以看出,i(t)由初值增长至约t=7时达到最大值,然后减少,t→∞,i→0,s(t)则单调减少,t→∞,s→. 并分析i(t),s(t)的一般变化规律.
t 0 1 2 3 4 5 6 7 8 i(t)
s(t)
t 9 10 15 20 25 30 35 40 45 i(t) 0 s(t)
3.相轨线分析
我们在数值计算和图形观察的基础上，利用相轨线讨论解i（t）,s（t）的
性质。
D = ｛（s ，i ）| s≥0，i≥0 ， s + i ≤1｝ （15） 在方程（14）中消去t d 并注意到σ的定义，可得
11i s
d d ??
=- ???
s σ 00|s s i i == （16） 所以：
11i s d d ??=- ???s σ ?00i 11s i s i s d d ??
=- ???
??s σ
利用积分特性容易求出方程(5)的解为：
000
1
()ln
s
i s i s s σ
=+-=
（17）
在定义域D 内,(17)式表示的曲线即为相轨线,如图3所示.其中箭头表示了随着时间t 的增加s(t)和i(t)的变化趋向。
图3
下面根据（14），（17）式和图3分析s(t),i(t)和r(t)的变化情况(t→∞时它们的极限值分别记作s ∞, i ∞和r ∞).
1. 不论初始条件s0,i0如何,病人将消失,即：0,t ?→?∞?→?
i 2. 最终未被感染的健康者的比例是 ,在(7)式中令i=0得到, 是方
000
1
ln
0s s i s s σ
∞
∞+-+
= 在(0,1/σ)内的根.在图形上 是相轨线与s 轴在(0,1/σ)内交点的横坐标 3.若0s >1/σ,则开始有
11i s d o d ??=-> ???s σ，i(t)先增加, 令11i s d d ??=- ???
s σ=0，可得当s=1/σ时,i(t)达到最大值：
0001
1ln )m i s i s σσ
=+-+（
然后s<1/σ时，有
11i s d o d ??=-< ???
s σ ，所以i(t)减小且趋于零,s(t)则单调减小至s ∞,如图3中由P1(0s ，0i )出发的轨线
4.若0s ≤1/σ,则恒有
110i s d d ??=-< ???
s σ，i(t)单调减小至零,s(t)单调减小至s ∞,如图3中由P2(s0,i0)出发的轨线
可以看出,如果仅当病人比例i(t)有一段增长的时期才认为传染病在蔓延,那么1/σ是一个阈值,当0s >1/σ(即σ>1/s0)时传染病就会蔓延.而减小传染期接触数σ,即提高阈值1/σ使得0s ≤1/σ(即σ ≤1/0s ),传染病就不会蔓延(健康者比例的初始值0s 是一定的,通常可认为0s 接近1)。
并且,即使0s >1/σ, σ减小时, s ∞增加(通过作图分析), m i 降低,也控制了蔓延的程度.我们注意到在σ=λμ中,人们的卫生水平越高,日接触率λ越小;医疗水平越高,日治愈率μ越大,于是σ越小,所以提高卫生水平和医疗水平有
助于控制传染病的蔓延.
从另一方面看, 1/s s σλμ=?是传染期内一个病人传染的健康者的平均数,称为交换数,其含义是一病人被s σ个健康者交换.所以当 01/s σ≤ 即01s σ≤时必有 .既然交换数不超过1,病人比例i(t)绝不会增加,传染病不会蔓延。
5. 群体免疫和预防：
根据对SIR 模型的分析,当01/s σ≤ 时传染病不会蔓延.所以为制止蔓延,除了提高卫生和医疗水平,使阈值1/σ变大以外,另一个途径是降低0s ,这可以通过比如预防接种使群体免疫的办法做到。
忽略病人比例的初始值0i 有001s r =-,于是传染病不会蔓延的条件01/s σ≤ 可以表为
01
1r σ
≥-
这就是说，只要通过群体免疫使初始时刻的移出者比例（即免疫比例）就可以制止传染病的蔓延。
这种办法生效的前提条件是免疫者要均匀分布在全体人口中，实际上这是很难做到的。据估计当时印度等国天花传染病的接触数 σ=5，至少要有80%的人接受免疫才行。据世界卫生组织报告，即使花费大量资金提高0r ，也因很难做到免疫者的均匀分布，使得天花直到1977年才在全世界根除。而有些传染病的σ更高，根除就更加困难。
6.模型验证:
上世纪初在印度孟买发生的一次瘟疫中几乎所有病人都死亡了。死亡相当于
移出传染系统，有关部门记录了每天移出者的人数，即有了r
t
d d 的实际数据，Kermack 等人用这组数据对SIR 模型作了验证。 首先，由方程（12），（14）可以得到
s r t d d
si si s d dt
λσμσ=-=-=- 1
s r d d s
σ?=-t 上式两边同时乘以d 可 ，两边积分得
0001s
r s r
s r d d s σ==-??0ln |s
s s r σ?=-0r s e s σ-?= 所以：
()0()r t s t s e σ-= (8)
再0(1)(1)r r
t
d i r s r s
e d σμμμ-?
==--=-- (9) 当 1/r σ≤ 时，取（13）式右端r e σ-Taylor 展开式的前3项得：
22000(1)2
r
t s r d r s s r d σμσ=--+- (10)
在初始值0r =0 下解高阶常微分方程得:
0201()(1)()2t r t s th s αμσα?σ??
=
-+-????
(11) 其中222000(1)2s s i ασσ=-+，01
s th σ?α
-=
从而容易由（10）式得出：
2220
2()
2
r t d t d s ch αμ
αμσ?=-
然后取定参数 s0, σ等，画出（11）式的图形，如图4中的曲线，实际数据在图中用圆点表示，可以看出，理论曲线与实际数据吻合得相当不错。
六、模型评价与推广
根据传染病的模型建立研究进而推广产生了传染病动力学模型。传染病动力学[1]是对进行理论性定量研究的一种重要方法,是根据种群生长的特性,疾病的发生及在种群内的传播,发展规律,以及与之有关的社会等因素,建立能反映传染病动力学特性的数学模型,通过对模型动力学性态的定性,定量分析和数值模拟,来分析疾病的发展过程,揭示流行规律,预测变化趋势,分析疾病流行的原因和关键。
对于2003年发生的SARS疫情,国内外学者建立了大量的动力学模型研究其传播规律和趋势,研究各种隔离预防措施的强度对控制流行的作用,为决策部门提供参考.有关SARS传播动力学研究多数采用的是SIR或SEIR模型.评价措施效果或拟合实际流行数据时,往往通过改变接触率和感染效率两个参数的值来实现.石耀霖[2]建了SARS传播的系统动力学模型,以越南的数据为参考,进行了Monte Carlo实验,初步结果表明,感染率及其随时间的变化是影响SARS传播的最重要
因素.蔡全才[3]建立了可定量评价SARS干预措施效果的传播动力学模型,并对北京的数据进行了较好的拟合。
七、参考文献
[1]耀霖.SARS传染扩散的动力学随机模型[J].科学通报,2003,48(13)1373-1377
[2]唐焕文,贺明峰.数学建模引论.北京.高等教育出版社.