基于SEIR 的傳染病傳播模型

朱 柳 鄭婉婷 張貴層

(北方工業(yè)大學(xué),北京100000)

縱觀歷史,一旦出現(xiàn)較大型的傳染病傳播情況,對(duì)世界的經(jīng)濟(jì)發(fā)展和人民生活都會(huì)帶來了很大影響。因此想要預(yù)測(cè)和控制傳染病的蔓延,我們就必須了解一般傳染病毒的傳播機(jī)制,定量地研究傳染病的傳播規(guī)律,找出其傳播的數(shù)學(xué)模型。

本文結(jié)合實(shí)際影響因素建立了一個(gè)最接近現(xiàn)實(shí)情況的傳播模型,并通過建立的預(yù)測(cè)模型對(duì)病毒傳播未來一段時(shí)間發(fā)展態(tài)勢(shì)進(jìn)行預(yù)測(cè)。

1 傳播模型的建立

1.1 傳染病傳播機(jī)制

針對(duì)一般傳染病的傳播情況,把人群大致分為四類:(1)健康者S(即沒有被感染的人群,有被感染的幾率)；(2)潛伏者E(已經(jīng)被感染的人群,但是沒有出現(xiàn)癥狀。根據(jù)查閱相關(guān)資料[1]潛伏者也具有一定的傳染性)；(3)感染者I(即已經(jīng)被感染并且出現(xiàn)一定的癥狀)；(4)排除者R(即康復(fù)人群和死亡人群,不再參與到感染的過程)。通過相關(guān)資料[2],顯示潛伏期大約為1—14天,本文對(duì)感染者和潛伏者的定義為,兩天以下的潛伏期人群稱為感染者,三天及以上的潛伏期人群稱為潛伏者。

設(shè)t 時(shí)刻SEIR 四類人口數(shù)分別記為s(t)、e(t)、i(t)和r(t),某國(guó)的人口總數(shù)為N。

α1為感染者I 的單位時(shí)間傳播率,即單位時(shí)間內(nèi)一個(gè)感染者會(huì)導(dǎo)致α1S(t)/N 個(gè)健康人帶上病毒。其中帶病毒的人群有k1的比例抵抗力較差從而會(huì)直接患病成為感染者I,有k2比例因?yàn)樽陨碛休^強(qiáng)的抵抗力,沒有立即患病,但是在之后也有很大患病幾率,成為潛伏者E。

α2為潛伏者E 的單位時(shí)間傳播率,即單位時(shí)間內(nèi)一個(gè)潛伏者會(huì)導(dǎo)致α2S(t)/N 個(gè)健康人帶上病毒,成為潛伏者E,因?yàn)闈摲叩膫魅韭时雀腥菊叩膫魅韭矢?大部分潛伏者傳染別人使別人變成潛伏者E,而不會(huì)使得其直接變成感染者I。

γ1為感染者被嚴(yán)格隔離的比例,γ2為潛伏者中被嚴(yán)格隔離的比例。隨著政府干涉以及人民的防范意識(shí)增強(qiáng),隔離率隨時(shí)間成線性遞增關(guān)系。

β1為單位時(shí)間潛伏者E 變?yōu)楦腥菊逫 的比例,β2為單位時(shí)間潛伏者E 變?yōu)榻】嫡逽 的比例,變成健康者S 后仍然有被感染的風(fēng)險(xiǎn)。

l1為單位時(shí)間內(nèi)感染者I 的治愈率,l2為單位時(shí)間內(nèi)感染者的死亡率。為直觀看出其傳播機(jī)制,本文做出了病毒的傳播圖。

圖1 傳染病傳播圖

1.2 病毒傳播模型微分方程

根據(jù)傳播機(jī)制建立關(guān)于s(t),e(t),i(t)和r(t)的方程:

1.2.1 四種人群人數(shù)總和為N,保持不變:

1.2.2 等式左邊表示單位時(shí)間內(nèi)健康者S 的變化,等式右邊從左向右第一部分是正,表示單位時(shí)間內(nèi)潛伏者會(huì)有β2比例變成健康者。第二部分是負(fù),表達(dá)單位時(shí)間未被隔離的感染者(比例為1-γ1)會(huì)感染健康者使得健康者變?yōu)楦腥菊呋驖摲?傳播率為α1。第三部分是負(fù),表達(dá)單位時(shí)間未隔離的潛伏者(比例為1-γ2)會(huì)感染健康者使得健康者變?yōu)闈摲?傳播率為α2。

1.2.3 等式左邊表示單位時(shí)間內(nèi)潛伏者E 的變化,等式右邊從左向右第一部分是正,表示單位時(shí)間未被隔離的感染者(比例為1-γ1)會(huì)感染健康者使得健康者變?yōu)闈摲?比例為k2,傳播率為α1。第二部分是正,表達(dá)單位時(shí)間未被隔離的潛伏者(比例為1-γ2)會(huì)感染健康者使得健康者變?yōu)闈摲?傳播率為α2。第三部分為負(fù),表示單位時(shí)間內(nèi)潛伏者會(huì)有β1比例變成感染者。第四部分為負(fù),表示單位時(shí)間內(nèi)潛伏者會(huì)有β2比例變成健康者。

1.2.4 等式左邊表示單位時(shí)間內(nèi)感染者I 的變化,等式右邊從左向右第一部分是正,表示單位時(shí)間未隔離的感染者(比例為1-γ1)會(huì)感染健康者使得健康者變?yōu)楦腥菊?比例為k1,感染率為α1。第二部分是正,表達(dá)單位時(shí)間內(nèi)潛伏者會(huì)有β1比例變成感染者。第三部分為負(fù),表示單位時(shí)間內(nèi)有l(wèi)1比例的感染者會(huì)被治愈成為排除者R,第四部分為負(fù),表示單位時(shí)間內(nèi)有l(wèi)2比例的感染者會(huì)死亡成為排除者R。

1.2.5 等式左邊表單位時(shí)間內(nèi)排除者R 的變化,等式右邊第一部分為正,表示單位時(shí)間內(nèi)有l(wèi)1比例的感染者會(huì)被治愈成為排除者R,第二部分為負(fù),表示單位時(shí)間內(nèi)有l(wèi)2比例的感染者會(huì)死亡成為排除者R

1.2.6 單位時(shí)間內(nèi)感染者會(huì)感染健康人群使得比例為k1的健康者變成感染者,比例為k2的健康者變成潛伏者,它們之和為1。

聯(lián)立(1)-(6)式進(jìn)一步整理可得以下方程:

設(shè)健康者S、潛伏人群E、感染者I 和康復(fù)人群R 初始值分別為s1,e1,i1,r1。

1.3 系數(shù)優(yōu)化模型的建立

通過現(xiàn)有的數(shù)據(jù),可以擬合出感染者I 以及排除者R 隨時(shí)間的變化曲線。在模型中有兩個(gè)重要的參數(shù):單位時(shí)間傳播率α1,以及單位時(shí)間內(nèi)潛伏者向感染者轉(zhuǎn)化的比率β1。通過現(xiàn)有資料,大致給出一個(gè)范圍,即為優(yōu)化模型的約束條件。通過建立的系數(shù)優(yōu)化模型,可以用智能算法擬合出最符合真實(shí)數(shù)據(jù)的曲線,以得到最佳系數(shù)。

2 算例分析

2.1 傳播模型系數(shù)的估計(jì)

人口總N:經(jīng)過網(wǎng)站[3]查詢某國(guó)2020 人口總數(shù)為143,964,709。

2.1.1 確定系數(shù)α1α2

α1:為未被隔離的感染者的傳播率,是一個(gè)重要參數(shù),對(duì)其進(jìn)行優(yōu)化得出其值,α2為未被隔離的潛伏者的傳播率,由相關(guān)文獻(xiàn)[4]稱,4 成左右的潛伏者有感染能力,在此設(shè)α2=0.4×α1,α1范圍為:[0.01-1]。

2.1.2 確定系數(shù)γ1γ2

γ1γ2分別為感染者的平均隔離率和潛伏者的平均隔離率,沒有確定統(tǒng)計(jì)的值,隔離率隨著時(shí)間線性遞增,假設(shè):隔離率,γ2=γ1=0.2+0.01×t。

2.1.3 確定系數(shù)β1β2

β1為單位時(shí)間潛伏者變成感染者的比例,β2為單位時(shí)間潛伏者變?yōu)榻】嫡叩谋壤刹殚喯嚓P(guān)資料可得,最終潛伏者中有45%的比例變?yōu)楦腥菊?那么55%的比例會(huì)變成健康者,基于目前的流行病學(xué)調(diào)查,潛伏期1-14 天,多為3-7 天,范圍為:[0.01-0.45]。

2.1.4 確定系數(shù)l1l2

根據(jù)真實(shí)數(shù)據(jù)用Matlab2018a 可繪出:

圖2 單位時(shí)間死亡率和治愈率曲線

記3.29 為第一天,深色線條代表單位時(shí)間治愈率,淺色線條代表單位時(shí)間死亡率,單位時(shí)間死亡率基本保持一樣,取其平均值代表其值,算得l2單位時(shí)間死亡率為0.001。深色線有一定波動(dòng),用cftool 工具進(jìn)行擬合。得到最擬合方程式。

圖3 單位時(shí)間治愈率擬合曲線圖

l1單位時(shí)間治愈率=f(t)=p1×t3+p2×t2+p3×t+p4 p1=2.447e-06 P2=-0.000115 P3=0.001307 P4=0.009866

2.1.5 確定系數(shù)k1k2

k1為被感染者感染的健康者變?yōu)楦腥菊叩谋壤?k2為被感染者感染的健康者變?yōu)闈摲叩母怕省摲跒?～14 天,在此定義潛伏期2 天以下的就視為感染者,三天及以上視為潛伏者。

根據(jù)現(xiàn)有文獻(xiàn)研究[4],可以算出潛伏期兩天以下的占比約為45%,即k1為45%,潛伏期三天及以上的占比約為k2為55%。

2.2 最佳系數(shù)優(yōu)化解

用MATLAB2018b 編程求解 由粒子群算[5]法擬合出的數(shù)據(jù)為

α1=0.2923 β1=0.40

擬合值和真實(shí)值效果圖如圖4:

圖4 擬合值真實(shí)值比較圖

2.3 R0 的最終確定

基本傳染數(shù)(basic reproduction number,R0)是一個(gè)流行病學(xué)術(shù)語,在沒有外力介入、所有人都沒有免疫力的情況下,一個(gè)感染者在具有傳染性的這段時(shí)間內(nèi),平均可以傳染給多少人。

在本文的傳播模型中,把隔離率設(shè)置為0 來模擬沒有外力介入的情況。將用粒子群算法擬合出來的單位時(shí)間傳播率α1=0.2923,和潛伏者轉(zhuǎn)化率β1=0.40 帶入到傳播模型。

圖5 有干涉和無干涉病毒傳染比較圖

根據(jù)相關(guān)文獻(xiàn)[6]可得:

R0=(1+rTL)(1+rTI)

其中r 表示指數(shù)增長(zhǎng)的增長(zhǎng)率,若用b(t)表示第t 天的新增被感染數(shù),則:

b(t)=b(t-△t)er△t

TL表示感染者潛伏期的平均長(zhǎng)度,SI 表示一個(gè)感染者被感染的時(shí)間和他感染的下一個(gè)人被感染的時(shí)間的間隔,TI=SI-TL,這三個(gè)數(shù)據(jù)都來自對(duì)現(xiàn)實(shí)病人情況的觀察,通過查詢文獻(xiàn)[7]可以估計(jì)SI 為8.4 天,TL為6 天。

最后剩一個(gè)未知數(shù)據(jù)r 為指數(shù)增長(zhǎng)率,本文用擬合的無干預(yù)(即隔離率為零)數(shù)據(jù)來擬合r。

所以r=Ln1.1787=0.1644

最終求得R0=(1+rTL)(1+rTI)=(1+0.1644×6)(1+0.1644×2.4)=2.77

3 病毒傳播數(shù)據(jù)的預(yù)測(cè)

傳播模型對(duì)病毒傳播數(shù)據(jù)的預(yù)測(cè):

用已做好的模型進(jìn)行預(yù)測(cè),記3.29 日為第一天,前面的傳播模型有現(xiàn)存感染者I,和R 的數(shù)據(jù),用微分方程很容易分別求得累計(jì)治愈人數(shù)以及累計(jì)死亡人數(shù)、新增確診人數(shù)為:

當(dāng)日存在感染者+當(dāng)日累計(jì)治愈+當(dāng)日累計(jì)死亡- (前一日日存在感染者+前一日日累計(jì)治愈+前一日累計(jì)死亡),用MATLAB 擬合圖線(3.29 日為第一天)可得:

圖6 病毒傳播情況預(yù)測(cè)圖

列出相應(yīng)表格可得出:

表1 病毒傳播情況預(yù)測(cè)表

同時(shí)可以看出在5.14-5.16,將會(huì)迎來病毒傳播的拐點(diǎn),感染者人數(shù)會(huì)下降,新增確診感染者數(shù)也會(huì)下降。

4 模型結(jié)果分析

通過現(xiàn)有的文獻(xiàn)可得一般的R0為2-3,解得的R0大致為2.77,所以傳播模型的建立應(yīng)該具有合理性。通過4.30 日之前的數(shù)據(jù),預(yù)測(cè)出5.1 到5.20 號(hào)的數(shù)據(jù)。根據(jù)現(xiàn)有的5.1,5.2 日的真實(shí)數(shù)據(jù)。比較其誤差,累出誤差比較表可得。(表2)

表2 誤差比較表

根據(jù)5.1、5.2 日預(yù)測(cè)和真實(shí)值比較可以得出,誤差相對(duì)小,證明模型合理,結(jié)果合理。