传染病控制模型

传染病控制模型 封面 传染病传播控制模型 摘要 本文建立了生物动力学中的传染病动力学模型, 把本研究系统的人群分为七类人群:易感人群，潜伏者人群，自由带菌者人群，疑似人群，确诊人群，治愈人群，死亡人群。并且在分析这七类人群的传染病传播网络图的基础上,建立了六个微分方程,通过这一组微分方程,分析求解传染病传播的控制状态.并对模型进行了数据模拟和灵敏度分析,分析结果与实际情况比较符合. 在问题二,三,四的求解中,我们注意到当初始的传染病潜伏者的基数不同时,潜伏者与易感人群的有效接触率是不同的,因此,我们分别引入不同的有效接触率计算,,并通过作图进行比较分析. 关键词 传染病控制 传染病动力学 微分方程 灵敏度分析 数据模拟 一、问题提出 已知某种不完全确知的具有传染性病毒的潜伏期为d1~d2天，病患者的治愈时间为d3天。该病毒可通过直接接触、口腔飞沫进行传播、扩散，该人群的人均每天接触人数为r。为了控制病毒的扩散与传播将该人群分为五类:确诊患者、疑似患者、治愈者、死亡和正常人，可控制参数是隔离强度p(潜伏期内的患者被隔离的百分数)。 要求: 1(在合理的假设下建立该病毒扩散与传播的控制模型; 2(利用所建立的模型针对如下数据进行模拟 条件1:d1=1, d2=14, d3=30, r=20, 条件2:已经知道的初始发病人数为900、疑似患者为2000 条件3:隔离措施强度p=60% 条件4:患者2天后入院治疗，疑似患者2天后被隔离,试给出患者人数随时间变化的曲线图，并明确标识图中的一些特殊点的具体数据，分析结果的合理性。 3.若将2中的条件4改为条件:患者1.5天后入院治疗，疑似患者1.5天后被隔离,模拟结果有何变化, 4(若仅将2中的条件3改为条件:隔离措施强度p=40%，模拟结果有何变化, 二、问题分析 分析一:该问题属于生物动力学中的传染病动力学模型，我们将在前人有关 [1],[2],[3]的研究的基础上，将该模型进一步的改造与拓展，使得题目中的问题得以解决; 分析二:针对问题一，我们考虑到病患者在潜伏期阶段具有较小的病毒传染能力，病患者治愈后，大部分具有免疫能力，有免疫能力的治愈者将移出传染病的传播网络图，小部分治愈者会二次感染，再次成为病毒感染者。 分析三:我们根据问题，假设该地区的总人口不变，把该地区的人群分为七类人群:易感人群，潜伏者人群，自由带菌者人群，疑似人群，确诊人群，治愈人群，死亡人群。这七类人群的传染病传播网络图(图1)如下: k 8kp(1), 6Kt()1St() Et() It() pk 3 k 7k2k1 Pt() k9 Kt()kkRt()Dt() 2108 图 1 1 三、模型假设 假设一:自由带菌者，潜伏者人群有都具传播病毒的能力; 假设二:隔离人群完全断绝与外界接触，不再具有传染性; 假设三:自由带菌者一经发现就被确定为确诊病人，而不必隔离观察;而被自由带菌者感染的人归为潜伏者人群，未被发现则归为自由带菌者; 假设四:不考虑潜伏者人群，疑似人群的死亡情况，只考虑自由带菌者人群确诊人群的死亡情况; 和 假设五:该地区的总人口不变，该地区的自然死亡率等于该地区的出生率; 假设六:病患者治愈后，大部分具有免疫能力，有免疫能力的治愈者将移出传染病的传播网络图，小部分治愈者会二次感染，再次成为病毒感染者。 四、符号说明(列表) —第天易感人群人数 S(t)t —第天潜伏人群人数 Et()t It()—第天自由携带病毒人群人数 t Pt()—第天疑似人群人数 t Dt()—第天确诊病人人数 t Rt()—第天痊愈者人数 t Kt()—第天自由病毒携带者中因病死亡人数 t1 Kt()—第天确诊病人中因病死亡人数 t2 kk~—各过程比例系数 110 ,k—潜伏者感染能力系数(也即) 14 ,k—自由病毒携带者感染能力系数(也即) 25 r—人群中每人日接触人数 —控制常数 p 五、模型建立与求解 5.1 模型建立 为了考虑隔离对传播的影响，我们将人口分为易感者类、潜伏者类、自由携带病毒者、疑似病人、确诊病人、死亡者、痊愈者，这几类人相互转化。 考虑易感人群:部分易感者被潜伏者和自由携带病毒者感染成潜伏者，痊愈 2 者痊愈后回到易感人群，部分易感人群因症状与该病症状相似而被当做疑似病人隔离离开易感人群，部分疑似病人被证明并非该病患者后会解除隔离回到易感人群，所以易感人群的人数变化的微分方程为: dStSt()(),，,,，,,kRtkPtkPtEtItr()()()[()()] (1) 12312，，，dtStEtItRt()()()() St(),,EtIt()()，其中表示实际具有感染能力的人数，表12StEtItRt()()()()，，，示健康者在可自由活动的人群中的分布，因为考虑到具有感染能力的人群每人日 rr接触人数也为，但是人并非全是健康人，其中也可能有自由病毒携带者和潜伏者，具有感染能力的人的接触并不能导致健康人群的减少，只有当具有感染能力的人与健康者接触才考虑健康者人数的减少。不过该分布只符合初始时潜伏人群和自由病毒携带者的人数相对健康人群(易感人群)比例较大的情况，如果潜伏人群和自由病毒携带者比较少的时候，该分布就接近1，是不符合的。因为考虑到初始时潜伏人群和自由病毒携带者比较少的时候，它们不是均匀分布在自由人 群中的。 考虑潜伏人群:部分易感人群被感染导致潜伏人群增加，潜伏者被隔离导致潜伏者减少，潜伏者变为病毒自由携带者导致潜伏者减少，所以潜伏者的人数变化的微分方程为: dEtSt()(),，,,,,,[()()]()(1)()EtItrpEtkpEt (2) 126，，，dtStEtItRt()()()() kpEt(1)(),其中表示潜伏者转化为病毒自由携带者的人数。 6 考虑疑似人群:部分潜伏人群被隔离到疑似人群，部分易感人群因为症状与该病症状相似也被隔离到疑似人群，部分疑似人群被确定并非由该病引起症状后会被放回到易感人群，部分疑似人群会被确诊后被移到确诊病人隔离区，所以疑似人群人数变化的微分方程为: dPt() (3) ,，,,pEtkStkPtkPt()()()()329dt kPt()kPt()其中表示为误诊的正常人，表示转为确诊病人的疑似病人数。 92 考虑自由病毒携带者:部分潜伏者转化为自由病毒携带者，部分自由病毒携带者被确诊后隔离，部分病毒自由携带者死亡，所以自由病毒携带者的人数变化微分方程为: dIt() (4) ,,,,kpEtkItkIt(1)()()()687dt kIt()kIt()其中表示死亡的自由病毒携带者人数，表示确诊后被隔离的自87 由病毒携带者人数。 考虑确证病人数:部分疑似病人确诊为病人，确证的病人要么治疗后痊愈要么治疗无效后死亡，所以确诊病人数变化的微分方程为: 3 dDt() (5) ,，,,kPtkItkDtkRt()()()()97810dt kRt()kDt()其中，表示确证后治疗无效死亡的人数，为去确诊后痊愈的人108 数。 考虑痊愈者:痊愈者痊愈后回到易感人群，确诊的病人治疗后成为痊愈人员，所以痊愈者的人数变化的微分方程为: dRt() (6) ,,kRtkRt()()101dt 综上，该传染病传染过程中各状态量的微分方程组如下: dStSt()(),,，,,，,,kRtkPtkPtEtItr()()()[()()]12312,dtStEtItRt()()()()，，，, dEtSt()(),,，,,,[()()]()(1)()EtItrpEtkpEt,,126,，，，dtStEtItRt()()()(), dPt(),,，,,pEtkStkPtkPt()()()()329,dt(7) , dIt(),,,,,kpEtkItkIt(1)()()()687,dt,dDt(),,，,,kPtkItkDtkRt()()()()97810,dt ,dRt(),,,kRtkRt()()101dt, 5.2 模型求解 下面采用上述模型(7)对本题进行求解。 首先确定各比例系数。 考虑痊愈者恢复后回到易感人群后部分痊愈者已经具有抵抗该染病的能力，不会再次被感染，即不属于易感人群，可以算作退出了该病的传播系统，但仍有 k,0.1部分痊愈者会再次被感染，在这里我们赋值 1 被误诊的疑似病人被确定误诊需要经过一段时间的观察才能确定误。诊，考 dd~虑到潜伏期是从天，且疑似人群部分直接来自易感人群，所以赋值12 kdd,，0.6431/()。 212 据国家公布的数据计算，当增加一个确诊病人时大约有1.3个新增疑似病人来自易感人群和潜伏人群，而来自易感人群的比例为63.41%，又由于确诊病人被隔离时间30天才可治愈，故平均每天新增确诊人数为D/30，所以我们赋值kd,1.3*0.6431/。 33 k,0.1kk表示潜伏者的感染能力，在这里我们赋值，表示自由病毒携带445 k,0.8者的感染能力，在这里我们赋值，我们赋值不同的理由是考虑到潜伏阶5 4 段，身体内病毒数量低于自由携带病毒阶段，潜伏阶段感染能力相对较弱。 潜伏者中一部分被隔离，另一部分将在自由环境中生活，最后发展成为确诊 ()/2dd，病人，我们设平均潜伏期为。根据国家公布数据，每天的确诊病12 例中有80%是从疑似病转来的，另外20%是新收治的，因此我们赋值 ，，k8=0.14/d3，k9=0.3659/3。 k6=(2/(d1+d2))*(20/100)k7=1/3 治愈者治愈时间是d3，所以我们可赋值k10=1/d3。 综上即为参数的赋值，下面我们利用题目中条件，应用MATLAB软件求解模型。 由于题目没有给出初始时潜伏人群和自由病毒携带者的人数，本文就根据模型作出假设:初始易感人群人数为1200000，潜伏者人数为100000，自由病毒携带者人数为40000。这是相对符合潜伏者和自由病毒携带者人数相对易感人群比例大的情况。 5.2.1 问题二的求解 ddd,,,1,14,30当，，初始发病人数为900，疑似发病人数为2000，r,20123 患者两天后住院治疗，疑似患者2天后被隔离，隔离措施强度P=60%时图像图 2 如下，考虑到疑似患者2天后被隔离，k9=0.3659/1，病人治愈时间也因为住院推 d,32迟而导致治愈期也会推迟2天，: 3 图 2 5.2.2 问题三的求解 ddd,,,1,14,30当，，初始发病人数为900，疑似发病人数为2000，r,20123 患者1.5天后住院治疗，疑似患者1.5天后被隔离，考虑到疑似患者2天后被隔离，k9=0.3659/1.5，病人治愈时间也因为住院推迟而导致治愈期也会推迟1.5天，d,31.5隔离措施强度P=60%时图像(图3)如下: 3 5 图 3 5.2.3 问题四的求解 ddd,,,1,14,30当，，初始发病人数为900，疑似发病人数为2000，r,20123 患者两天后住院治疗，疑似患者2天后被隔离，考虑到疑似患者2天后被隔离， d,32k9=0.3659/1，病人治愈时间也因为住院推迟而导致治愈期也会推迟2天，， 3 隔离措施强度P=40%时图像(图4)如下: 图 4 程序见附录程序一。 5.3进一步分析求解 因为上述模型只是符合潜伏者和自由病毒携带者人数相对易感人群比例大的情况。但在实际中潜伏者和自由病毒携带者人数相对是比较少的，所以就对模 [3]型进行了一下的优化。我们通过相关查到有效接触率是服从指数分布的而 ,0.1303t,(t),0.002，0.2496e且是递减的，其公式为: (8) 所以易感人群的人数变化的微分方程为: 6 dSt() (9) ,，,,，,,,kRtkPtkPtEtItrt()()()[()()]()12312dt 然后对模型求解，其比例系数不变，现在作出的假设潜伏者和自由病毒携带者人数可以相对少一些。假设初始易感人群人数为100000，潜伏者人数为2000，自由病毒携带者人数为800。 5.3.1 问题二的求解 ddd,,,1,14,30当，，初始发病人数为900，疑似发病人数为r,20123 2000，患者两天后住院治疗，疑似患者2天后被隔离，隔离措施强度P=60%时图像(图5)如下，考虑到疑似患者2天后被隔离，k9=0.3659/1，病人治愈时间也 d,32因为住院推迟而导致治愈期也会推迟2天，: 3 图5 5.3.2 问题三的求解 ddd,,,1,14,30当，，初始发病人数为900，疑似发病人数为2000，r,20123 患者1.5天后住院治疗，疑似患者1.5天后被隔离，考虑到疑似患者2天后被隔离，k9=0.3659/1.5，病人治愈时间也因为住院推迟而导致治愈期也会推迟1.5天，d,31.5隔离措施强度P=60%时图像(图6)如下: 3 图 6 7 5.3.3 问题四的求解 ddd,,,1,14,30当，，初始发病人数为900，疑似发病人数为2000，r,20123 患者两天后住院治疗，疑似患者2天后被隔离，考虑到疑似患者2天后被隔离， d,32k9=0.3659/1，病人治愈时间也因为住院推迟而导致治愈期也会推迟2天，，3隔离措施强度P=40%时图像(图 7)如下: 图7 程序见附录程序二。 六、模型分析 6.1 模型稳定性分析 6.2 模型灵敏度分析 当隔离措施强度增大时，在其它条件不变的情况下，发病的高峰就会减少，同时高峰也会提前到达。 当对患者提前治疗和对疑似患者提前隔离时，在其它条件不变的情况下，发病的高峰会相对减少，同时高峰会推迟到达。 所以，随着时间的推移，对病情进一步的了解，适当增加隔离措施强度，还有对患者提前治疗和对疑似患者提前隔离，这样就可以更加有效地控制病情。 七、参考文献. [1] 徐翠翠,赫孝良.一类具有潜伏期和隔离项的SEIQR流行病模型[J].商丘师范学院学报，2010,26 (6) [2] 韩晓娜.SARS流行病学传播动力模型研究.中国人民解放军军事医学科学院硕士学位论文，2006.5 [3] 周义仓,赫孝良.数学建模实验[M].西安:西安交通大学出版社，2007.8 八、附录 程序一 function dy=chuanranbing(t,y) d1=1; d2=14; d3=32; 8 r=20; p=0.4; k1=0.1; k2=0.6431/(d1+d2); k3=1.3*0.6431/d3; k4=0.1; k5=0.8; k6=(2/(d1+d2))*(20/100); k7=1/3; k8=0.14/d3; k9=0.3659/1; k10=1/d3; dy=zeros(6,1); dy(1)=k1*y(6)+k2*y(4)-k3*y(5)-y(1)*r*(k5*y(3)+k4*y(2))/(y(1)+y(2)+y(3)+y(6)); dy(2)=y(1)*r*(k5*y(3)+k4*y(2))/(y(1)+y(2)+y(3)+y(6))-p*y(2)-k6*(1-p)*y(2); dy(3)=k6*(1-p)*y(2)-k7*y(3)-k8*y(3); dy(4)=k3*y(5)-k2*y(4)+p*y(2)-k9*y(4); dy(5)=k9*y(4)+k7*y(3)-k10*y(5)-k8*y(5); dy(6)=k10*y(5)-k1*y(6); [T,Y]=ode45('chuanranbing',[0 360],[1200000 100000 40000 2000 900 0]); plot(T,Y(:,5),'r') text(0,900,'(0,900)','color','r') text(13.9,981226,'(13.9,981226)','color','r') text(360,310065,'(310065)','color','r') hold on 程序二 function dy=chuanranbing(t,y) d1=1; d2=14; d3=32; r=20; p=0.6; k1=0.1; k2=0.6431/(d1+d2); k3=1.3*0.6431/d3; k4=0.1; k5=0.8; k6=(2/(d1+d2))*(20/100); k7=1/3; k8=0.14/d3; k9=0.3659/1; k10=1/d3; 9 dy=zeros(6,1); dy(1)=k1*y(6)+k2*y(4)-k3*y(5)-(0.002+0.2496*exp(-0.1303*t))*r*(k5*y(3)+k4*y(2)); dy(2)=(0.002+0.2496*exp(-0.1303*t))*r*(k5*y(3)+k4*y(2))-p*y(2)-k6*(1-p)*y(2); dy(3)=k6*(1-p)*y(2)-k7*y(3)-k8*y(3); dy(4)=k3*y(5)-k2*y(4)+p*y(2)-k9*y(4); dy(5)=k9*y(4)+k7*y(3)-k10*y(5)-k8*y(5); dy(6)=k10*y(5)-k1*y(6); [T,Y]=ode45('chuanranbing',[0 360],[100000 2000 800 2000 900 0]); plot(T,Y(:,5),'r') text(0,900,'(0,900)','color','r') text(16.8,16844,'(16.8,16844)','color','r') text(360,235,'(360,235)','color','r') hold on 10