SARS传播模型

1、SARS传播模型摘 要：本文中我们对北京地区4月20日-6月8日的SARS疫情数据进行了分析处理，把北京地区SARS疫情分为两个时期：感染期（4月20日-5月16日）和恢复期（5月17日-6月8日）。由于医务人员人群是SARS病的高发人群，所以我们在本文的模型中把病人人群分为医务人员病人人群和非医务人员病人人群。通过分析文中附件1的数据，我们建立了两个时期SARS传播的微分方程模型，并得到了模型的解，感染期模型的解为：，恢复期模型的解为：。从模型解的曲线与实际数据的比较，我们发现该模型的解与实际数据是非常吻合的。关键词：SARS，三次样条插值，高次多项式拟合0. 引 言SARS（Severe

2、Acute Respiratory Syndrome，严重急性呼吸道综合症, 俗称：非典型肺炎）是一种新的、传染性很强的疾病，它在我国部分地区的暴发与蔓延，严重威胁了人民的健康与生命安全，影响我国的社会稳定与经济发展。我们从整个抗击SARS的斗争中得到了许多重要的经验和教训，同时也认识到定量地研究传染病的传播规律、预测传染病发展趋势的必要性和重要性。1. 问题分析由于SARS主要是通过近距离的飞沫传播，与病人有过密切接触的人群就很可能被感染，成为SARS病人，所以我们在模型里假设健康者只要与病人接触，则感染为病人，且由于医务人员每天都直接与病人接触，所以医务人员人群是SARS病的高发人群，其主

3、要是被住院的病人传染。而普通病人即非医务人员的病人，主要是由一些病人在发病后未及时被隔离治疗而与健康人接触，并使其感染病毒，因此我们在模型里把病人分为医务人员病人和普通病人。同时临床统计数据表明SARS病的潜伏期为2-12天，一般在4-5天，治愈后的病人没有出现再次患病现象，所以我们也假定治愈后的病人具有免疫力。虽然SARS病在2002年底到2003年初就在我国各地市广泛传播，但是疫情数据并没有精确的得到统计，从4月20日后，国家卫生部才在专门的网站上发布各地的SARS疫情数据，本文中我们只对北京地区4月20日-6月8日的数据进行分析。通过对本文附录1中数据的处理，我们发现北京市4月20日-6

4、月8日的疫情可分为两个时期：感染期（4月20日-5月16日）和恢复期（5月17日-6月8日）（见附录2）。并且通过拟合附件1中所给的数据，我们得到了两个时期内医务人员及普通人员的日感染率k，日死亡率d和日治愈率z，建立了SARS传染病的一个微分方程模型。2. 模型假设1) 本文主要考察的是北京地区的疫情变化, 在疾病传播期内该地区的人口总数N=10000000不变，即为常量。2) 政府发布的北京地区的疫情数据真实可信。3) SARS传播途径都视为与病源的直接接触, 每个病人与健康者接触时，都使健康者感染变成病人,且病人发病后马上被隔离，即住院治疗 。4) 将人群分为三类：健康者（易受感染者）、

5、病人（已被感染者）、退出者（死亡者和治愈者）,在t时刻三者在总人数中的比例分别为s(t)、i(t)、r(t)，且s(t)+i(t)+r(t)=1。其中病人又分为医务人员病人和非医务人员即普通病人。5) 假设不考虑SARS传播期内的北京地区的人口出生率和自然死亡率。对于由SARS引起的死亡人数全部归为“退出者”, 治愈者不会再次被感染，具有免疫力，归为“退出者”。6) 假设每日新增的普通病人与当日健康人的人数成正比，比例系数为k1,称为普通病人的日感染率。每日新增的医务人员病人与当日住院的病人人数成正比，比例系数为k2,称为医务人员的日感染率。每日新增的病人死亡人数与当日已确诊的病人数成正比，比

6、例系数为d,称为病人的日死亡率。每日新增治愈的病人人数与当日已确诊的病人数成正比，比例系数为z,称为病人的日治愈率。由于k1，k2,d,z具有很强的实际意义，其具体的数值我们可以通过对给定的数据拟合得到。3. 符号说明，,，i0,s0分别为病人数，健康人数占总人口数的比例初始值。4. 模型建立从模型假设中的6)可知普通病人的日感染率k1, 医务人员的日感染率k2, 病人的日死亡率d, 病人的日治愈率z, 我们可以通过对附录1中的实际数据处理来得到。为了使得到的模型比较简单，我们考虑k1,k2,d,z都是常数，这时k1、k2、d、z分别表示普通病人的日感染率、医务人员的日感染率、病人的日死亡率、

7、病人的日治愈率的平均值，对感染期与恢复期我们可以建立下面相同的模型，这两个时期的差别主要在于普通病人的日感染率k1是否等于0(见附录2)。由上面的分析及假设，我们得到了下面的SARS病传播的微分方程模型： (1)5. 模型求解51 对感染期的模型求解我们考虑4月20日至5月16日，由附录1中所给的数据，经过拟合得到k1=5.3527401482552×10-6，k2=0.01742890531138，d=0.00453527706374，z=0.02388800922302。且以4月20日的疫情数据作为初始值，即。对于方程组（1）我们很难得到它的解析解，但由方程组（1）的前两个方程我

8、们消去dt，可以把（1）转化为下面的方程组： （2）于是令，其中把s看作关于i的函数得，即 。由于，再加上初始条件s（0）、i（0），从而我们可以得到方程组（2）的解：,其中，。表示t时刻，健康人数与病人人数的比值，它反映了传染病的严重程度。当t时刻u(t)越小，即病人占总人数的比例i(t)越大，健康人数的比例s(t)越小，则SARS疫情严重。反之，t时刻u(t)越大，则SARS疫情得到有效控制。 从4月20日至5月16日的SARS数据可得到u(t)的变化情况(如图1.)。从图1.中我们看到4月20日至5月16日，的值随着t的变化在急剧减少，主要是由于s(t)减少的同时i(t)也在快速增加，而

9、5月17日以后u(t)又开始上升，这是由于从5月17日北京市的疫情得到了有效控制，新增的SARS病例被控制在1例以内，健康人数不再减少。相反，病人每天逐渐被治愈，i(t)也在减少。图1.通过给模型中的解代入具体的数据，我们可以得到4月20日至5月16日病人数所占的总人数的比例i(t)的变化情况(如图2.)。可以看出模型解与实际数据的值非常吻合，这里我们所用的u(t)是健康人数与病人数的比值的三次样条插值函数，u(t)与实际数据越吻合则我们的解的曲线也就更加吻合于实际数据。图2.52 对恢复期的模型求解考虑4月20日至5月16日，由附录1中所给的数据，经过拟合得到k1=0，k2=2.386572

10、17981×10-4，d=0.00117550039945，z=0.04137494891164。且以5月16日的疫情数据作为初始值，即。将其代入5.1的方程组（1）可得恢复期的模型： (3)通过消去方程组（3）中的dt，则方程组（3）转化为 （4）其解为：，其中。令，则方程组（4）的解可转化为：。我们仍然用5.1中得到的u(t)，即健康人数与病人数的比值的三次样条插值函数，则恢复期模型的解的图形如下：图3.53 北京市4月20日至6月8日SARS的疫情传播情况由5.1和5.2的结果，我们得到了北京市4月20日至6月8日SARS的疫情传播过程中，病人数占总人数的比例的变化趋势图（如图

11、4.）。图4.6. 模型推广由于每一天k1,k2,d,z的值是不同的，即它们是随着时间而变化的值，所以我们可以考虑k1,k2,d,z都是关于时间t的函数，即在t时刻医务人员的日感染率k2(t)、普通病人的日感染率k1(t)、病人的日死亡率d(t)、病人的日治愈率z(t)，再利用实际数据经过多项式拟合得到k1(t)、k2(t)、d(t)、z(t)的多项式函数，或者用三次样条函数进行插值。这样使得模型的解与实际数据更加的吻合。同时我们也可以注意到SARS病的有效控制，很大程度上是由于传染病人在发病后马上被隔离，得到治疗。但是现实生活中不可能做到这样，在病人发病到被隔离治疗之间的这段时间 L内，病人

12、就会与健康人接触而感染健康人。因此L的值在建立SARS数学模型中非常重要。所以本模型也可以在这些方面进行改进。在得到的模型的解后，我们为了看到模型的解与实际数据的接近程度，通过三次样条函数进行插值u(t)的值，然后代入解中来分析病人数占总人数的比例i(t)随时间的变化情况。对于u(t)，我们也可以用高次多项式或分段高次多项式来处理，从下列图中我们发现选用的多项式次数越高，模型的解的值与实际数据的接近程度就更好，但是多项式次数越高就会使数据处理量加大，并且得到的解不稳定，出现病态数据。图5.7. 参考文献1 姜启源编，数学模型，北京：高等教育出版社，2003年4月。2 Wliliam F.luc

13、as主编，微分方程模型,长沙:国防科技大学出版社,1998年5月。3 王高雄等编，常微分方程，北京：高等教育出版社，1999年3月。4 龚建华等著，SARS疫情控制的模拟分析，遥感学报，2003,7(4):260-265。5 Shoichiro Nakamura著，科学计算引论-基于MATLAB的数值分析，北京：电子工业出版社，2002年6月。6 傅鹂等编著，数学实验，北京：科学出版社，2001年3月。7 叶其孝主编，大学生数学建模竞赛辅导教材（一），长沙：湖南教育出版社，1998年4月。附录1：北京市疫情的数据（ 数据来源：中华人民共和国卫生部广东非典型肺炎防治网中国疾病预防控制中心 ）日

14、期已确诊病例累计医务人员现有疑似病例死亡累计治愈出院累计4月20日3392440218334月21日4827861025434月22日5889966628464月23日69312678235554月24日77414386339644月25日87716095442734月26日988167109348764月27日1114187125556784月28日1199206127559 784月29日1347255135866834月30日1440270140875905月1日15532881415821005月2日16363001468911095月3日17413191493961155月4日180

15、332615371001185月5日189733515101031215月6日196034815231071345月7日204935815141101415月8日213637014861121525月9日217737214251141685月10日222737513971161755月11日226537914111201865月12日230438413781292085月13日234738413381342445月14日237038713081392525月15日238838913171402575月16日240539112651412735月17日242039312501453075月18日

16、243439412501473325月19日243739612491503495月20日244439512251543955月21日244439512211564475月22日245639512051585285月23日246539511791605825月24日249039611341636675月25日249939611051677045月26日250439710691687475月27日251239810051728285月28日25143989411758665月29日25173988031769285月30日252039876017710065月31日2521399747181108

17、76月1日252239973918111246月2日252239973418111576月3日252239972418111896月4日252239971818112636月5日252239971618113216月6日252239971318314036月7日252339966818314466月8日25223995501841543附录2从附录1的数据中，我们知道确诊的医务人员的人数从4月20日开始一直在增长，而且4月20日至5月8日增长速度很快，从5月9日至6月8日增长速度趋于缓和（如图6.）。4月20日至5月8日之间每日新增的医务人员病人数变化很大，这也反映出北京市SARS疫情这时期十

18、分复杂，有很多原本处在潜伏期内的SARS感染者发作被确诊，因而被隔离住院治疗，同时很多医务人员在与SARS病患者接触时被感染。5月9日至6月8日每日新增的医务人员病人数降为几例，说明这个时期北京的SARS疫情已经得较好的控制（如图7.），被确诊的SARS病人及住院的病人明显降低，同时也反映出国家采取的干预政策起了很好的效果，人们的防范意识得到了增强。 图6. 图7.图8.中我们可以看出从4月20日至5月16日非医务人员的病人数一直在增长，在5月16日出现了普通人员感染人数的最大值，从5月17日至6月8日普通人员的病人数开始下降（如图8.）。同时我们通过图9.知道在5月16日以后新增的普通人员的

19、病人数为负值，这主要是由于普通人员没有新的病人增加，相反每天有一些病人被治愈而导致。因此我们在模型中把北京市4月20日至6月8日的疫情分成两个时期考虑，即4月20日至5月16日称为感染期，5月17日至6月8日称为恢复期。在恢复期内，我们假设这时期内普通人员病人的感染率为0。 图8. 图9.图10.中每日新增的治愈病人数在5月15日之前变化很小，5月15日以后快速上升，图11.中每日死亡的病人数5月15日之前变化很大，5月15日之后显著下降，基本控制在1个人以内。因此可以看出北京市的SARS疫情从5月15日之后得到了有效的控制。 图10. 图11.附录3（问题处理的Matlab源程序）% 下列程

20、序计算医务人员的日感染率k2， %普通病人的日感染率k1，死亡率d，治愈率z %clear allformat longN=10000000;%北京SARS确诊，疑似，治愈，死亡人数的数据（4月20日-6月8日） %qz=339,482,588,693,774,877,988,1114,1199,1347,1440,1553,1636,1741,1803,1897,1960,2049,2136,2177,2227,2265,2304,2347,2370,2388,2405,2420,2434,2437,2444,2444,2456,2465,2490,2499,2504,2512,2514,2

21、517,2520,2521,2522,2522,2522,2522,2522,2522,2523,2522; dead=18,25,28,35,39,42,48,56,59,66,75,82,91,96,100,103,107,110,112,114,116,120,129,134,139,140,141,145,147,150,154,156,158,160,163,167,168,172,175,176,177,181,181,181,181,181,181,183,183,184; cure=33,43,46,55,64,73,76,78,78,83,90,100,109,115,118

22、,121,134,141,152,168,175,186,208,244,252,257,273,307,332,349,395,447,528,582,667,704,747,828,866,928,1006,1087,1124,1157,1189,1263,1321,1403,1446,1543;yw=24,78,99,126,143,160,167,187,206,255,270,288,300,319,326,335,348,358,370,372,375,379,384,384,387,389,391,393,394,396,395,395,395,395,396,396,397,3

23、98,398,398,398,399,399,399,399,399,399,399,399,399;pt=qz-yw-dead-cure;good=N-qz;infect=qz-dead-cure;%计算普通病人的日感染率disp('普通病人4月20日至5月16日，5月16日至6月8日的感染率分别为k11,k12')n=length(pt);k11=0;k12=0;j=0;for i=1:n-1 r(i)=pt(i+1)-pt(i); if(r(i)>0) k11=k11+r(i)/good(i+1); j=j+1; end end k11=k11/j k12%已确诊的

24、普通病人数的数据图形描绘%i=1:1:length(pt);%set(gcf,'position',300 200 560 420);%plot(pt,'ko')%hold on%plot(i,pt,'k-')%xlabel('时间(单位：天)');%ylabel('普通病人数');%title('普通病人数随时间的变化情况');%新增的普通病人数的数据图形描绘%i=1:1:n-1;%set(gcf,'position',300 200 560 420);%plot(r,'k

25、o')%hold on%plot(i,r,'k-')%xlabel('时间(单位：天)');%ylabel('每日新增的普通病人数');%title('新增的普通病人数随时间的变化情况');%计算医务人员的感染率disp('医务人员病人4月20日至5月16日，5月16日至6月8日的感染率分别为k21,k22')n=length(yw);k21=0;k22=0;for i=1:n-1 r(i)=yw(i+1)-yw(i); if(i<=25) k21=k21+r(i)/infect(i+1); else

26、 k22=k22+r(i)/infect(i+1); endend k21=k21/25 k22=k22/(n-26)%已确诊的医务人员病人人数的数据图形描绘%i=1:1:n;%set(gcf,'position',300 200 560 420);%plot(yw,'ko')%hold on%plot(i,yw,'k-')%xlabel('时间(单位：天)');%ylabel('已确诊的医务人员病人的人数');%title('已确诊的医务人员病人的人数随时间的变化情况');%新增的医务人员病人数随

27、时间的变化情况图形描绘%i=1:1:n-1;%set(gcf,'position',300 200 560 420);%plot(r,'ko')%hold on%plot(i,r,'k-')%xlabel('时间(单位：天)');%ylabel('每日新增的医务人员病人数');%title('新增的医务人员病人数随时间的变化情况');%计算病人的日死亡率disp('病人4月20日至5月16日，5月16日至6月8日的死亡率分别为d1,d2')n=length(dead);d1=0;d2

28、=0;for i=1:n-1 r(i)=(dead(i+1)-dead(i); if(i<=25) d1=d1+r(i)/infect(i+1); else d2=d2+r(i)/infect(i+1); endend d1=d1/25 d2=d2/(n-26)%已确诊的死亡病人数随时间的变化情况图形描绘%i=1:1:n;%set(gcf,'position',300 200 560 420);%plot(dead,'k-')%hold on%plot(i,dead,'ko')%xlabel('时间(单位：天)');%yla

29、bel('死亡病人数');%title('死亡病人数随时间的变化情况');%新增死亡的病人数随时间的变化情况图形描绘%i=1:1:n-1;%set(gcf,'position',300 200 560 420);%plot(r,'k-')%hold on%plot(i,r,'ko')%xlabel('时间(单位：天)');%ylabel('每日新增的死亡病人数');%title('新增的死亡病人数随时间的变化情况');%计算病人的日治愈率disp('病人4月2

30、0日至5月16日，5月16日至6月8日的治愈率分别为z1,z2')n=length(cure);z1=0;z2=0;for i=1:n-1 r(i)=(cure(i+1)-cure(i); if(i<=25) z1=z1+r(i)/infect(i+1); else z2=z2+r(i)/infect(i+1); endend z1=z1/25 z2=z2/(n-26)%治愈病人数随时间的变化情况图形描绘%i=1:1:n;%set(gcf,'position',300 200 560 420);%plot(cure,'k-')%hold on%pl

31、ot(i,cure,'ko')%xlabel('时间(单位：天)');%ylabel('死亡的病人数');%title('病人死亡人数随时间的变化情况');%新增治愈的病人数随时间的变化情况图形描绘%i=1:1:n-1;%set(gcf,'position',300 200 560 420);%plot(r,'k-')%hold on%plot(i,r,'ko')%xlabel('时间(单位：天)');%ylabel('每日新增的治愈病人数');%ti

32、tle('新增的治愈病人数随时间的变化情况');% % 下面的程序求模型的解i(t)的数值解 % %计算健康人数与病人数的比值n=length(qz);for i=1:n infect1(i)=infect(i)/N;%计算每日病人数占总人数的比值 t(i)=good(i)/infect(i);%计算健康人数与病人数的比值 end%对健康人数与病人数的比值做多项式拟合i=1:1:n;u=interp1(i,t,'spline'); %三次样条插值%f=polyfit(i,t,3) %高次多项式拟合%u=polyval(f,i);%4月20日-6月8日模型解的统一表示for i=1:n if (i&l