傳染病控制模型

1、傳染病控制模型封面?zhèn)魅静鞑タ刂颇Ｐ驼?本文建立了生物動(dòng)力學(xué)中的傳染病動(dòng)力學(xué)模型 , 把本研究系統(tǒng)的人群分 為七類(lèi)人群 : 易感人群，潛伏者人群，自由帶菌者人群，疑似人群，確診人群，治 愈人群，死亡人群。并且在分析這七類(lèi)人群的傳染病傳播網(wǎng)絡(luò)圖的基礎(chǔ)上 , 建立了 六個(gè)微分方程 , 通過(guò)這一組微分方程 ,分析求解傳染病傳播的控制狀態(tài) .并對(duì)模型進(jìn) 行了數(shù)據(jù)模擬和靈敏度分析 , 分析結(jié)果與實(shí)際情況比較符合 .在問(wèn)題二,三,四的求解中 ,我們注意到當(dāng)初始的傳染病潛伏者的基數(shù)不同時(shí) ,潛 伏者與易感人群的有效接觸率是不同的 ,因此, 我們分別引入不同的有效接觸率計(jì)算 公式 , 并通過(guò)作圖進(jìn)行比較分析

2、 .關(guān)鍵詞 傳染病控制 傳染病動(dòng)力學(xué) 微分方程 靈敏度分析 數(shù)據(jù)模擬一、問(wèn)題提出 已知某種不完全確知的具有傳染性病毒的潛伏期為 d1d2 天，病患者的治愈時(shí) 間為 d3 天。該病毒可通過(guò)直接接觸、口腔飛沫進(jìn)行傳播、擴(kuò)散，該人群的人均每 天接觸人數(shù)為 r 。為了控制病毒的擴(kuò)散與傳播將該人群分為五類(lèi) : 確診患者、疑似 患者、治愈者、死亡和正常人，可控制參數(shù)是隔離措施強(qiáng)度 p( 潛伏期內(nèi)的患者被 隔離的百分?jǐn)?shù) ) 。 要求:1( 在合理的假設(shè)下建立該病毒擴(kuò)散與傳播的控制模型 ;2( 利用所建立的模型針對(duì)如下數(shù)據(jù)進(jìn)行模擬條件 1:d1=1, d2=14, d3=30, r=20,條件 2: 已經(jīng)知道

3、的初始發(fā)病人數(shù)為 900、疑似患者為 2000 條件 3: 隔離措施強(qiáng)度 p=60%條件 4: 患者 2天后入院治療，疑似患者 2天后被隔離 ,試給出患者人數(shù)隨時(shí)間 變化的曲線圖，并明確標(biāo)識(shí)圖中的一些特殊點(diǎn)的具體數(shù)據(jù)，分析結(jié)果的合理性。 3. 若將 2 中的條件 4 改為條件 : 患者 1.5 天后入院治療，疑似患者 1.5 天后被隔離 ,模 擬結(jié)果有何變化 ,4(若僅將2中的條件3改為條件：隔離措施強(qiáng)度p=40%模擬結(jié)果有何變化,二、問(wèn)題分析分析一： 該問(wèn)題屬于生物動(dòng)力學(xué)中的傳染病動(dòng)力學(xué)模型，我們將在前人有關(guān)1,2,3 的研究的基礎(chǔ)上，將該模型進(jìn)一步的改造與拓展，使得題目中的 問(wèn)題得以解決

4、;分析二 ： 針對(duì)問(wèn)題一，我們考慮到病患者在潛伏期階段具有較小的病毒傳染能 力，病患者治愈后，大部分具有免疫能力，有免疫能力的治愈者將移出傳染病的傳 播網(wǎng)絡(luò)圖，小部分治愈者會(huì)二次感染，再次成為病毒感染者。分析三 ： 我們根據(jù)問(wèn)題，假設(shè)該地區(qū)的總?cè)丝诓蛔?，把該地區(qū)的人群分為七類(lèi) 人群 ： 易感人群，潛伏者人群，自由帶菌者人群，疑似人群，確診人群，治愈人 群，死亡人群。這七類(lèi)人群的傳染病傳播網(wǎng)絡(luò)圖(圖 1)如下：k 8kp(1), 6Kt()1St() Et() It()pk 3k 7k2k1Pt()k9Kt()kkRt()Dt() 2108圖 11三、模型假設(shè)假設(shè)一: 自由帶菌者，潛伏者人群有都

5、具傳播病毒的能力 ;假設(shè)二: 隔離人群完全斷絕與外界接觸，不再具有傳染性 ;假設(shè)三:自由帶菌者一經(jīng)發(fā)現(xiàn)就被確定為確診病人，而不必隔離觀察 ; 而被自由 帶菌者感染的人歸為潛伏者人群，未被發(fā)現(xiàn)則歸為自由帶菌者 ;假設(shè)四: 不考慮潛伏者人群，疑似人群的死亡情況，只考慮自由帶菌者人群確 診人群的死亡情況 ; 和假設(shè)五: 該地區(qū)的總?cè)丝诓蛔?，該地區(qū)的自然死亡率等于該地區(qū)的出生率 ;假設(shè)六: 病患者治愈后，大部分具有免疫能力，有免疫能力的治愈者將移出傳 染病的傳播網(wǎng)絡(luò)圖，小部分治愈者會(huì)二次感染，再次成為病毒感染者。四、符號(hào)說(shuō)明 ( 列表)第天易感人群人數(shù) S(t)t第天潛伏人群人數(shù) Et()tIt()

6、第天自由攜帶病毒人群人數(shù) tPt() 第天疑似人群人數(shù) tDt() 第天確診病人人數(shù) tRt() 第天痊愈者人數(shù) tKt() 第天自由病毒攜帶者中因病死亡人數(shù) t1Kt() 第天確診病人中因病死亡人數(shù) t2kk各過(guò)程比例系數(shù)110,k 潛伏者感染能力系數(shù) (也即) 14,k 自由病毒攜帶者感染能力系數(shù) (也即) 25r 人群中每人日接觸人數(shù)控制常數(shù) p五、模型建立與求解5.1 模型建立為了考慮隔離對(duì)傳播的影響，我們將人口分為易感者類(lèi)、潛伏者類(lèi)、自由攜帶 病毒者、疑似病人、確診病人、死亡者、痊愈者，這幾類(lèi)人相互轉(zhuǎn)化?？紤]易感人群 : 部分易感者被潛伏者和自由攜帶病毒者感染成潛伏者，痊愈者痊愈后回

7、到易感人群，部分易感人群因癥狀與該病癥狀相似而被當(dāng)做疑似病 人隔離離開(kāi)易感人群，部分疑似病人被證明并非該病患者后會(huì)解除隔離回到易感人 群，所以易感人群的人數(shù)變化的微分方程為 :dStSt()(),，,kRtkPtkPtEtItr()()()()() (1) 12312dtStEtItRt()()()()St(),EtIt()()，其中表示實(shí)際具有感染能力的人數(shù)，表12StEtItRt()()()()，示健康者在可自由活動(dòng)的人群中的分布，因?yàn)榭紤]到具有感染能力的人群每人日rr 接觸人數(shù)也為，但是人并非全是健康人，其中也可能有自由病毒攜帶者和潛伏者，具有感染能力的人的接觸并不能導(dǎo)致健康人群的減少，

8、只有當(dāng)具有感染能力的人與健康者接觸才考慮健康者人數(shù)的減少。不過(guò)該分布只符合初始時(shí)潛伏人群和自由病毒攜帶者的人數(shù)相對(duì)健康人群 （易感人群 ）比例較大的情況，如果潛伏人群和自由病毒攜帶者比較少的時(shí)候，該分布就接近 1，是不符合的。因?yàn)榭紤]到初始時(shí)潛伏人群和自由病毒攜帶者比較少的時(shí)候，它們不是均勻分布在自由人 群中的?？紤]潛伏人群 : 部分易感人群被感染導(dǎo)致潛伏人群增加，潛伏者被隔離導(dǎo)致潛 伏者減少，潛伏者變?yōu)椴《咀杂蓴y帶者導(dǎo)致潛伏者減少，所以潛伏者的人數(shù)變化的 微分方程為 :dEtSt()(), ，,()()()(1)()EtItrpEtkpEt (2) 126 ， dtStEtItRt()()(

9、)()kpEt(1)(), 其中表示潛伏者轉(zhuǎn)化為病毒自由攜帶者的人數(shù)。 6考慮疑似人群 : 部分潛伏人群被隔離到疑似人群，部分易感人群因?yàn)榘Y狀與該 病癥狀相似也被隔離到疑似人群，部分疑似人群被確定并非由該病引起癥狀后會(huì)被 放回到易感人群，部分疑似人群會(huì)被確診后被移到確診病人隔離區(qū)，所以疑似人群 人數(shù)變化的微分方程為 :dPt() (3) , ，,pEtkStkPtkPt()()()()329dt9287kPt()kPt() 其中表示為誤診的正常人，表示轉(zhuǎn)為確診病人的疑似病人數(shù)?？紤]自由病毒攜帶者 : 部分潛伏者轉(zhuǎn)化為自由病毒攜帶者，部分自由病毒攜帶 者被確診后隔離，部分病毒自由攜帶者死亡，所以

10、自由病毒攜帶者的人數(shù)變化微分 方程為:dIt() (4) ,kpEtkItkIt(1)()()()687dtkIt()kIt() 其中表示死亡的自由病毒攜帶者人數(shù)，表示確診后被隔離的自 由病毒攜帶者人數(shù)?？紤]確證病人數(shù) : 部分疑似病人確診為病人，確證的病人要么治療后痊愈要么 治療無(wú)效后死亡，所以確診病人數(shù)變化的微分方程為 :3dDt() (5) , ，,kPtkItkDtkRt()()()()97810dtkRt()kDt() 其中，表示確證后治療無(wú)效死亡的人數(shù)，為去確診后痊愈的人數(shù)?？紤]痊愈者 : 痊愈者痊愈后回到易感人群，確診的病人治療后成為痊愈人員， 所以痊愈者的人數(shù)變化的微分方程為

11、:dRt() (6) ,kRtkRt()()101dt綜上，該傳染病傳染過(guò)程中各狀態(tài)量的微分方程組如下 :dStSt()(),，, ， ,kRtkPtkPtEtItr()()()()()12312,dtStEtItRt()()()() ， ,dEtSt()(),，,()()()(1)()EtItrpEtkpEt,126,，dtStEtItRt()()()(),dPt(), ，,pEtkStkPtkPt()()()()329,dt(7) , dIt(),kpEtkItkIt(1)()()()687,dt,dDt(),，,kPtkItkDtkRt()()()()97810,dt,dRt(),kRt

12、kRt()()101dt,5.2 模型求解下面采用上述模型 (7) 對(duì)本題進(jìn)行求解。首先確定各比例系數(shù)。 考慮痊愈者恢復(fù)后回到易感人群后部分痊愈者已經(jīng)具有抵抗該染病的能力，不 會(huì)再次被感染，即不屬于易感人群，可以算作退出了該病的傳播系統(tǒng)，但仍有k,0.1 部分痊愈者會(huì)再次被感染，在這里我們賦值 1 被誤診的疑似病人被確定誤診需要經(jīng)過(guò)一段時(shí)間的觀察才能確定誤。診，考 dd-慮到潛伏期是從天，且疑似人群部分直接來(lái)自易感人群，所以賦值12kdd, ， 0.6431/() 。 212 據(jù)國(guó)家公布的數(shù)據(jù)計(jì)算，當(dāng)增加一個(gè)確診病人時(shí)大約有 1.3 個(gè)新增疑似病人來(lái) 自易感人群和潛伏人群，而來(lái)自易感人群的比例

13、為63.41%，又由于確診病人被隔離時(shí)間30天才可治愈，故平均每天新增確診人數(shù)為D/30，所以我們賦值kd,1.3*0.6431/ 。 33k,0.1kk 表示潛伏者的感染能力，在這里我們賦值，表示自由病毒攜帶 445k,0.8 者的感染能力，在這里我們賦值，我們賦值不同的理由是考慮到潛伏階54 段，身體內(nèi)病毒數(shù)量低于自由攜帶病毒階段，潛伏階段感染能力相對(duì)較弱。 潛伏者中一部分被隔離，另一部分將在自由環(huán)境中生活，最后發(fā)展成為確診 ()/2dd ，病人，我們?cè)O(shè)平均潛伏期為。根據(jù)國(guó)家公布數(shù)據(jù)，每天報(bào)告的確診病12例中有 80%是從疑似病轉(zhuǎn)來(lái)的，另外 20%是新收治的，因此我們賦值， k8=0.14

14、/d3 ，k9=0.3659/3 。 k6=(2/(d1+d2)*(20/100)k7=1/3治愈者治愈時(shí)間是d3,所以我們可賦值k10=1/d3。綜上即為參數(shù)的賦值，下面我們利用題目中條件，應(yīng)用MATLAB件求解模型。由于題目沒(méi)有給出初始時(shí)潛伏人群和自由病毒攜帶者的人數(shù)，本文就根據(jù)模型 作出假設(shè) :初始易感人群人數(shù)為 1200000，潛伏者人數(shù)為 100000，自由病毒攜帶者 人數(shù)為 40000。這是相對(duì)符合潛伏者和自由病毒攜帶者人數(shù)相對(duì)易感人群比例大的 情況。5.2.1 問(wèn)題二的求解ddd,1,14,30 當(dāng)，初始發(fā)病人數(shù)為 900，疑似發(fā)病人數(shù)為 2000， r,20123 患者兩天后住

15、院治療，疑似患者 2天后被隔離，隔離措施強(qiáng)度 P=60%時(shí)圖像圖 2 如下，考慮到疑似患者 2 天后被隔離， k9=0.3659/1 ，病人治愈時(shí)間也因?yàn)樽≡?推d,32 遲而導(dǎo)致治愈期也會(huì)推遲 2 天， : 310065ioyBE圖2522問(wèn)題三的求解ddd,1,14,30 當(dāng)，初始發(fā)病人數(shù)為 900,疑似發(fā)病人數(shù)為2000, r,20123患者1.5天后住院治療，疑似患者1.5天后被隔離，考慮到疑似患者 2天后被隔離，k9=0.3659/1.5，病人治愈時(shí)間也因?yàn)樽≡和七t而導(dǎo)致治愈期也會(huì)推遲1.5天，d,31.5隔離措施強(qiáng)度P=60%寸圖像（圖3）如下：3523問(wèn)題四的求解ddd,1,14

16、,30 當(dāng)，初始發(fā)病人數(shù)為 900,疑似發(fā)病人數(shù)為2000, r,20123患者兩天后住院治療，疑似患者 2天后被隔離，考慮到疑似患者2天后被隔離，d,32k9=0.3659/1，病人治愈時(shí)間也因?yàn)樽≡和七t而導(dǎo)致治愈期也會(huì)推遲2天,3隔離措施強(qiáng)度P=40%時(shí)圖像(圖4)如下:10010D圖4程序見(jiàn)附錄程序一。5.3進(jìn)一步分析求解因?yàn)樯鲜瞿Ｐ椭皇欠蠞摲吆妥杂刹《緮y帶者人數(shù)相對(duì)易感人群比例大的情 況。但在實(shí)際中潛伏者和自由病毒攜帶者人數(shù)相對(duì)是比較少的，所以就對(duì)模3型進(jìn)行了一下的優(yōu)化。我們通過(guò)相關(guān)資料查到有效接觸率是服從指數(shù)分布 的而所以易感人群的人數(shù)變化的微分方程為：6dSt() (9) ,

17、, , ,kRtkPtkPtEtltrt()()()()()()12312dt然后對(duì)模型求解，其比例系數(shù)不變，現(xiàn)在作出的假設(shè)潛伏者和自由病毒攜帶者 人數(shù)可以相對(duì)少一些。假設(shè)初始易感人群人數(shù)為100000,潛伏者人數(shù)為2000,自由病毒攜帶者人數(shù)為800。5.3.1問(wèn)題二的求解ddd,1,14,30 當(dāng)，初始發(fā)病人數(shù)為900,疑似發(fā)病人數(shù)為r,201232000,患者兩天后住院治療，疑似患者 2天后被隔離，隔離措施強(qiáng)度 P=60%時(shí)圖像(圖5)如下，考慮到疑似患者2天后被隔離，k9=0.3659/1，病人治愈時(shí)間也圖55.3.2問(wèn)題三的求解ddd,1,14,30 當(dāng)，初始發(fā)病人數(shù)為 900,疑似

18、發(fā)病人數(shù)為2000, r,20123患者1.5天后住院治療，疑似患者1.5天后被隔離，考慮到疑似患者 2天后被 隔離，k9=0.3659/1.5，病人治愈時(shí)間也因?yàn)樽≡和七t而導(dǎo)致治愈期也會(huì)推遲 1.5 天，d,31.5隔離措施強(qiáng)度P=60%寸圖像（圖6）如下：375.3.3問(wèn)題四的求解ddd,1,14,30 當(dāng)，初始發(fā)病人數(shù)為 900,疑似發(fā)病人數(shù)為2000，r,20123 患者兩天后住院治療，疑似患者 2天后被隔離，考慮到疑似患者2天后被隔 離，d,32k9=0.3659/1，病人治愈時(shí)間也因?yàn)樽≡和七t而導(dǎo)致治愈期也會(huì)推遲2天，3隔離措施強(qiáng)度P=40%寸圖像（圖7）如下:程序見(jiàn)附錄程序二六、

19、模型分析6.1模型穩(wěn)定性分析6.2模型靈敏度分析當(dāng)隔離措施強(qiáng)度增大時(shí)，在其它條件不變的情況下，發(fā)病的高峰就會(huì)減少，同 時(shí)高峰也會(huì)提前到達(dá)。當(dāng)對(duì)患者提前治療和對(duì)疑似患者提前隔離時(shí)，在其它條件不變的情況下，發(fā)病 的高峰會(huì)相對(duì)減少，同時(shí)高峰會(huì)推遲到達(dá)。所以，隨著時(shí)間的推移，對(duì)病情進(jìn)一步的了解，適當(dāng)增加隔離措施強(qiáng)度，還有 對(duì)患者提前治療和對(duì)疑似患者提前隔離，這樣就可以更加有效地控制病情。七、參考文獻(xiàn).1 徐翠翠,赫孝良.一類(lèi)具有潛伏期和隔離項(xiàng)的SEIQR流行病模型J.商丘師 范學(xué)院學(xué)報(bào)，2010,26 (6)2 韓曉娜.SARS流行病學(xué)傳播動(dòng)力模型研究.中國(guó)人民解放軍軍事醫(yī)學(xué)科學(xué)院 碩士學(xué)位論文，20

20、06.53 周義倉(cāng), 赫孝良.數(shù)學(xué)建模實(shí)驗(yàn) M. 西安: 西安交通大學(xué)出版社， 2007.8八、附錄程序一function dy=chuanranbing(t,y)d1=1;d2=14;d3=32;8r=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(

21、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=