SARS传播的数学模型论文

摘要

本文分析了题目所提供的早期SARS传播模型的合理性与实用性，认为该模型可以预测疫情发展的大致趋势，但是存在一定的不足.第一，混淆了累计患病人数与累计确诊人数的概念；第二，借助其他地区数据进行预测，后期预测结果不够准确；第三，模型的参数L、K的设定缺乏依据，具有一定的主观性.

针对早期模型的不足，在系统分析了SARS的传播机理后，把SARS的传播过程划分为：征兆期，爆发期，高峰期和衰退期4个阶段.将每个阶段影响SARS传播的因素参数化，在传染病SIR模型的基础上，改进得到SARS传播模型.采用离散化的方法对本模型求数值解得到：北京SARS疫情的预测持续时间为106天，预测SARS患者累计2514人，与实际情况比较吻合.

应用SARS传播模型，对隔离时间及隔离措施强度的效果进行分析，得出结论：“早发现，早隔离”能有效减少累计患病人数；“严格隔离”能有效缩短疫情持续时间.

在建立模型的过程中发现，需要认清SARS传播机理，获得真实有效的数据.而题目所提供的累计确诊人数并不等于同期累计患病人数，这给模型的建立带来不小的困难.

本文分析了海外来京旅游人数受SARS的影响，建立时间序列半参数回归模型进行了预测，估算出SARS会对北京入境旅游业造成23.22亿元人民币损失，并预计北京海外旅游人数在10月以前能恢复正常.

最后给当地报刊写了一篇短文，介绍了建立传染病数学模型的重要性.

4． SARS传播的SIR改进模型

4.1模型的假设

1．SARS的持续期不太长，可以忽略在SARS持续期内的城市人口的自然出生率和自然死亡率.

2．被SARS感染后经治疗康复的人群在SARS流行期不会被再次感染.

3．病人被严格隔离、治愈或者死亡后，不再有感染作用.

4．不考虑人口的流动，仅仅在一个城市范围内研究SARS疫情的发展过程.

4.2模型的符号定义

:易感类人群占城市人口总数的比例.

:传染类人群占城市人口总数的比例.

:排除类人群占城市人口总数的比例.

:SARS患者的就诊率

：单位时间内一个传染者与他人的接触率.

：平均传染期.

4.3传播机理分析 

针对早期模型的不足，需要在模型的合理性和实用性方面进行改进.考虑在经典传染病模型SIR的基础上，通过机理分析，用实际因素来描述SARS的传播过程.

为了简化模型，这里不考虑人口的流动带来的影响，仅仅在一个封闭城市中研究SARS的传播机理.那么，整个社会人群可以分为3类：

S类：称为易感类，该类成员没有染上传染病，但缺乏免疫能力，可以被染上传染病.

I类：称为传染类，该类成员已经染上传染病，而且可以传染给S类成员.

R类：称为排除类或恢复类，R类成员或者是I类成员被严格隔离、治愈，或者死亡等.I类成员转化为R类后，立刻失去传染能力.

S(t)、I(t)、R(t)分别表示t时刻上述3类成员占城市人口总数的比例.

对于传播过程有3条基本假设：

：人口总数为常数N，N足够大，可以把变量S(t)、I(t)、R(t)视为连续变量，还可进一步假定为连续可微变量.

：人群中3类成员均匀分布，传播方式为接触性传播.单位时间内一个传染者与他人的接触率为，则一个传播者在单位时间内与S类成员的接触率为，因此，单位时间内I类成员与S类成员的接触总数为，这就是单位时间内I类成员增加的数量，称为发病率，它是S(t)和I(t)的双线性函数.

：传播者的被控制数正比于传染者的数量，比例系数为，称为被控制率，则平均传染期为.为一个传染者在其传播期内与其他成员的接触总数，称为接触数.

那么SARS的传播流程如图3：

图3   SARS传播流程图

在这个模型中，排除类就是已确诊SARS患者累计数，而是全社会累计SARS患者数，包括已确诊的和未被发现的两部分.

4.4模型的建立

有了以上的机理分析，建立起针对SARS的改进SIR模型：

该模型中参数和在疫情发展的各个阶段受实际因素影响，会有比较明显的变化，现分析如下：

参数表示单位时间内一个传染者与他人的接触率，其与全社会的警觉程度和、公众采取的各种措施有关，例如，佩戴口罩，减少停留在公共场所的时间，喷洒消毒药剂，提高隔离强度等都能有效地降低接触率的值.

一般认为，的数值随着SARS发展的4个阶段不断变化.在SARS初期，由于潜伏期的存在和社会对SARS病毒传播的速度认识不足，和公众并未引起重视，故维持在一个较高的数值；进入爆发期后，公众发现感染者不断增加，恐慌情绪增加，随即采取多种措施，使得到一定的控制，但效果不明显，此处假设呈线性形式缓慢衰减；在高峰期，当高强度的控制措施实施后，病毒传播的有效接触率明显减少，可以认为按天数呈指数形式衰减；此后进入衰减期，就维持在一个较低值附近.

参数表示传播者的被控制率.称为平均传染期，表示一个传播者在被隔离或者死亡之前具有传播能力的平均时间.一般认为，SARS患者经过传染期L过后，将隔离治疗或者死亡，从I类成员变为R类，失去传播能力.

L与采取的措施密切相关，例如，尽量早地发现病患，对疑似病例提前进行隔离，“早发现，早隔离” ；提供更广范围的医疗手段，使更多的人接受有效的治疗等，都可以有效地降低平均传染期L的长度.因此这里将L直接抽象为每一时期SARS患者的就诊率的函数.

平均传染期L应随的变化而变化.但是在初期，由于对SARS的认识不足，并没有采取有效控制措施， L的变化很小可以近似看作定值，这里我们取SARS病毒最长潜伏期（约19天）为这个定值；在爆发期，有效控制措施的逐步加强，使SARS患者的就诊率逐渐增加，而平均传染期L会逐渐减小并趋于一个定值，这里我们将SARS病毒平均潜伏期（约7天）定为L的最小值；在此后的高峰期以及衰减期，由于控制措施都保持在一定水平，L的值会维持在7天左右.

4.5针对北京疫情求解模型

首先采用数学推导的方法，确定参数和，并证明模型有唯一解.

确定和的关系

令，方程组中得：

在病情刚开始时，，由于是单调减少的，且最终趋近于0，则当时，单调减少趋近于0；当时，先单调增加达到最大值，然后单调减少趋近于0.

容易知道，当时，才满足SARS的传播规律，所以参数和的取值必须满足这个条件.

证明模型有唯一解

在初值条件下解微分方程组：

得到关系式：

令，由得

因为，所以令

则                 ， 

当时，由于在范围内有根，因而在内有根.

当时，因为

当时，，所以，因而在内也有根.

注意到当时，，故在内有唯一根.

所以，在内有唯一解.

划分SARS传播的4个阶段

由于SARS的传播经历了4个阶段，所以，要以具体的指标划分这4个阶段.因为在4个阶段中，日发病率是一个区分每个阶段特点的关键特征，所以以日发病率作为划分的指标.从第一个患者出现日开始：

征兆期：日发病率在10（人/天）以下.北京疫情期的前40天.

爆发期：从日发病率10（人/天）到日发病率最大，即时.北京疫情期的第40天到第74天.

高峰期：从日发病率最大到患者数量最大，即时.北京疫情期的第74天到第79天.

衰退期：患者数量最大点以后.北京疫情期第79天以后.

确定和

根据北京最终SARS患者总数2521人以及北京人口总数（约14000000人），得，所以.

因为平均传染期，而是SARS患者就诊率的函数，且，所以，这里设计L函数为：

由的控制措施决定，它的变化反映了控制措施的力度.根据实际情况，推导出：

而接触率与全社会的警觉程度和公众采取的各种措施有关，根据实际情况确定为：

确定出所有的参数后，做出北京各时期累计全社会SARS患者数和各时期累计确诊SARS患者数预测图（图4）以及北京市预测确诊SARS患者累计和实际确诊SARS患者累计对比图（图5）.同时得到：北京SARS疫情的预测持续时间为106天，预测SARS患者累计2514人.(计算程序见附件1：SIR模型程序)

图4   北京市预测非典病人累计总数和预测非典病人确诊病例累计对比图

图5  北京市预测确诊病例累计和实际确诊病例累计对比图

附件

附件1：SIR模型程序

function f=sor

S(1)=14000000;I(1)=1;R(1)=0;na=0.126;F=19;L=19;JU=19;M(1)=1;

for i=2:74                                             %初期与爆发期

   if i>=40&i<74

      JU=JU-0.25;

   end

   if i>=40

      na=na-0.01/35;                                   %爆发期缓慢减少

   end

   S(i)=S(i-1)-na*S(i-1)*I(i-1)/14000000;              %求解S，I，R

   if i>L+2

      R(i)=S(i-L-2)-S(i-L-1);

   else

      R(i)=0;

   end

   if i>=51

      F=F-0.5;L=fix(F);

      if F==L

         R(i)=S(i-L-3)-S(i-L-1);

      end

   end

   I(i)=I(i-1)+na*S(i-1)*I(i-1)/14000000-R(i);

   t=log(abs(14000001-S(i)))/log(10);

   o(i)=abs(14000001-S(i));p=log(o(i)-o(i-1))/log(10);

   plot(i+JU-19,t,'sr'),hold on,plot(i+JU-19,t,'sr'),

   if p>0.1

      plot(i+JU-19,p,'or'),plot(i+JU-19,p,'or')

   end

end

h=na;g(1)=17;n=1;

for i=75:139                              %高峰期与衰减期

   n=n+1;

   if i<80

      na=h-log(i-73)/35;                  %在高峰期 急剧减少，此后为一定值

   end

   R(i)=S(i-L-2)-S(i-L-1);                %求解S，I，R

   S(i)=S(i-1)-na*S(i-1)*I(i-1)/14000000;

   I(i)=I(i-1)+na*S(i-1)*I(i-1)/14000000-R(i);

   t=log(abs(14000001-S(i)))/log(10);o(i)=abs(14000001-S(i));

   p=log(o(i)-o(i-1)+1)/log(10);

   plot(i+JU-19,t,'sr'),plot(i+JU-19,p,'or')

end

附件2：时间序列程序

function sjxl

X(1,:)=[9.4 11.3 16.8 19.8 20.3 18.8 20.9 24.9 24.7 24.3 19.4 18.6];

X(2,:)=[9.6 11.7 15.8 19.9 19.5 17.8 17.8 23.3 21.4 24.5 20.1 15.9];

X(3,:)=[10.1 12.9 17.7 21.0 21.0 20.4 21.9 25.8 29.3 29.8 23.6 16.5];

X(4,:)=[11.4 26.0 19.6 25.9 27.6 24.3 23.0 27.8 27.3 28.5 32.8 18.5];

X(5,:)=[11.5 26.4 20.4 26.1 28.9 28.0 25.2 30.8 28.7 28.1 22.2 20.7];

X(6,:)=[13.7 29.7 23.1 28.9 29.0 27.4 26.0 32.2 31.4 32.6 29.2 22.9];

for i=1:12

   y(i)=0;

   for j=1:6

      y(i)=y(i)+X(j,i);

   end

   y(i)=y(i)/6;

end

for i=1:6                           %计算、和

   for j=1:12

      Z(i,j)=X(i,j)-y(j);

   end

end

for i=1:6

   p(i)=i-3.5; 

end

b=0;q=0;

for i=1:6

   for j=1:12

      b=b+p(i)*Z(i,j);

   end

   q=q+p(i)^2;

end

b=b/q/12; 

for i=1:12

   g(i)=y(i)-b*7/2;

end

si=0;

for i=1:6

   for j=1:12

      si=(Z(i,j)-b*p(i))^2;

   end

end

si=sqrt(si/59);

for i=1:12                          %预测2003年无SARS情况下每个月的旅游人数

   X7(i)=7*b+g(i)+rand*si;

end