SARS传播的数学模型论文

PAGEPAGE7SARS传播的数学模型摘要本文分析了题目所提供的早期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模型的假设SARS的持续期不太长，可以忽略在SARS持续期内的城市人口的自然出生率和自然死亡率.被SARS感染后经治疗康复的人群在SARS流行期不会被再次感染.病人被严格隔离、治愈或者死亡后，不再有感染作用.不考虑人口的流动，仅仅在一个城市范围内研究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：图3SARS传播流程图在这个模型中，排除类就是已确诊SARS患者累计数，而是全社会累计SARS患者数，包括已确诊的和未被发现的两部分.4.4模型的建立有了以上的机理分析，建立起针对SARS的改进SIR模型：该模型中参数和在疫情发展的各个阶段受实际因素影响，会有比较明显的变化，现分析如下：eq\o\ac(○,1)参数表示单位时间内一个传染者与他人的接触率，其与全社会的警觉程度和政府、公众采取的各种措施有关，例如，佩戴口罩，减少停留在公共场所的时间，喷洒消毒药剂，提高隔离强度等都能有效地降低接触率的值.一般认为，的数值随着SARS发展的4个阶段不断变化.在SARS初期，由于潜伏期的存在和社会对SARS病毒传播的速度认识不足，政府和公众并未引起重视，故维持在一个较高的数值；进入爆发期后，公众发现感染者不断增加，恐慌情绪增加，随即采取多种措施，使得到一定的控制，但效果不明显，此处假设呈线性形式缓慢衰减；在高峰期，当高强度的控制措施实施后，病毒传播的有效接触率明显减少，可以认为按天数呈指数形式衰减；此后进入衰减期，就维持在一个较低值附近.eq\o\ac(○,2)参数表示传播者的被控制率.称为平均传染期，表示一个传播者在被隔离或者死亡之前具有传播能力的平均时间.一般认为，SARS患者经过传染期L过后，将隔离治疗或者死亡，从I类成员变为R类，失去传播能力.L与政府采取的措施密切相关，例如，尽量早地发现病患，对疑似病例提前进行隔离，“早发现，早隔离”；提供更广范围的医疗手段，使更多的人接受有效的治疗等，都可以有效地降低平均传染期L的长度.因此这里将L直接抽象为每一时期SARS患者的就诊率的函数.平均传染期L应随的变化而变化.但是在初期，由于政府对SARS的认识不足，并没有采取有效控制措施，L的变化很小可以近似看作定值，这里我们取SARS病毒最长潜伏期（约19天）为这个定值；在爆发期，有效控制措施的逐步加强，使SARS患者的就诊率逐渐增加，而平均传染期L会逐渐减小并趋于一个定值，这里我们将SARS病毒平均潜伏期（约7天）定为L的最小值；在此后的高峰期以及衰减期，由于控制措施都保持在一定水平，L的值会维持在7天左右.4.5针对北京疫情求解模型首先采用数学推导的方法，确定参数和，并证明模型有唯一解.eq\o\ac(○,1)确定和的关系令，方程组中得：在病情刚开始时，，由于是单调减少的，且最终趋近于0，则当时，单调减少趋近于0；当时，先单调增加达到最大值，然后单调减少趋近于0.容易知道，当时，才满足SARS的传播规律，所以参数和的取值必须满足这个条件.eq\o\ac(○,2)证明模型有唯一解在初值条件下解微分方程组：得到关系式：令，由eq\o\ac(○,1)得因为，所以令则，当时，由于在范围内有根，因而在内有根.当时，因为当时，，所以，因而在内也有根.注意到当时，，故在内有唯一根.所以，在内有唯一解.eq\o\ac(○,3)划分SARS传播的4个阶段由于SARS的传播经历了4个阶段，所以，要以具体的指标划分这4个阶段.因为在4个阶段中，日发病率是一个区分每个阶段特点的关键特征，所以以日发病率作为划分的指标.从第一个患者出现日开始：征兆期：日发病率在10（人/天）以下.北京疫情期的前40天.爆发期：从日发病率10（人/天）到日发病率最大，即时.北京疫情期的第40天到第74天.高峰期：从日发病率最大到患者数量最大，即时.北京疫情期的第74天到第79天.衰退期：患者数量最大点以后.北京疫情期第79天以后.eq\o\ac(○,4)确定和根据北京最终SARS患者总数2521人以及北京人口总数（约14000000人），得，所以.因为平均传染期，而是SARS患者就诊率的函数，且，所以，这里设计L函数为：由政府的控制措施决定，它的变化反映了政府控制措施的力度.根据实际情况，推导出：而接触率与全社会的警觉程度和公众采取的各种措施有关，根据实际情况确定为：确定出所有的参数后，做出北京各时期累计全社会SARS患者数和各时期累计确诊SARS患者数预测图（图4）以及北京市预测确诊SARS患者累计和实际确诊SARS患者累计对比图（图5）.同时得到：北京SARS疫情的预测持续时间为106天，预测SARS患者累计2514人.(计算程序见附件1：SIR模型程序)图4北京市预测非典病人累计总数和预测非典病人确诊病例累计对比图图5北京市预测确诊病例累计和实际确诊病例累计对比图附件附件1：SIR模型程序functionf=sorS(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&i<74JU=JU-0.25;endifi>=40na=na-0.01/35;%爆发期缓慢减少endS(i)=S(i-1)-na*S(i-1)*I(i-1)/14000000;%求解S，I，Rifi>L+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'),ifp>0.1plot(i+JU-19,p,'or'),plot(i+JU-19,p,'or')endendh=na;g(1)=17;n=1;fori=75:139%高峰期与衰减期n=n+1;ifi<80na=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附件2：时间序列程序functionsjxlX(1,:)=[9.411.316.819.820.318.820.924.924.724.319.418.6];X(2,:)=[9.611.715.819.919.517.817.823.321.424.520.115.9];X(3,:)=[10.112.917.721.021.020.421.925.829.329.823.616.5];X(4,:)=[11.426.019.625.927.624.323.027.827.328.532.818.5];X(5,:)=[11.526.420.426.128.928.025.230.828.728.122.220.7];X(6,:)=[13.729.723.128.929.027.426.032.231.432.629.222.9];fori=1:12y(i)=0;forj=1:6y(i)=y(i)+X(j,i);endy(i)=y(i)/6;endfori=1:6%计算、和forj=1:12Z(i,j)=X(i,j)-y(j);endendfori=1:6p(i)=i-3.5;endb=0;q=0;for