基于SIR模型對埃博拉病毒感染人數(shù)增長率的確定

劉雅倩，劉元志，朱家明，章啟明

(1.安徽財經(jīng)大學(xué) 金融學(xué)院，安徽 蚌埠233030；2.安徽財經(jīng)大學(xué) 統(tǒng)計與應(yīng)用數(shù)學(xué)學(xué)院，安徽 蚌埠 233030；3.安徽財經(jīng)大學(xué) 經(jīng)濟學(xué)院；安徽 蚌埠 233030)

基于SIR模型對埃博拉病毒感染人數(shù)增長率的確定

劉雅倩1，劉元志3，朱家明2，章啟明3

(1.安徽財經(jīng)大學(xué) 金融學(xué)院，安徽 蚌埠233030；2.安徽財經(jīng)大學(xué) 統(tǒng)計與應(yīng)用數(shù)學(xué)學(xué)院，安徽 蚌埠 233030；3.安徽財經(jīng)大學(xué) 經(jīng)濟學(xué)院；安徽 蚌埠 233030)

針對埃博拉病毒感染人數(shù)增長率確定的問題，建立SIR改進模型，結(jié)合歷史數(shù)據(jù)，定量預(yù)測出未來感染人數(shù)增長率值，結(jié)合所建立的SIR改進模型，分析出當(dāng)接觸率為大于死亡率時，應(yīng)該研發(fā)新的藥品，藥品的生產(chǎn)計劃要滿足藥品的需求，通過分析不同階段的藥物需求，制定未來3年各生產(chǎn)地的生產(chǎn)計劃.

埃博拉；SIR傳染病模型；藥品生產(chǎn)計劃

2014年2月埃博拉[1]開始爆發(fā)于西非，截至2014年10月31日，世界衛(wèi)生組織關(guān)于埃博拉疫情報告稱，幾內(nèi)亞、利比里亞、塞拉利昂、美國、西班牙、馬里以及已結(jié)束疫情的尼日利亞與塞內(nèi)加爾累計出現(xiàn)埃博拉確診、疑似和可能感染病例13567例，其中4951人死亡.同時，埃博拉病毒的爆發(fā)也帶來了民航禁飛、人力損失、經(jīng)濟受創(chuàng)、糧食短缺等多方面的影響.如何有效地抗擊埃博拉病毒已經(jīng)成為全世界關(guān)注的熱點，然而，如果能夠掌握埃博拉病毒的發(fā)病特征并及時地研發(fā)出藥品，則會很好地控制埃博拉病情的蔓延.本文就是基于SIR傳染病模型，針對西非三國(幾內(nèi)亞、利比里亞與塞拉利昂)，思考如何預(yù)測埃博拉病毒感染人數(shù)的增長率與確定藥品生產(chǎn)計劃，從而迅速有效地控制埃博拉病毒蔓延(相關(guān)數(shù)據(jù)見2014年第三屆“認(rèn)證杯”數(shù)學(xué)中國數(shù)學(xué)建模國際賽試題C題[2]).

1 埃博拉病毒感染人數(shù)增長率的預(yù)測

1.1 研究思路

由于目前生產(chǎn)的藥品只能治愈早期病人，所以根據(jù)現(xiàn)實情況，我們對病人身體狀況進行進一步劃分，將病人分為早期病人與晚期病人，早期病人可以治愈，晚期病人只有面臨死亡.早期病人病愈之后成為治愈病人，查閱資料可知，當(dāng)埃博拉病人接受治療之后，對病毒存在免疫能力，不會再患病，所以患病病愈之后的健康者視為系統(tǒng)中的移出者.同時，那些晚期病人由于無法治愈成為死亡人，也將移除系統(tǒng).根據(jù)不同人群的身體狀況，我們將人群健康狀況劃分為以下5種：沒有免疫能力的健康人，有免疫能力的健康人，早期病人，晚期病人與死亡人，并將這5種狀態(tài)之間的轉(zhuǎn)移情況表示為如下圖形.

圖1 5種身體狀況狀態(tài)轉(zhuǎn)移圖

假定總?cè)藬?shù)N不變，將人群劃分為5類，即健康者、早期病人、晚期病人、死亡病人移出者與治愈病人移出者，時刻t5類人在總?cè)藬?shù)中的占比記為s(t)，e(t)，l(t)，d(t)和r(t)；

假設(shè)病人的接觸率為λ，死亡率為μ，β為早期病人到晚期病人的人數(shù)比率，α為早期病人治愈率(因為根據(jù)題意，我們假設(shè)只有早期病人可以治愈，晚期病人只有死亡).

可以得知s(t)+e(t)+l(t)+r(t)+c(t)=1，同時，對于死亡病人移出者而言應(yīng)有：

再記初始時刻的健康者和病人的比例分別為s0(s0>0)和i0(i0>0)，則可以表達出如下的方程：

1.2 參數(shù)值的確定

要對這個模型進行求解，必須能夠先得到這個模型中的未知參數(shù)的值[3].我們要知道的確定值得未知參數(shù)有：接觸率為λ，死亡率為μ，β為早期病人到晚期病人的人數(shù)比率，α為早期病人治愈率.其中μ，β和α都可以通過查閱資料的方式得到一個估計值.

λ值是最難確定的，而λ值是取決于一個地區(qū)的衛(wèi)生醫(yī)療水平.在當(dāng)我們設(shè)定一個地區(qū)或者國家時，其衛(wèi)生醫(yī)療水平在短時間內(nèi)不會變動，那么λ值也會相應(yīng)固定下來.這里我們使用最簡單的傳染病模型——SIR模型[4]對λ進行估值，由于使用數(shù)據(jù)一樣，都是相同國家或地區(qū)在相同時間的人口總數(shù)，病人累計數(shù)與死亡累計人數(shù)，所以通過SIR模型計算出來的λ值能相對近似擬合出模型中的λ值.

1.2.1λ值的確定——SIR模型

假定總?cè)藬?shù)N是定值，將人群劃分為三類，即健康者、病人、死亡病人的移出者，時刻t三類人在總?cè)藬?shù)中的占比記作s(t)，i(t)和r(t)，可以得到如下方程.

首先以幾內(nèi)亞對于建立的模型，設(shè)定μ=0.012524954，s0=0.999995625，i0=0.000004375.在這個模型中只有λ為參變量，不斷對λ值進行調(diào)試，使得預(yù)測模型的曲線變動情況盡可能地擬合出真實情況，得出最優(yōu)的擬合情況如下.

從圖2可以看出從250d開始，擬合曲線與真實值曲線出現(xiàn)明顯偏差，出現(xiàn)這樣偏差的原因可能是：在前期埃博拉病毒爆發(fā)時，社會各個部門還沒有時間做出有效措施，在后期，社會各部門的一些保護防范措施起到一定的作用，從而使埃博拉病毒的感染速率有所降低.所以這里以第250d為分割點，以第1～250d作為第一階段，250～311d作為第二階段，對圖形進行分段擬合，第一階段擬合情況較好，見圖3，對第二階段重新擬合，擬合結(jié)果如下圖.

圖2 病人比例的擬合結(jié)果 圖3 第二階段的患病人數(shù)的擬合結(jié)果

計算出250d之前的λ=0.028，250～311d，λ=0.018.由于λ值是由當(dāng)?shù)蒯t(yī)療衛(wèi)生水平?jīng)Q定的，并且為了真實擬合出現(xiàn)階段的λ，所以選取以離現(xiàn)階段時間最近的λ作為預(yù)測值，則使用λ=0.018.

同理，得到利比里亞與塞拉利昂兩個階段的λ值，并將計算出的λ值與查閱文獻獲得的值進行對比，可以看出計算結(jié)果相對比較合理正確.

表1 計算出λ值與查閱文獻獲得λ值的對比

1.3 參數(shù)值的確定

根據(jù)前面的分析，我們可以確定μ，β，基于這個方程，對早期病人與晚期病人的人數(shù)分別進行求解，發(fā)現(xiàn)當(dāng)早期病人治愈率α與接觸率λ不相等時，我們得到圖4和圖5，其中綠色的曲線表示的是晚期病人所占比率變動情況，藍色曲線表示的是早期病人所占比率的變動情況，可以看出隨著藥物的使用，晚期病人人數(shù)先快速上升最后慢慢下降直至趨于零，早期病人人數(shù)會一直下降直至趨于零；當(dāng)早期病人治愈率α與接觸率λ相等時，我們得到圖5，可以看出隨著藥物的使用，晚期病人人數(shù)先快速上升最后維持在一個水平，早期病人人數(shù)會一直下降直至維持于一個水平.

圖4 第一種晚期病人與早期病人人數(shù)變動情況 圖5 第二種晚期病人與早期病人人數(shù)變動情況

2 藥品生產(chǎn)計劃的確定

2.1 是否研發(fā)藥品

首先，疾病本身的特性是決定是否需要生產(chǎn)新藥品的關(guān)鍵.如果一種疾病是具有顯著周期性或者說，這個疾病不會根本消失，那么就有必要準(zhǔn)備生產(chǎn)藥品以防備疾病隨時來襲.

其次，根據(jù)研究一所建立的模型可以分析出σ=1是一個閾值，因為當(dāng)σ>1時，i(t)的增減性取決于i0的大小，但是其極限值隨著σ的增加而增加，當(dāng)σ≤1時，病人比例i(t)越來越小，最終趨近于0，這是由于傳染期內(nèi)經(jīng)有效解除從而使健康者變成的病人數(shù)不超過原來病人數(shù)的緣故.所以當(dāng)σ>1時，即為當(dāng)3d累計接觸率為λ大于3d累計死亡率ω，應(yīng)該研發(fā)新的藥品.

2.2 藥品生產(chǎn)計劃的確定

根據(jù)一般藥物的規(guī)模化生產(chǎn)進程，假設(shè)該藥物在半年后可以達到較大的生產(chǎn)能力，可以滿足同時期的藥物需求量，其生產(chǎn)能力的增長是線性的.

結(jié)合一般疫苗的保質(zhì)期，同時考慮節(jié)約成本，藥物3個月進行一次運輸，但其中第一個月剛開始生產(chǎn)藥品時應(yīng)及時運輸.即第一個月為周期一，此后周期長度都是3個月.

因為該藥物只能治愈早期病人，因此以早期病人數(shù)量為研究對象，作為藥物需求的依據(jù)，從而制定生產(chǎn)計劃.為平衡生產(chǎn)壓力，每個產(chǎn)地生產(chǎn)相同數(shù)量的藥物.

2.2.1 模型的建立

藥品的生產(chǎn)計劃要滿足藥品的需求，通過分析不同階段的藥物需求，制定未來3年5個生產(chǎn)國的生產(chǎn)計劃[5].

⑴第一年

生產(chǎn)能力需要經(jīng)歷逐步增長的過程，在初期沒有大規(guī)模生產(chǎn)的能力，因此，第一年投放藥物的目的是控制埃博拉疫情，使其不爆發(fā).第一年的第一個月，即生產(chǎn)計劃的第1期，新藥物初步進行生產(chǎn)，沒有藥物投放到災(zāi)區(qū)，災(zāi)區(qū)疫情按照目前形勢發(fā)展，第1期生產(chǎn)的藥物用來滿足災(zāi)區(qū)第2期的需求，則有：

其中，i為生產(chǎn)計劃的期數(shù)，Di為第i期的藥物需求量，k表示生產(chǎn)初期藥物生產(chǎn)力隨時間增長的能力，t為時間(d)，ei為第i期的早期病人比例，P為非洲3國總?cè)丝冢?是治愈所有患病人的治愈率，α2是使埃博拉疫情不至于爆發(fā)時的治愈率.

同理，對于生產(chǎn)計劃的第二期：

其中x表示，在第二期結(jié)束前x天，可以停止生產(chǎn)擴張，此時的生產(chǎn)能力已經(jīng)可以滿足需求.

半年后，從第三期開始，生產(chǎn)能力已經(jīng)足夠滿足需求，此時藥物的供需平衡滿足：

⑵第二年

生產(chǎn)力經(jīng)過第一年發(fā)展，已經(jīng)可以滿足非洲三國的需求，第二年的目的是提供足夠的藥物，治療所有患病者，此時：

⑶第三年

第三年的目的是在治療所有患病者，為健康人提供疫苗，徹底根除埃博拉病毒.提供疫苗的數(shù)量與當(dāng)時的患病人數(shù)成正比，我們合理假設(shè)所有與病人接觸過的人都要接種疫苗，則此比例應(yīng)該等于本文改進SIR模型中的接觸率，此時：

2.2.2 模型的解決

首先，根據(jù)研究一所建立的SIR改進模型，通過MATLAB編程求解，得到在此藥物供給計劃下，非洲3國未來3年12期的早期病人數(shù)量之和.

如圖6，第一年前期期埃博拉仍在以較快速度傳播，是因為剛研制出藥物，還沒有投入使用；第一年后期埃博拉病毒得到抑制，是因為藥物投入使用，使早期病人數(shù)量逐漸減少并趨于穩(wěn)定，第二年、第三年均能向災(zāi)區(qū)供給充足的藥物，因此早期病人數(shù)量逐漸減少直至趨近于0.

根據(jù)前面建立的公式，結(jié)合未來3年早期病人預(yù)測值，生產(chǎn)力發(fā)展時線性增長；生產(chǎn)力充足，每天生產(chǎn)的藥物數(shù)量剛好滿足每天需求.此時，可以用生產(chǎn)地每天的產(chǎn)量表示每日實際發(fā)揮的生產(chǎn)力的變化.

圖6 藥物供給計劃下的早期病人預(yù)測 圖7 每天運輸藥物時實際發(fā)揮的生產(chǎn)力(所有產(chǎn)地生產(chǎn)量之和)

由于實際生活中不是每天都要運輸藥物，而是按周期運輸，因此在每一期的生產(chǎn)中采用平均速度生產(chǎn)即可，可以得到每個生產(chǎn)地的生產(chǎn)計劃.

表2 每個生產(chǎn)地每期面臨的總需求和生產(chǎn)速度(劑量/d)

3 總 結(jié)

以西非3國(幾內(nèi)亞、利比里亞與塞拉利昂)作為分析對象，建立SIR改進模型，對埃博拉病毒感染人數(shù)增長率進行預(yù)測.基于建立的SIR改進模型，分析出當(dāng)接觸率大于死亡率，應(yīng)該研發(fā)新的藥品，并預(yù)測出每一期的生產(chǎn)中每個生產(chǎn)地的生產(chǎn)計劃.本文將SIR模型與現(xiàn)實意義結(jié)合，對進一步對抗埃博拉病毒蔓延有一定的參考價值.

[1] 劉安立.埃博拉全解讀[J].大自然探索,2014(10):8-21.

[2] 第三屆數(shù)學(xué)中國數(shù)學(xué)建模國際賽：http://mcm.tzmcm.cn/[EB/OL].2014-12-01

[3] 何艷輝,唐三一.經(jīng)典SIR模型辨識和參數(shù)估計問題[J].應(yīng)用數(shù)學(xué)和力學(xué),2013(03):2-4

[4] 劉冠軍.SIR傳染病模型的辨識分析及穩(wěn)定性分析[D].信陽:信陽師范學(xué)院碩士學(xué)位論文,2011.

[5] 王漢斌,于軍生.基于關(guān)鍵鏈理論的多項目生產(chǎn)計劃制訂與控制模型構(gòu)建[J].中國高新技術(shù)企業(yè),2010(06):55-57.

[責(zé)任編輯：王軍]

Determination on the infections growth rate of the Ebola virus based on the SIR model

LIU Yaqian1, LIU Yuanzhi3,ZHU Jiaming2, ZHANG Qiming3,

(1.School of Finance Anhui University of Finance and Economics;2.School of Statistics and Appl. Math;3. School of Economic,Bengbu 233030，China)

As for determining the infections growth rate of the Ebola virus, we establish SIR improvement model. We predict the growth rate value of infecting since the outbreak of Ebola after our quantitative analysis. Using the established SIR improved model, we come to the conclusion that when the contact rate is greater than the death rate, it is time to develop a new drug. At the same time, after analyzing the drugs demand in different stage, we predict the deadline for research and development and the production of drugs.

Ebola; SIR epidemic model; drug research and development

2014-11-06

國家自然科學(xué)基金資助項目(11301001)；安徽財經(jīng)大學(xué)教研項目(acjyzd201429)

朱家明(1973-),男,安徽泗縣人,安徽財經(jīng)大學(xué)副教授,碩士,安徽財經(jīng)大學(xué)數(shù)學(xué)建模實驗室主任,主要從事應(yīng)用數(shù)學(xué)與數(shù)學(xué)建模的研究.

R373.32

A

1672-3600(2015)06-0025-05