孫 海， 熊思燦， 吳志強(qiáng)

(1.樂(lè)山師范學(xué)院數(shù)學(xué)與信息科學(xué)學(xué)院，四川樂(lè)山 614000;2.東華理工大學(xué)理學(xué)院，江西撫州 344000)

甲型H1N1流感是一種由A型甲流感病毒引起的豬呼吸系統(tǒng)疾病，該病毒可在豬群中造成流感的爆發(fā)。目前，此種病毒已在人群中大量爆發(fā)，它的潛伏期較流感、禽流感潛伏期長(zhǎng)，病毒可能在人體潛伏7天后才表現(xiàn)出病癥，感染后一般會(huì)在一周或一周多后就發(fā)病。所以分析感染人群的變化規(guī)律具有十分重要的意義。

霍闊等(2010)采用四階Runge-kutta方法，利用一般SIR模型(姜啟源等，2003;劉來(lái)福等，1997;李偉等，2004;李建奎等，2004;宇永仁等，2005;王汝發(fā)等，2004;周永衛(wèi)等，2007)對(duì)甲型 H1N1流感傳播進(jìn)行了研究，本文在此模型的基礎(chǔ)上進(jìn)行了改進(jìn)，增加了人群種類(lèi):潛伏者和接種疫苗者，且對(duì)他們間的轉(zhuǎn)化系數(shù)也做了相應(yīng)的改變，增加了每日潛伏者中轉(zhuǎn)化為病毒攜帶者的比例和易感染者的接種疫苗比例，這樣就使得模型更接近實(shí)際情況。

1 改進(jìn)SIR模型的建立和求解

1.1 模型的假設(shè)與符號(hào)的約定

1.1.1 模型假設(shè)

(1)假設(shè)人群總數(shù)不變，并且不考慮這段時(shí)間內(nèi)的人口出生率和自然死亡率及因其他疾病死亡和國(guó)內(nèi)人口流動(dòng)情況;

(2)假設(shè)以天作為時(shí)間計(jì)量單位;

(3)假設(shè)甲型H1N1流感潛伏期為7 d，在發(fā)病前1 d就具有傳染性;

(4)假設(shè)由于一旦發(fā)現(xiàn)流感，就會(huì)加強(qiáng)防控，而且防控措施到位，凡感染者均能去醫(yī)院就醫(yī)，并成為確診者;

(5)假設(shè)治愈恢復(fù)者對(duì)甲型H1N1病毒有免疫力，不再感染;

(6)假設(shè)甲型H1N1流感的死亡率為0;

(7)假設(shè)所有的統(tǒng)計(jì)數(shù)據(jù)真實(shí)，沒(méi)有遺漏現(xiàn)象。

1.1.2 符號(hào)約定

S(t)為t時(shí)刻易感染者在人群中所占的比例;I(t)為t時(shí)刻感染者在人群中所占的比例;R(t)為截止到t時(shí)刻累計(jì)治愈恢復(fù)者在人群中所占的比例;V(t)為截止到t時(shí)刻累計(jì)接種疫苗者在人群中所占的比例;L(t)為潛伏者在人群中所占的比例;r為每天已確認(rèn)成為治愈者的比例;p為每個(gè)攜帶甲流病毒者的日平均傳染率;q為每日潛伏者中轉(zhuǎn)化為病毒攜帶者的比例;l為每日潛伏者發(fā)病成為感染者的比例;v為易感染者的接種疫苗比例。

1.2 模型的建立

由以上假設(shè)和符號(hào)的約定，可以得到改進(jìn)的SIR模型甲型H1N1流感各類(lèi)人群間的關(guān)系轉(zhuǎn)化如圖1所示:

圖1 改進(jìn)的SIR模型的各類(lèi)人群間的關(guān)系圖Fig.1 Relationship graph of all kinds of people based on the improved SIR model

由圖1可以看出，不同人群間增加了更多的轉(zhuǎn)換過(guò)程:易感染者可以向潛伏者或接種疫苗者轉(zhuǎn)化;潛伏者可以向感染者轉(zhuǎn)化;感染者可以向治愈者轉(zhuǎn)化，使得模型與實(shí)際情況更加符合，并由此我們可以得到改進(jìn)的SIR模型的甲型H1N1流感防控模型的常微分方程如式1所示:

1.3 模型的求解

設(shè)日平均傳染率p為1.5，易感染者的接種疫苗比例v為0.002。由于假設(shè)甲型H1N1流感潛伏期為7 d，故每日潛伏者成為傳染者的比例l為1/7，每日潛伏者中轉(zhuǎn)化為病毒攜帶者的比例q為1/6，本文中 r取估計(jì)值 0.309 4。

對(duì)改進(jìn)的SIR模型，利用一般的解微分方程的方法難以求解，本文采用德國(guó)學(xué)者(Felhberg，1968)提出的四階五級(jí)RKF算法。假設(shè)當(dāng)前的步長(zhǎng)為hk，則可以定義下面的6個(gè)ki變量:

式中，tk為當(dāng)前計(jì)算時(shí)刻，中間參數(shù)αi，βij及其他參數(shù)由表1(四階五級(jí)RKF算法系數(shù)表)給出，其中αi，βij參數(shù)對(duì)又稱為 Dormand-Prince對(duì)。這時(shí)下一步的狀態(tài)向量可以由下式求出:

表1 四階五級(jí)RKF算法系數(shù)表Table 1 The coefficient of four order and five stage RKF method

利用上述方法對(duì)模型進(jìn)行求解，得到改進(jìn)SIR模型 S(t)，I(t)，R(t)，V(t)，L(t)隨時(shí)間變化的曲線圖(圖2)。

由圖2可以看出，S(t)曲線隨時(shí)間遞減，表明易感染者受到傳染或接種疫苗而成為潛伏者或接種疫苗者，所以數(shù)量不斷減少;I(t)曲線隨時(shí)間呈先遞增后遞減的趨勢(shì)，這是由于易感染者轉(zhuǎn)化為感染者，所以感染者數(shù)量逐漸增大，增大到一定值后，隨著中后期防控措施的加強(qiáng)，感染者最終將全部轉(zhuǎn)化為治愈者，直至為0;L(t)曲線隨時(shí)間呈先遞增后遞減的趨勢(shì)，這是由于潛伏者轉(zhuǎn)化為確診者，所以確診者數(shù)量逐漸增大，增大到一定值后，隨著中后期防控措施的加強(qiáng)，確診者最終將全部轉(zhuǎn)化為治愈恢復(fù)者，直至為0;R(t)曲線隨時(shí)間呈遞增趨勢(shì)，這是由于確診者不斷治愈成為治愈恢復(fù)者，到最后趨于一個(gè)穩(wěn)定值。V(t)曲線隨時(shí)間呈遞增趨勢(shì)，這是由于易感染者不斷接種疫苗而成為接種疫苗者。

1.4 參數(shù)修正

對(duì)日平均傳染率對(duì)參數(shù)p值進(jìn)行修正，再對(duì)SIR模型和改進(jìn)的SIR模型進(jìn)行求解。當(dāng)p為1.4，1.3，1時(shí)，SIR模型和改進(jìn)的 SIR 模型的 S(t)，I(t)，R(t)，V(t)，L(t)隨時(shí)間變化的曲線圖如圖3～5所示。

由圖3可以看出，隨著p值的不斷減小，SIR模型的S(t)曲線，I(t)曲線，R(t)曲線隨時(shí)間的變化趨勢(shì)不大，并且最后易感染者和感染者最后幾乎都轉(zhuǎn)化成了治愈恢復(fù)者，這顯然與實(shí)際情況不符合，改進(jìn)的SIR模型的S(t)曲線，I(t)曲線，R(t)曲線，V(t)曲線，L(t)曲線隨時(shí)間呈遞增或遞減趨勢(shì)越來(lái)越明顯，這是因?yàn)?隨著接觸率的不斷減少，以及后期防控措施的加強(qiáng)，易感染者數(shù)量在不斷減少，接種疫苗者一直隨時(shí)間呈遞增趨勢(shì)，這是由于易感染者不斷接種疫苗，確診者和潛伏者已與模型無(wú)關(guān)，治愈恢復(fù)者變?yōu)?。這說(shuō)明所建的模型與實(shí)際情況越來(lái)越接近，模型越來(lái)越合理。

圖2 防控模型的預(yù)測(cè)曲線圖Fig.2 The prediction of prevention model

圖3 SIR曲線圖及改進(jìn)SIR曲線圖Fig.3 The SIR graphics and improved SIR graphics

2 基于改進(jìn)SIR模型的甲型H1N1流感防控的檢驗(yàn)

2.1 I(t)預(yù)測(cè)值與實(shí)際值的對(duì)比

當(dāng) p 為1.5，1.4，1.3，1 時(shí)，得到 SIR 模型的預(yù)測(cè)I(t)、改進(jìn)的SIR模型的預(yù)測(cè)I(t)和實(shí)際I(t)隨時(shí)間變化的曲線圖(圖4)。

由圖4可以看出，當(dāng)參數(shù)p逐步減小時(shí)，改進(jìn)的SIR模型的預(yù)測(cè)I(t)值和實(shí)際I(t)值越來(lái)越吻合。當(dāng)p值為1時(shí)，SIR模型的預(yù)測(cè)I(t)值卻與實(shí)際預(yù)測(cè)I(t)值存在較大的偏差，而改進(jìn)的SIR模型的預(yù)測(cè)I(t)值和實(shí)際I(t)值幾乎吻合，這進(jìn)一步說(shuō)明了改進(jìn)的SIR模型的合理性。

2.2 R(t)預(yù)測(cè)值與實(shí)際值的對(duì)比

當(dāng) p 為 1.5，1.4，1.3，1 時(shí)，得到 SIR 模型的預(yù)測(cè)R(t)、改進(jìn)的SIR模型的預(yù)測(cè)R(t)和實(shí)際R(t)隨時(shí)間變化的曲線圖(圖5)。

圖4 日平均感染率取不同值時(shí)感染者比例的預(yù)測(cè)值、改進(jìn)后預(yù)測(cè)值和實(shí)際值的對(duì)比圖Fig.4 The comparison of the infected rate of the prediction，the improved prediction and the actual with different daily average infection rate

由圖5可以看出，當(dāng)參數(shù)p逐步減小時(shí)，改進(jìn)的SIR模型的預(yù)測(cè)R(t)值和實(shí)際R(t)值越來(lái)越吻合。當(dāng)p值為1時(shí)，SIR模型的預(yù)測(cè)R(t)值卻與實(shí)際預(yù)測(cè)R(t)值存在較大的偏差，而改進(jìn)的SIR模型的預(yù)測(cè)R(t)值和實(shí)際的R(t)值幾乎完全吻合，這也進(jìn)一步說(shuō)明了改進(jìn)的SIR模型的合理性。

3 結(jié)論

較之一般的SIR模型，本文的人群分類(lèi)更準(zhǔn)確，因素考慮更全面。建立模型時(shí)，增加了對(duì)潛伏者、接種疫苗者的考慮，并將甲型H1N1流感患者在發(fā)病前一天即具有傳染性這一因素也考慮了進(jìn)去，使得所建模型更合理，所得結(jié)果與實(shí)際情況更一致。但由于沒(méi)有隱性傳染者相關(guān)數(shù)據(jù)的支持，所以無(wú)法獲得隱性傳染者與其他人群間的轉(zhuǎn)換關(guān)系，故未將其列入所建數(shù)學(xué)模型中，這也成了模型的一大缺點(diǎn)。

圖5 日平均感染率取不同值時(shí)治愈者比例的預(yù)測(cè)值、改進(jìn)后預(yù)測(cè)值和實(shí)際值的對(duì)比圖Fig.5 The comparison of the cured rate of the prediction，the improved prediction and the actual with different daily average infection rate

霍闊，李世霖.2010.甲型H1N1流感傳播的SIR模型研究［J］.湖南工業(yè)大學(xué)學(xué)報(bào)，24(4):40-42.

姜啟源，謝金星，葉俊.2003.數(shù)學(xué)模型［M］.北京:高等教育出版社.

劉來(lái)福，曾文藝.1997.數(shù)學(xué)模型與數(shù)學(xué)建模［M］.北京:北京師范大學(xué)出版社.

李偉.2004.關(guān)于SARS病毒傳播的數(shù)學(xué)模型［J］.畢節(jié)師范高等?？茖W(xué)校學(xué)報(bào)，22(2):46-52.

李建奎，劉天喜.2004.SARS病毒傳染的數(shù)學(xué)模型［J］.科技情報(bào)開(kāi)發(fā)與經(jīng)濟(jì)，4(2):163-164.

宇永仁.2005.SARS傳播模型及其對(duì)經(jīng)濟(jì)的影響［J］.遼寧大學(xué)學(xué)報(bào)，32(1):1-2.

王汝發(fā).2004.SARS傳播的數(shù)學(xué)模型分析［J］.數(shù)理醫(yī)藥學(xué)雜志，17(2):99-100.

周永衛(wèi)，范賀花.2007.傳染病數(shù)學(xué)模型的馬爾可夫骨架過(guò)程建模［J］.安陽(yáng)師范學(xué)院學(xué)報(bào)，(2):33-35.

Fehlberg E.1968.Classical fifth-，sixth-，seventh-，and eighth-order Runge-Kutta formulas with stepsize control［R］.NASA Technical Report，287:87.