SARS传播的数学模型论文

基于微分方程对象建模及实现（二）摘要：SARS时疫对中国社会发展产生了重大影响，本论文以传统的微分方程为理论基础，以2003年6月以前的有关SARS的数据为参考资料，着重从数学的角度研究和预测其发展趋势。本文分析了题目所提供的早期SARS传播模型的合理性与实用性，认为该模型可以预测疫情发展的大致趋势，但是存在一定的不足.第一，混淆了累计患病人数与累计确诊人数的概念；第二，借助其他地区数据进行预测，后期预测结果不够准确；第三，模型的参数L、K的设定缺乏依据，具有一定的主观性.针对早期模型的不足，在系统分析了SARS的传播机理后，把SARS的传播过程划分为：征兆期，爆发期，高峰期和衰退期4个阶段.将每个阶段影响SARS传播的因素参数化，在传染病SIR模型的基础上，改进得到SARS传播模型.采用离散化的方法对本模型求数值解得到：北京SARS疫情的预测持续时间为106天，预测SARS患者累计2514人，与实际情况比较吻合.关键字：微分方程概率平均SARS传播的SIR改进模型1.模型的假设1.假设SARS的传播方式为接触性传播，不与患病者接触就不会被传染2.假设人们被感染后需先进入潜伏期，在潜伏期内不具备传染性3.假设SARS患者被发现后就立即被隔离，被隔离者不具备传染性，SARS患者只在被发现前可以传染他人4.假设SARS康复者不会被再次感染，并且不具备传染性5.假设预情时间短不考虑在SARS传播期间人口的自然出生和自然死亡6.所研究地区的人口总量一定，不考虑该段时间内人口的迁入迁出2．模型的符号说明)(tS:易感类人群占城市人口总数的比例.)(tI:传染类人群占城市人口总数的比例.)(tR:排除类人群占城市人口总数的比例.)(t:SARS患者的就诊率患者总数时刻全社会患者数时刻被隔离的SARSSARStt：单位时间内一个传染者与他人的接触率.L：平均传染期.3．传播机理分析针对早期模型的不足，需要在模型的合理性和实用性方面进行改进.考虑在经典传染病模型SIR的基础上，通过机理分析，用实际因素来描述SARS的传播过程.为了简化模型，这里不考虑人口的流动带来的影响，仅仅在一个封闭城市中1研究SARS的传播机理.那么，整个社会人群可以分为3类：S类：称为易感类，该类成员没有染上传染病，但缺乏免疫能力，可以被染上传染病.I类：称为传染类，该类成员已经染上传染病，而且可以传染给S类成员.R类：称为排除类或恢复类，R类成员或者是I类成员被严格隔离、治愈，或者死亡等.I类成员转化为R类后，立刻失去传染能力.S(t)、I(t)、R(t)分别表示t时刻上述3类成员占城市人口总数的比例.对于传播过程有3条基本假设：1A：人口总数为常数N，N足够大，可以把变量S(t)、I(t)、R(t)视为连续变量，还可进一步假定为连续可微变量.2A：人群中3类成员均匀分布，传播方式为接触性传播.单位时间内一个传染者与他人的接触率为，则一个传播者在单位时间内与S类成员的接触率为)(tS，因此，单位时间内I类成员与S类成员的接触总数为)()(tItSN，这就是单位时间内I类成员增加的数量，称为发病率，它是S(t)和I(t)的双线性函数.3A：传播者的被控制数正比于传染者的数量)(tNI，比例系数为v，v称为被控制率，则平均传染期为vL/1.v/为一个传染者在其传播期内与其他成员的接触总数，称为接触数.那么SARS的传播流程图：)()()(tNRtNItNSvNSNSI排除类传染类易感类控制传染在这个模型中，排除类)(tNR就是已确诊SARS患者累计数，而)](1[tSN是全社会累计SARS患者数，包括已确诊的和未被发现的两部分.4模型的建立有了以上的机理分析，建立起针对SARS的改进SIR模型：00,01(2)(1)000RISSRIvIdtdRvISIdtdISIdtdS该模型中参数和v在疫情发展的各个阶段受实际因素影响，会有比较明显的变化，现分析如下：○1参数表示单位时间内一个传染者与他人的接触率，其与全社会的警觉程度和政府、公众采取的各种措施有关，从而有效地降低接触率的值.2一般认为，的数值随着SARS发展的4个阶段不断变化.在SARS初期，由于潜伏期的存在和社会对SARS病毒传播的速度认识不足，政府和公众并未引起重视，故维持在一个较高的数值；进入爆发期后，公众发现感染者不断增加，恐慌情绪增加，随即采取多种措施，使得到一定的控制，但效果不明显，此处假设呈线性形式缓慢衰减；在高峰期，当高强度的控制措施实施后，病毒传播的有效接触率明显减少，可以认为按天数呈指数形式衰减；此后进入衰减期，就维持在一个较低值附近.○2参数v表示传播者的被控制率.vL/1称为平均传染期，表示一个传播者在被隔离或者死亡之前具有传播能力的平均时间.一般认为，SARS患者经过传染期L过后，将隔离治疗或者死亡，从I类成员变为R类，失去传播能力.L与政府采取的措施密切相关，例如，对疑似病例提前进行隔离，“早发现，早隔离”；提供更广范围的医疗手段，使更多的人接受有效的治疗等，都可以有效地降低平均传染期L的长度.因此这里将L直接抽象为每一时期SARS患者的就诊率)(t的函数.平均传染期L应随)(t的变化而变化.但是在初期，由于政府对SARS的认识不足，并没有采取有效控制措施，L的变化很小可以近似看作定值，这里我们取SARS病毒最长潜伏期（约19天）为这个定值；在爆发期，有效控制措施的逐步加强，使SARS患者的就诊率)(t逐渐增加，而平均传染期L会逐渐减小并趋于一个定值，这里我们将SARS病毒平均潜伏期（约7天）定为L的最小值；在此后的高峰期以及衰减期，由于控制措施都保持在一定水平，L的值会维持在7天左右.5．针对北京疫情求解模型首先采用数学推导的方法，确定参数和v，并证明模型有唯一解.○1确定和v的关系令v，方程组中)1()2(得：SdSdI11在病情刚开始时，011SdSdI，由于)(tS是单调减少的，且)(tI最终趋近于0，则当1S时，)(tI单调减少趋近于0；当1S时，)(tI先单调增加达到最大值，然后单调减少趋近于0.容易知道，当1S时，才满足SARS的传播规律，所以参数和v的取值必须满足这个条件.○2证明模型有唯一解在初值条件下解微分方程组：111000RSISdSdI得到关系式：)ln(11)(00SSSRtI3令t，由○1得)ln(11000SSSR因为0S，所以令)ln(11)(00SxxRxf则)(lim0xfx，01)(0000ISRSf当10S时，由于0)(xf在),0(0S范围内有根，因而在)1,0(内有根.当10S时，因为xxxf1)('当1x时，0)('xf，所以0)()1(00ISff，因而0)(xf在)1,0(内也有根.注意到当10x时，0)('xf，故0)(xf在)1,0(内有唯一根.所以，S在)1,0(内有唯一解.○3划分SARS传播的4个阶段由于SARS的传播经历了4个阶段，所以，要以具体的指标划分这4个阶段.因为在4个阶段中，日发病率)()()(tItSNt是一个区分每个阶段特点的关键特征，所以以日发病率作为划分的指标.从第一个患者出现日开始：征兆期：日发病率在10（人/天）以下.北京疫情期的前40天.爆发期：从日发病率10（人/天）到日发病率最大，即0dtd时.北京疫情期的第40天到第74天.高峰期：从日发病率最大到患者数量最大，即0dtdI时.北京疫情期的第74天到第79天.衰退期：患者数量最大点以后.北京疫情期第79天以后.○4确定和v根据北京最终SARS患者总数2521人以及北京人口总数（约14000000人），得19998.01400000025211S，所以11v.因为平均传染期vL1，而L是SARS患者就诊率)(t的函数，且]19,7[L，所以，这里设计L函数为：)(17teL)(t由政府的控制措施决定，它的变化反映了政府控制措施的力度.根据实际情况，推导出：474t174t40)178.340(log40t00)(10tt而接触率与全社会的警觉程度和公众采取的各种措施有关，根据实际情况确定为：79t0672.079t7433ln116.074t403400126.040t0126.0tt确定出所有的参数后，做出北京各时期累计全社会SARS患者数和各时期累计确诊SARS患者数预测图（图1）以及北京市预测确诊SARS患者累计和实际诊SARS患者累计对比图（图2）.同时得到：北京SARS疫情的预测持续时间为106天，预测SARS患者累计2514人.图1北京市预测非典病人累计总数和预测非典病人确诊病例累计对比图5图2北京市预测确诊病例累计和实际确诊病例累计对比图北京市预测确诊病例累计和实际确诊病例累计的程序代码functionf=sirS(1)=14000000;I(1)=1;R(1)=0;na=0.126;F=19;L=19;JU=19;M(1)=1;fori=2:74%初期与爆发期ifi=40&i74JU=JU-0.25;endifi=40na=na-0.01/35;%爆发期缓慢减少endS(i)=S(i-1)-na*S(i-1)*I(i-1)/14000000;%求解S，I，RifiL+2R(i)=S(i-L-2)-S(i-L-1);elseR(i)=0;endifi=51F=F-0.5;L=fix(F);ifF==LR(i)=S(i-L-3)-S(i-L-1);endendI(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'),holdon,plot(i+JU-19,t,'sr'),ifp0.1plot(i+JU-19,p,'or'),plot(i+JU-19,p,'or')endend6h=na;g(1)=17;n=1;fori=75:139%高峰期与衰减期n=n+1;ifi80na=h-log(i-73)/35;%在高峰期急剧减少，此后为一定值endR(i)=S(i-L-2)-S(i-L-1);%求解S，I，RS(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