SARS傳播模型

1、SARS傳播模型摘 要：本文中我們對北京地區(qū)4月20日-6月8日的SARS疫情數(shù)據(jù)進(jìn)行了分析處理，把北京地區(qū)SARS疫情分為兩個時期：感染期（4月20日-5月16日）和恢復(fù)期（5月17日-6月8日）。由于醫(yī)務(wù)人員人群是SARS病的高發(fā)人群，所以我們在本文的模型中把病人人群分為醫(yī)務(wù)人員病人人群和非醫(yī)務(wù)人員病人人群。通過分析文中附件1的數(shù)據(jù)，我們建立了兩個時期SARS傳播的微分方程模型，并得到了模型的解，感染期模型的解為：，恢復(fù)期模型的解為：。從模型解的曲線與實際數(shù)據(jù)的比較，我們發(fā)現(xiàn)該模型的解與實際數(shù)據(jù)是非常吻合的。關(guān)鍵詞：SARS，三次樣條插值，高次多項式擬合0. 引 言SARS（Severe

2、Acute Respiratory Syndrome，嚴(yán)重急性呼吸道綜合癥, 俗稱：非典型肺炎）是一種新的、傳染性很強(qiáng)的疾病，它在我國部分地區(qū)的暴發(fā)與蔓延，嚴(yán)重威脅了人民的健康與生命安全，影響我國的社會穩(wěn)定與經(jīng)濟(jì)發(fā)展。我們從整個抗擊SARS的斗爭中得到了許多重要的經(jīng)驗和教訓(xùn)，同時也認(rèn)識到定量地研究傳染病的傳播規(guī)律、預(yù)測傳染病發(fā)展趨勢的必要性和重要性。1. 問題分析由于SARS主要是通過近距離的飛沫傳播，與病人有過密切接觸的人群就很可能被感染，成為SARS病人，所以我們在模型里假設(shè)健康者只要與病人接觸，則感染為病人，且由于醫(yī)務(wù)人員每天都直接與病人接觸，所以醫(yī)務(wù)人員人群是SARS病的高發(fā)人群，其主

3、要是被住院的病人傳染。而普通病人即非醫(yī)務(wù)人員的病人，主要是由一些病人在發(fā)病后未及時被隔離治療而與健康人接觸，并使其感染病毒，因此我們在模型里把病人分為醫(yī)務(wù)人員病人和普通病人。同時臨床統(tǒng)計數(shù)據(jù)表明SARS病的潛伏期為2-12天，一般在4-5天，治愈后的病人沒有出現(xiàn)再次患病現(xiàn)象，所以我們也假定治愈后的病人具有免疫力。雖然SARS病在2002年底到2003年初就在我國各地市廣泛傳播，但是疫情數(shù)據(jù)并沒有精確的得到統(tǒng)計，從4月20日后，國家衛(wèi)生部才在專門的網(wǎng)站上發(fā)布各地的SARS疫情數(shù)據(jù)，本文中我們只對北京地區(qū)4月20日-6月8日的數(shù)據(jù)進(jìn)行分析。通過對本文附錄1中數(shù)據(jù)的處理，我們發(fā)現(xiàn)北京市4月20日-6

4、月8日的疫情可分為兩個時期：感染期（4月20日-5月16日）和恢復(fù)期（5月17日-6月8日）（見附錄2）。并且通過擬合附件1中所給的數(shù)據(jù)，我們得到了兩個時期內(nèi)醫(yī)務(wù)人員及普通人員的日感染率k，日死亡率d和日治愈率z，建立了SARS傳染病的一個微分方程模型。2. 模型假設(shè)1) 本文主要考察的是北京地區(qū)的疫情變化, 在疾病傳播期內(nèi)該地區(qū)的人口總數(shù)N=10000000不變，即為常量。2) 政府發(fā)布的北京地區(qū)的疫情數(shù)據(jù)真實可信。3) SARS傳播途徑都視為與病源的直接接觸, 每個病人與健康者接觸時，都使健康者感染變成病人,且病人發(fā)病后馬上被隔離，即住院治療 。4) 將人群分為三類：健康者（易受感染者）、

5、病人（已被感染者）、退出者（死亡者和治愈者）,在t時刻三者在總?cè)藬?shù)中的比例分別為s(t)、i(t)、r(t)，且s(t)+i(t)+r(t)=1。其中病人又分為醫(yī)務(wù)人員病人和非醫(yī)務(wù)人員即普通病人。5) 假設(shè)不考慮SARS傳播期內(nèi)的北京地區(qū)的人口出生率和自然死亡率。對于由SARS引起的死亡人數(shù)全部歸為“退出者”, 治愈者不會再次被感染，具有免疫力，歸為“退出者”。6) 假設(shè)每日新增的普通病人與當(dāng)日健康人的人數(shù)成正比，比例系數(shù)為k1,稱為普通病人的日感染率。每日新增的醫(yī)務(wù)人員病人與當(dāng)日住院的病人人數(shù)成正比，比例系數(shù)為k2,稱為醫(yī)務(wù)人員的日感染率。每日新增的病人死亡人數(shù)與當(dāng)日已確診的病人數(shù)成正比，比

6、例系數(shù)為d,稱為病人的日死亡率。每日新增治愈的病人人數(shù)與當(dāng)日已確診的病人數(shù)成正比，比例系數(shù)為z,稱為病人的日治愈率。由于k1，k2,d,z具有很強(qiáng)的實際意義，其具體的數(shù)值我們可以通過對給定的數(shù)據(jù)擬合得到。3. 符號說明，,，i0,s0分別為病人數(shù)，健康人數(shù)占總?cè)丝跀?shù)的比例初始值。4. 模型建立從模型假設(shè)中的6)可知普通病人的日感染率k1, 醫(yī)務(wù)人員的日感染率k2, 病人的日死亡率d, 病人的日治愈率z, 我們可以通過對附錄1中的實際數(shù)據(jù)處理來得到。為了使得到的模型比較簡單，我們考慮k1,k2,d,z都是常數(shù)，這時k1、k2、d、z分別表示普通病人的日感染率、醫(yī)務(wù)人員的日感染率、病人的日死亡率、

7、病人的日治愈率的平均值，對感染期與恢復(fù)期我們可以建立下面相同的模型，這兩個時期的差別主要在于普通病人的日感染率k1是否等于0(見附錄2)。由上面的分析及假設(shè)，我們得到了下面的SARS病傳播的微分方程模型： (1)5. 模型求解51 對感染期的模型求解我們考慮4月20日至5月16日，由附錄1中所給的數(shù)據(jù)，經(jīng)過擬合得到k1=5.3527401482552×10-6，k2=0.01742890531138，d=0.00453527706374，z=0.02388800922302。且以4月20日的疫情數(shù)據(jù)作為初始值，即。對于方程組（1）我們很難得到它的解析解，但由方程組（1）的前兩個方程我

8、們消去dt，可以把（1）轉(zhuǎn)化為下面的方程組： （2）于是令，其中把s看作關(guān)于i的函數(shù)得，即 。由于，再加上初始條件s（0）、i（0），從而我們可以得到方程組（2）的解：,其中，。表示t時刻，健康人數(shù)與病人人數(shù)的比值，它反映了傳染病的嚴(yán)重程度。當(dāng)t時刻u(t)越小，即病人占總?cè)藬?shù)的比例i(t)越大，健康人數(shù)的比例s(t)越小，則SARS疫情嚴(yán)重。反之，t時刻u(t)越大，則SARS疫情得到有效控制。 從4月20日至5月16日的SARS數(shù)據(jù)可得到u(t)的變化情況(如圖1.)。從圖1.中我們看到4月20日至5月16日，的值隨著t的變化在急劇減少，主要是由于s(t)減少的同時i(t)也在快速增加，而

9、5月17日以后u(t)又開始上升，這是由于從5月17日北京市的疫情得到了有效控制，新增的SARS病例被控制在1例以內(nèi)，健康人數(shù)不再減少。相反，病人每天逐漸被治愈，i(t)也在減少。圖1.通過給模型中的解代入具體的數(shù)據(jù)，我們可以得到4月20日至5月16日病人數(shù)所占的總?cè)藬?shù)的比例i(t)的變化情況(如圖2.)。可以看出模型解與實際數(shù)據(jù)的值非常吻合，這里我們所用的u(t)是健康人數(shù)與病人數(shù)的比值的三次樣條插值函數(shù)，u(t)與實際數(shù)據(jù)越吻合則我們的解的曲線也就更加吻合于實際數(shù)據(jù)。圖2.52 對恢復(fù)期的模型求解考慮4月20日至5月16日，由附錄1中所給的數(shù)據(jù)，經(jīng)過擬合得到k1=0，k2=2.386572

10、17981×10-4，d=0.00117550039945，z=0.04137494891164。且以5月16日的疫情數(shù)據(jù)作為初始值，即。將其代入5.1的方程組（1）可得恢復(fù)期的模型： (3)通過消去方程組（3）中的dt，則方程組（3）轉(zhuǎn)化為 （4）其解為：，其中。令，則方程組（4）的解可轉(zhuǎn)化為：。我們?nèi)匀挥?.1中得到的u(t)，即健康人數(shù)與病人數(shù)的比值的三次樣條插值函數(shù)，則恢復(fù)期模型的解的圖形如下：圖3.53 北京市4月20日至6月8日SARS的疫情傳播情況由5.1和5.2的結(jié)果，我們得到了北京市4月20日至6月8日SARS的疫情傳播過程中，病人數(shù)占總?cè)藬?shù)的比例的變化趨勢圖（如圖

11、4.）。圖4.6. 模型推廣由于每一天k1,k2,d,z的值是不同的，即它們是隨著時間而變化的值，所以我們可以考慮k1,k2,d,z都是關(guān)于時間t的函數(shù)，即在t時刻醫(yī)務(wù)人員的日感染率k2(t)、普通病人的日感染率k1(t)、病人的日死亡率d(t)、病人的日治愈率z(t)，再利用實際數(shù)據(jù)經(jīng)過多項式擬合得到k1(t)、k2(t)、d(t)、z(t)的多項式函數(shù)，或者用三次樣條函數(shù)進(jìn)行插值。這樣使得模型的解與實際數(shù)據(jù)更加的吻合。同時我們也可以注意到SARS病的有效控制，很大程度上是由于傳染病人在發(fā)病后馬上被隔離，得到治療。但是現(xiàn)實生活中不可能做到這樣，在病人發(fā)病到被隔離治療之間的這段時間 L內(nèi)，病人

12、就會與健康人接觸而感染健康人。因此L的值在建立SARS數(shù)學(xué)模型中非常重要。所以本模型也可以在這些方面進(jìn)行改進(jìn)。在得到的模型的解后，我們?yōu)榱丝吹侥Ｐ偷慕馀c實際數(shù)據(jù)的接近程度，通過三次樣條函數(shù)進(jìn)行插值u(t)的值，然后代入解中來分析病人數(shù)占總?cè)藬?shù)的比例i(t)隨時間的變化情況。對于u(t)，我們也可以用高次多項式或分段高次多項式來處理，從下列圖中我們發(fā)現(xiàn)選用的多項式次數(shù)越高，模型的解的值與實際數(shù)據(jù)的接近程度就更好，但是多項式次數(shù)越高就會使數(shù)據(jù)處理量加大，并且得到的解不穩(wěn)定，出現(xiàn)病態(tài)數(shù)據(jù)。圖5.7. 參考文獻(xiàn)1 姜啟源編，數(shù)學(xué)模型，北京：高等教育出版社，2003年4月。2 Wliliam F.luc

13、as主編，微分方程模型,長沙:國防科技大學(xué)出版社,1998年5月。3 王高雄等編，常微分方程，北京：高等教育出版社，1999年3月。4 龔建華等著，SARS疫情控制的模擬分析，遙感學(xué)報，2003,7(4):260-265。5 Shoichiro Nakamura著，科學(xué)計算引論-基于MATLAB的數(shù)值分析，北京：電子工業(yè)出版社，2002年6月。6 傅鸝等編著，數(shù)學(xué)實驗，北京：科學(xué)出版社，2001年3月。7 葉其孝主編，大學(xué)生數(shù)學(xué)建模競賽輔導(dǎo)教材（一），長沙：湖南教育出版社，1998年4月。附錄1：北京市疫情的數(shù)據(jù)（ 數(shù)據(jù)來源：中華人民共和國衛(wèi)生部廣東非典型肺炎防治網(wǎng)中國疾病預(yù)防控制中心 ）日

14、期已確診病例累計醫(yī)務(wù)人員現(xiàn)有疑似病例死亡累計治愈出院累計4月20日3392440218334月21日4827861025434月22日5889966628464月23日69312678235554月24日77414386339644月25日87716095442734月26日988167109348764月27日1114187125556784月28日1199206127559 784月29日1347255135866834月30日1440270140875905月1日15532881415821005月2日16363001468911095月3日17413191493961155月4日180

15、332615371001185月5日189733515101031215月6日196034815231071345月7日204935815141101415月8日213637014861121525月9日217737214251141685月10日222737513971161755月11日226537914111201865月12日230438413781292085月13日234738413381342445月14日237038713081392525月15日238838913171402575月16日240539112651412735月17日242039312501453075月18日

16、243439412501473325月19日243739612491503495月20日244439512251543955月21日244439512211564475月22日245639512051585285月23日246539511791605825月24日249039611341636675月25日249939611051677045月26日250439710691687475月27日251239810051728285月28日25143989411758665月29日25173988031769285月30日252039876017710065月31日2521399747181108

17、76月1日252239973918111246月2日252239973418111576月3日252239972418111896月4日252239971818112636月5日252239971618113216月6日252239971318314036月7日252339966818314466月8日25223995501841543附錄2從附錄1的數(shù)據(jù)中，我們知道確診的醫(yī)務(wù)人員的人數(shù)從4月20日開始一直在增長，而且4月20日至5月8日增長速度很快，從5月9日至6月8日增長速度趨于緩和（如圖6.）。4月20日至5月8日之間每日新增的醫(yī)務(wù)人員病人數(shù)變化很大，這也反映出北京市SARS疫情這時期十

18、分復(fù)雜，有很多原本處在潛伏期內(nèi)的SARS感染者發(fā)作被確診，因而被隔離住院治療，同時很多醫(yī)務(wù)人員在與SARS病患者接觸時被感染。5月9日至6月8日每日新增的醫(yī)務(wù)人員病人數(shù)降為幾例，說明這個時期北京的SARS疫情已經(jīng)得較好的控制（如圖7.），被確診的SARS病人及住院的病人明顯降低，同時也反映出國家采取的干預(yù)政策起了很好的效果，人們的防范意識得到了增強(qiáng)。 圖6. 圖7.圖8.中我們可以看出從4月20日至5月16日非醫(yī)務(wù)人員的病人數(shù)一直在增長，在5月16日出現(xiàn)了普通人員感染人數(shù)的最大值，從5月17日至6月8日普通人員的病人數(shù)開始下降（如圖8.）。同時我們通過圖9.知道在5月16日以后新增的普通人員的

19、病人數(shù)為負(fù)值，這主要是由于普通人員沒有新的病人增加，相反每天有一些病人被治愈而導(dǎo)致。因此我們在模型中把北京市4月20日至6月8日的疫情分成兩個時期考慮，即4月20日至5月16日稱為感染期，5月17日至6月8日稱為恢復(fù)期。在恢復(fù)期內(nèi)，我們假設(shè)這時期內(nèi)普通人員病人的感染率為0。 圖8. 圖9.圖10.中每日新增的治愈病人數(shù)在5月15日之前變化很小，5月15日以后快速上升，圖11.中每日死亡的病人數(shù)5月15日之前變化很大，5月15日之后顯著下降，基本控制在1個人以內(nèi)。因此可以看出北京市的SARS疫情從5月15日之后得到了有效的控制。 圖10. 圖11.附錄3（問題處理的Matlab源程序）% 下列程

20、序計算醫(yī)務(wù)人員的日感染率k2， %普通病人的日感染率k1，死亡率d，治愈率z %clear allformat longN=10000000;%北京SARS確診，疑似，治愈，死亡人數(shù)的數(shù)據(jù)（4月20日-6月8日） %qz=339,482,588,693,774,877,988,1114,1199,1347,1440,1553,1636,1741,1803,1897,1960,2049,2136,2177,2227,2265,2304,2347,2370,2388,2405,2420,2434,2437,2444,2444,2456,2465,2490,2499,2504,2512,2514,2

21、517,2520,2521,2522,2522,2522,2522,2522,2522,2523,2522; dead=18,25,28,35,39,42,48,56,59,66,75,82,91,96,100,103,107,110,112,114,116,120,129,134,139,140,141,145,147,150,154,156,158,160,163,167,168,172,175,176,177,181,181,181,181,181,181,183,183,184; cure=33,43,46,55,64,73,76,78,78,83,90,100,109,115,118

22、,121,134,141,152,168,175,186,208,244,252,257,273,307,332,349,395,447,528,582,667,704,747,828,866,928,1006,1087,1124,1157,1189,1263,1321,1403,1446,1543;yw=24,78,99,126,143,160,167,187,206,255,270,288,300,319,326,335,348,358,370,372,375,379,384,384,387,389,391,393,394,396,395,395,395,395,396,396,397,3

23、98,398,398,398,399,399,399,399,399,399,399,399,399;pt=qz-yw-dead-cure;good=N-qz;infect=qz-dead-cure;%計算普通病人的日感染率disp('普通病人4月20日至5月16日，5月16日至6月8日的感染率分別為k11,k12')n=length(pt);k11=0;k12=0;j=0;for i=1:n-1 r(i)=pt(i+1)-pt(i); if(r(i)>0) k11=k11+r(i)/good(i+1); j=j+1; end end k11=k11/j k12%已確診的

24、普通病人數(shù)的數(shù)據(jù)圖形描繪%i=1:1:length(pt);%set(gcf,'position',300 200 560 420);%plot(pt,'ko')%hold on%plot(i,pt,'k-')%xlabel('時間(單位：天)');%ylabel('普通病人數(shù)');%title('普通病人數(shù)隨時間的變化情況');%新增的普通病人數(shù)的數(shù)據(jù)圖形描繪%i=1:1:n-1;%set(gcf,'position',300 200 560 420);%plot(r,'k

25、o')%hold on%plot(i,r,'k-')%xlabel('時間(單位：天)');%ylabel('每日新增的普通病人數(shù)');%title('新增的普通病人數(shù)隨時間的變化情況');%計算醫(yī)務(wù)人員的感染率disp('醫(yī)務(wù)人員病人4月20日至5月16日，5月16日至6月8日的感染率分別為k21,k22')n=length(yw);k21=0;k22=0;for i=1:n-1 r(i)=yw(i+1)-yw(i); if(i<=25) k21=k21+r(i)/infect(i+1); else

26、 k22=k22+r(i)/infect(i+1); endend k21=k21/25 k22=k22/(n-26)%已確診的醫(yī)務(wù)人員病人人數(shù)的數(shù)據(jù)圖形描繪%i=1:1:n;%set(gcf,'position',300 200 560 420);%plot(yw,'ko')%hold on%plot(i,yw,'k-')%xlabel('時間(單位：天)');%ylabel('已確診的醫(yī)務(wù)人員病人的人數(shù)');%title('已確診的醫(yī)務(wù)人員病人的人數(shù)隨時間的變化情況');%新增的醫(yī)務(wù)人員病人數(shù)隨

27、時間的變化情況圖形描繪%i=1:1:n-1;%set(gcf,'position',300 200 560 420);%plot(r,'ko')%hold on%plot(i,r,'k-')%xlabel('時間(單位：天)');%ylabel('每日新增的醫(yī)務(wù)人員病人數(shù)');%title('新增的醫(yī)務(wù)人員病人數(shù)隨時間的變化情況');%計算病人的日死亡率disp('病人4月20日至5月16日，5月16日至6月8日的死亡率分別為d1,d2')n=length(dead);d1=0;d2

28、=0;for i=1:n-1 r(i)=(dead(i+1)-dead(i); if(i<=25) d1=d1+r(i)/infect(i+1); else d2=d2+r(i)/infect(i+1); endend d1=d1/25 d2=d2/(n-26)%已確診的死亡病人數(shù)隨時間的變化情況圖形描繪%i=1:1:n;%set(gcf,'position',300 200 560 420);%plot(dead,'k-')%hold on%plot(i,dead,'ko')%xlabel('時間(單位：天)');%yla

29、bel('死亡病人數(shù)');%title('死亡病人數(shù)隨時間的變化情況');%新增死亡的病人數(shù)隨時間的變化情況圖形描繪%i=1:1:n-1;%set(gcf,'position',300 200 560 420);%plot(r,'k-')%hold on%plot(i,r,'ko')%xlabel('時間(單位：天)');%ylabel('每日新增的死亡病人數(shù)');%title('新增的死亡病人數(shù)隨時間的變化情況');%計算病人的日治愈率disp('病人4月2

30、0日至5月16日，5月16日至6月8日的治愈率分別為z1,z2')n=length(cure);z1=0;z2=0;for i=1:n-1 r(i)=(cure(i+1)-cure(i); if(i<=25) z1=z1+r(i)/infect(i+1); else z2=z2+r(i)/infect(i+1); endend z1=z1/25 z2=z2/(n-26)%治愈病人數(shù)隨時間的變化情況圖形描繪%i=1:1:n;%set(gcf,'position',300 200 560 420);%plot(cure,'k-')%hold on%pl

31、ot(i,cure,'ko')%xlabel('時間(單位：天)');%ylabel('死亡的病人數(shù)');%title('病人死亡人數(shù)隨時間的變化情況');%新增治愈的病人數(shù)隨時間的變化情況圖形描繪%i=1:1:n-1;%set(gcf,'position',300 200 560 420);%plot(r,'k-')%hold on%plot(i,r,'ko')%xlabel('時間(單位：天)');%ylabel('每日新增的治愈病人數(shù)');%ti

32、tle('新增的治愈病人數(shù)隨時間的變化情況');% % 下面的程序求模型的解i(t)的數(shù)值解 % %計算健康人數(shù)與病人數(shù)的比值n=length(qz);for i=1:n infect1(i)=infect(i)/N;%計算每日病人數(shù)占總?cè)藬?shù)的比值 t(i)=good(i)/infect(i);%計算健康人數(shù)與病人數(shù)的比值 end%對健康人數(shù)與病人數(shù)的比值做多項式擬合i=1:1:n;u=interp1(i,t,'spline'); %三次樣條插值%f=polyfit(i,t,3) %高次多項式擬合%u=polyval(f,i);%4月20日-6月8日模型解的統(tǒng)一表示for i=1:n if (i&l