SARS的传播模型第九组

1、SARS的传播摘要 SARS（Severe Acute Respiratory Syndrome，严重急性呼吸道综合症，俗称非典型肺炎）是21世纪第一个在世界范围内传播的传染病。SARS的爆发和蔓延给我国的经济发展和人民生活带来了很大的影响。为了能定量的研究传染病的传播规律，人们建立了各类模型来预测、控制疾病的发生发展。 对于问题1，给出了一个早期指数模型，它在短期内有着计算参数简单等合理性与实用性，但却存在着用短期数据分析预测后期疫情发展趋势的缺陷。基于此，我们考虑引进新的参数，建立更优的模型。对于问题2我们在早期模型的基础上进行改进，建立了SIR模型，考虑到疑

2、似患者的变化情况，分别建立了模型一和模型二。对模型进行了合理的理论证明和推导，所给出的理论证明结果符合模型所给数据，然后运用相轨分析法以及借助Matlab，excel软件，对附件中所提供的数据进行了拟合和分析。最后根据模型的分析结果对卫生部门所采取的措施如提前或延后5天采取严格的隔离措施，是有数学依据的。 对于问题3我们根据表格所给数据，使用二次回归的方法先对03年1月份至12月份的数据进行预测，再使用Matlab中的作图工具箱做出未受SARS的影响图和受SARS的影响图，最后进行比较分析，得出结论：SARS对北京的旅游业造成影响。 关键字：SIR模型、Matlab、excel、二次

3、回归方程、相轨线一、问题重述SARS（非典型肺炎）的爆发和蔓延给我国的经济发展和人民生活带来了很大影响，因此定量地研究传染病的传播规律，对预测和控制传染病蔓延起着很大的作用。现对SARS 的传播建立数学模型，具体要求如下：（1）对附件1所提供的一个早期的模型，评价其合理性和实用性。（2）建立你们自己的模型，说明为什么优于附件1中的模型；特别要说明怎样才能建立一个真正能够预测以及能为预防和控制提供可靠、足够的信息的模型，这样做的困难在哪里？对于卫生部门所采取的措施做出评论，如：提前或延后5天采取严格的隔离措施，对疫情传播所造成的影响做出估计。附件2提供的数据供参考。（3）收集SARS对经济某个方

4、面影响的数据，建立相应的数学模型并进行预测。附件3提供的数据供参考。（4）给当地报刊写一篇通俗短文，说明建立传染病数学模型的重要性。二、问题分析2.1问题1的分析：对本文提供的一个早期的模型进行评价，主要是分析模型所采用的方法，对其合理性进行肯定，同时要指出其存在的不足。再与实际联系，分析该模型的实用性。2.2问题2的分析：根据对早期模型的分析，对其不足指出进行反思，建立起新的改进模型SIR，首先要对SARS的传播机理进行深刻的了解和分析；其次是根据传染病的特点建立起合适的模型SIR；最后对模型进行求解。2.3问题3的分析： 由1997年至2002年的旅游人数预测未受SARS影响时2003各月

5、份的旅游人数，并建立二次回归方程进行回归分析，接着由实际2003年1月至8月的数据预测9月到12月的旅游人数，建立二次回归方程拟合，比较SARS影响前后的图像，分析出其对旅游业的影响。三、 模型假设1、 假设总人口数保持不变，不考虑人口的流动因素2、 忽略当地的自然出生率和自然死亡率。3、 假设被治愈的病人有免疫力，不再被感染。4、 假设移出者包括治愈免疫者和死亡人群四、定义与符号说明：总人口数；：健康人占总人口的比例；：病人占总人口的比例；：移出者占总人口的比例；：日接触率，即每个病人每天接触的健康的概率；：日治愈率，即每天被治愈的病人的概率：日死亡率， 即每天的死亡人数占病人总数的比例：传

6、染期接触数，即有；：疑似感染率。即每天感染为疑似病人的比例；：疑似病人占总人口的比例；：日转化率，即每天危险群体中的疑似病人被确诊为SARS患者的比例。五、模型建立与求解5.1问题1的求解1、早期模型的合理性评价：（1）该模型简单易行，方便对数据的拟合，并容易分析出所建模型与实际数据的误差，能够具有一定的合理性。（2）对广东、香港、以及北京的疫情发展趋势的比对可以看出因地区、人口等因素的影响疫情发展趋势存在很大的差异。（3）该模型选用公布数据直接拟合，从而预测后期疫情发展趋势，用短期数据来分析，这样建模具有一定的局限性，缺乏合理性。2、早期模型的实用性评价：（1）该模型反应出一般传染病模型的发

7、展趋势“快速蔓延期、相对稳定期、逐渐消亡期”，具有一定的实用性，而该模型对于SARS传播发展的初期的研究有参考价值。（2）模型的参数K的的选择没有给出客观的算法或依据，人工的调整数据具有一定的主观性。而平均传染期限L固定在20，显得片面缺乏可靠性，因为平均传染期限是会随着疫情的发展而变化的。（3）该模型只是考虑了健康者和感染者，并没有考虑到治愈者能否具有免疫力的情况，实用性不强。5.2问题2的求解 SARS的传播机理：1、总人数不变时，将社会人群分为三类，称为SIR模型。S类：称为健康人，该类成员没有染上传染病，但缺乏免疫能力，可以被染上传染病.I类：称为病人，该类成员已经染上传染病，而且可以

8、传染给S类成员.R类：称为移出者，R类成员或者是I类成员被严格隔离、治愈，或者死亡等.I类成员转化为R类后，立刻失去传染能力.、分别表示t时刻上述3类成员占城市人口总数的比例.2、SARS的传播过程：5.2.2模型的建立 模型一 感染为SARS患者情况由假设可知，每个病人每天可使个健康者变为病人，因为病人人数为，所以每天共有个健康者被感染，于是就是病人数的增加率，又因为每天被治愈率为，死亡率为，所以每天有个病人被治愈，有个病人死亡。那么病人的感染为 （1）显然有： （2） 对于病愈免疫的移出者而言应有： （3）由（1）（2）（3）可得SIR模型如下： （4）模型二 疑似患者的变化情况类似前面的

9、分析，得到疑似患者率模型： （5）5.2.3模型的求解1、参数的确定：对附表2中的数据有excel处理：表格见附件表一当天的病人总数=隔天的确诊病例-当天确诊病例 =， =， =0.055076 =0.002443 =0.038183 （处理数据见附件1）故可得+=0.0575192、的确定确定 从我们建立的模型是无法得到、的解析解。故求出他们的数值解。先通过实际统计数据算出每一天的、做出它们与时间的函数图象图1， <图1>:根据实际数据拟合的图象（画图程序见附件2） 当天病人变化<图1>然后我们再对取一组数，分别画出由通过模型解出的数值解随时间变化的图象图2，将这组图

10、象与由实际数据所得图象相比较,调试。我们发现当1.5时，理论图形与实际图形有最佳的吻合。图形如下：通过数值解作出的关于时间t 的变化（画图程序见附件3）<图2>分析两个图形可知，它们的高峰期、缓解期和平稳期曲线相当符合，具有相同的发展趋势。但是在0,10的SARS初期范围内，曲线变化不相同。这主要是因为在4月24日之前，没有相关数据的统计和报道，由于数据的不全，根据边界值画出来的曲线与通过数值解得到的曲线相比较，不能准确反映SARS产生初期时的趋势，所以边界值应该去掉，而通过数值解模拟的曲线可以得到之前的发展趋势。并且通过对SARS蔓延期特点的分析，<图2>在符合所给数

11、据反映的规律基础上，还能够模拟缺乏数据的SARS初始状态，所以曲线是合理的。（2）确定与确定时类似，先根据实际数据画出图形（画图程序见附件4）<图3>实际数据图形然后再对取一组数，分别画出通过模型解出的数值解随时间变化的图象，将这组图象与由实际数据所得图象相比较,调试。发现当1.0时，理论图形与实际图形有最佳的吻合。图形如下（画图程序见附件5）： <图4> 疑似病人变化在0,10的初期范围内，曲线趋势不同，原因同前。整个曲线反映了疑似患者在SARS的过程中的变化规律。5.24结果的分析与验证（一）讨论 的性质平面称为相平面，相轨线在相平面上的定义域为从模型（一）中消去，

12、利用的定义，可得 （6）由（6）式解得 （7）（二）对于合理确定的，我们可以画出图，图形如下：（画图程序见附件6） <图5> 图形（相轨线）由于在这个SARS病毒发展过程中，是变化的，故可以画出取不同值时的图形，如下取0.4192，0.2858、0.1858时的图形，（画图程序见附件7）<图6>图形（相轨线）分析（3）式和（7）式，可知：1 不论初始条件，如何，病人终会消失，即SARS最终会被消灭,亦即。从图形上看，相轨线终将与s轴相交（t充分大）。2 SARS疾病传染过程分析整个传染过程，随着政府和公众对SARS的重视程度的变化，可知接触数随着治愈率、死亡率和接触率的

13、不断变化而变化。（1）在SARS爆发的初期，由于潜伏期的存在，社会对SARS病毒传播的速度和危害程度认识不够，所以政府和公众没有引起重视。治愈率和死亡率很小，而接触率相对较大，所以很小。当，则开始增加，可认为是疾病蔓延阶段。（2）当=时，达到最大值 （9）对于我们确定的，可以求出0.8368，可认为是疾病传染到达了高峰期。（3）当<时，单调减小至零，单调减小至。这一时期病人比例绝不会增加，传染病不会蔓延，进入缓解期。3群体免疫和预防根据对模型的分析，当是传染病不会蔓延。所以为制止蔓延，除了提高卫生和医疗水平，使阈值1/变大以外，另一个途径是降低，这可以通过预防接种使群体免疫。第二个途径通

14、过预防接种使群众免疫，免疫后就不会被感染上病毒。按照我们人群的分类系统，将免疫人群归为退出者类，所以免疫人群的出现，不与模型的分类系统相矛盾。忽略病人比例的初始值，有=1-，于是SARS不再蔓延的条件可以表示为： （10）所以只要通过群体免疫使初始时刻的移出者比例满足（10），就可以制止SARS的蔓延。4数值验证与估量根据上面的分析，阻止SARS蔓延有两种手段，一是提高卫生水平和医疗水平，即降低日接触率，提高日治愈率，二是群体免疫，即提高移出者比例的初值。我们以最终未感染的健康者的比例和病人比例达到最大值，作为传染病蔓延程度的度量指标。给定不同的，用（）式计算，用（9）式计算1.00.30.3

15、0.980.020.03980.34490.60.30.50.980.020.19650.16350.50.51.00.980.020.81220.02000.40.51.250.980.020.91720.02001.00.30.30.700.020.08400.16850.60.30.50.700.020.30560.05180.50.51.00.700.020.65280.02000.40.51.250.700.020.67550.0200从计算得到的和可以看出：（1）对于一定的，降低，提高，使阈值1/变大，会使变大，变小。于是验证了群体免疫和预防中提出的提高卫生水平和医疗水平，可以使S

16、ARS最终的患者比例缩小，健康群体增加。（2）对于一定的，提高 ，会使变大，变小。所以实行群体免疫，降低受感染的基数，可以有效地减缓SARS蔓延的速度。在（8）式中略去很小的，即有 （11）对于表达式中的参数，已通过前面的参数分析得出，代入表达式，就可以对t时的患病率做预测，达到了预测的目的，满足题目的要求。5.25 建立可靠的预测模型存在的困难要建立一个真正能够预测以及能为预防和控制提供可靠、足够的信息的模型，还存在着一些困难：1、对于我们所建立的SIR模型，是基于附表二给出的数据，但是所给的数据不够齐全，只是给出了一个时间段的病情数据，用来分析整体不够全面，而且对于所给的数据缺乏可靠性。2

17、、对于我们所建立的模型，参数的选择是很关键的，由于数据不全，参数的选择也存在着一定的误差。5.26对卫生部措施的评估在模型中，的取值大小能充分反映接触率的变化。若采取的隔离措施提前T天，那么将相应减小，反之则增加。不妨将的值取为1.3和1.8，作出相应的图形7和图8（画图程序见附件8）。图7 疑似病人变化图8疑似病人变化由以上图形可见，T对SARS病人的增长有显著影响，因此，卫生部采取的提前或延后5天的隔离措施有其数学背景和科学依据。5.3问题3的求解5.31回归模型的建立与求解：1、根据附表3中1997年到2002年的旅游人数预测为受SARS影响时2003年各月份的旅游人数，建立回归模型得出

18、2003年各月份的预测值为1月2月3月4月5月 6月7月8月9月10月11月 12月15.236.4725.9732.1532.8331.6 29.3336.433.1432.8526.8527.79根据预测值的散点图建立二次回归方程画出回归图如下（程序见附件9） 由图可看出：未受SARS影响的情况下，2003年的旅游人数变化与往年的旅游趋势相同，先升后降，在夏季达到旅游的高峰。 2、根据实际2003年1月到8月的数据预测9月到12月的旅游人数1月2月3月4月5月 6月7月8月9月10月11月 12月15.4 17.123.511.61.782.618.816.215.4120.1926.17

19、33.37同样建立二次回归方程，得到如下的回归图（程序见附件10）; 由图可看出：在SARS的影响下，2003年的旅游趋势先降后升，低谷恰好是在SARS的高峰期。3、对比SARS影响前后的回归图像： 5.32结果分析： 对比SARS影响前后的回归图像可以发现：在一月份左右，SARS对旅游人数的影响不大，数据基本吻合，这是因为疫情刚发生时没有引起人们的注意，随着疫情的蔓延，社会处于敏感期，旅游人数逐渐减少，在疫情爆发的高峰期，由图像可以看到，在六月份左右两条曲线的差距达到最大，社会处于高度的警界期，出门旅游的人很少，而随着疫情不断的得到控制，旅游业也渐渐回暖，图像的差距也渐渐缩小。因此，SARS

20、直接影响了旅游人数,对旅游行业的影响是很大的。六、模型的评价与推广1、本模型适用预测SARS及其他类似的传染病（比如禽流感病毒）高峰期的来临。能够为预防和控制传染病提供可靠足够的信息。2、本模型有利于在当今医学领域中，分析各种传染病或森林、农业、科学上病虫害的变化规律，度量传染病蔓延的程序并探索制止蔓延手段。七、附件附件1日 期已确诊病例累计现有疑似病例死亡累计治愈出院累计当天退出者当天病例退出率治愈率死亡率转化率当天病人变化数4月20日3394021833174310.0394430.0232020.0162414月21日482610254365200.0115380.0057690.005

21、7690.2344261434月22日5886662846166190.0258480.014540.0113090.1591591064月23日6937823555136840.0190060.0131580.0058480.1342711054月24日7748633964127740.0155040.0116280.0038760.093859814月25日877954427398730.0103090.0034360.0068730.1079661034月26日98810934876109900.0101010.002020.0080810.1015551114月27日111412555

22、678310650.00281700.0028170.1003981264月28日1199127559781212100.0099170.0041320.0057850.066667854月29日1347135866831612910.0123930.0054220.0069710.1089841484月30日1440140875901713880.0122480.0072050.0050430.066051935月1日15531415821001814540.012380.006190.006190.0798591135月2日16361468911091115410.0071380.0038

23、940.0032450.05654835月3日1741149396115715920.0043970.0018840.0025130.0703281055月4日18031537100118616790.0035740.0017870.0017870.040338625月5日189715101031211717360.0097930.0074880.0023040.062252945月6日196015231071341018080.0055310.0038720.0016590.041366635月7日204915141101411318850.0068970.0058360.0010610.0

24、58785895月8日213614861121521819130.0094090.0083640.0010450.058546875月9日21771425114168919450.0046270.0035990.0010280.028772415月10日222713971161751519740.0075990.0055720.0020260.035791505月11日226514111201863119980.0155160.0110110.0045050.026931385月12日230413781292084120100.0203980.017910.0024880.028302395月

25、13日234713381342441319920.0065260.0040160.002510.032138435月14日23701308139252619970.0030050.0025040.0005010.017584235月15日238813171402571720080.0084660.0079680.0004980.013667185月16日240512651412733820060.0189430.0169490.0019940.013439175月17日242012501453072719820.0136230.0126140.0010090.012155月18日2434125

26、01473322019580.0102150.0086820.0015320.0112145月19日243712491503495019450.0257070.023650.0020570.00240235月20日244412251543955418950.0284960.0274410.0010550.00571475月21日244412211564478318530.0447920.0437130.001079005月22日245612051585285617790.0314780.0303540.0011240.009959125月23日246511791605828817480.050

27、3430.0486270.0017160.00763495月24日249011341636674116690.0245660.0221690.0023970.022046255月25日249911051677044416330.0269440.0263320.0006120.00814595月26日250410691687478515970.0532250.050720.0025050.00467755月27日251210051728284115140.0270810.0250990.0019820.0079685月28日25149411758666314760.0426830.0420050

28、.0006780.00212525月29日25178031769287914160.0557910.0550850.0007060.00373635月30日252076017710068513380.0635280.0605380.002990.00394735月31日252174718110873712540.0295060.02950600.00133916月1日252273918111243312170.0271160.0271160006月2日252273418111573211840.0270270.0270270006月3日252272418111897411520.0642360

29、.0642360006月4日252271818112635810780.0538030.0538030006月5日252271618113218410200.0823530.0803920.001961006月6日25227131831403439370.0458910.0458910006月7日25236681831446988930.1097420.1086230.00112006月8日252255018415431107950.1383650.1383650006月9日25224511841653966850.1401460.1372260.002920.00221716月10日2522

30、3511861747745900.1254240.1254240006月11日25232571861821565160.1085270.1065890.001938006月12日25231551871876684590.1481480.1481480006月13日2522711871944523910.1329920.1278770.005115006月14日252241891994213390.0619470.0619470006月15日252231892015393170.1230280.1198740.003155006月16日252131902053672780.2410070.241

31、0070006月17日252151902120352110.1658770.1611370.004739006月18日252141912154171760.0965910.0965910006月19日252131912171181590.1132080.1132080006月20日252131912189421410.2978720.2978720006月21日25212191223126990.2626260.26262600.516月22日25212191225720730.2739730.2739730006月23日25212191227700总数3.6811993.5248440.15

32、63552.443075平均值0.0575190.0550760.0024430.038173附件2当天病人变化（确定1）t=1:64;z=143 106 105 81 103 111 126 85 148 93 113 83 105 62 94 63 89 87 41 50 38 39 43 23 18 17 15 14 3 7 0 12 9 25 9 5 8 2 3 3 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0;p,S=polyfit(t,z,9)Y=polyconf(p,t,S) plot(t,z,'o',t,Y),g

33、rid,附件3it：clccleara=1.5;b=0.0575;fun2=(t,x)a*x(1)*x(2)-b*x(1);-a*x(1)*x(2);ts=0:1:70;x0=402/14000000,1-402/14000000;t,x=ode45(fun2,ts,x0); t,x 附件4确定2t=1:65;z=4026106667828639541093125512751358140814151468149315371510152315141486142513971411137813381308131712651250125012491225122112051179113411051069

34、1005941803760747739734724718716713668550451351257155714335433222;p,S=polyfit(t,z,9)Y=polyconf(p,t,S) plot(t,z,'b-',t,Y,'g-'),grid,附件5it：clccleara=1;b=0.038183;fun2=(t,x)a*x(1)*x(2)-b*x(1);-a*x(1)*x(2);ts=0:1:70;x0=402/14000000,1-402/14000000;t,x=ode45(fun2,ts,x0); t,x plot(t,x(:,1),grid,附件6相轨线：clcclearfun2=(t,x)1.5*x(1)*x(2)-0.0575*x(1);-1.5*x(1)*x(2);ts=0:1:70;x0=402/14000000,1-402/14000000;t,x=ode45(fun2,ts,x0); t,x plot(x(:,2),x(:,1),grid,附件7改变值s=0.000001:0.001:1;i=1-s+0.4192*log(s/0.9999