基于SEIR 的傳染病傳播模型

朱 柳 鄭婉婷 張貴層

(北方工業大學,北京100000)

縱觀歷史,一旦出現較大型的傳染病傳播情況,對世界的經濟發展和人民生活都會帶來了很大影響。因此想要預測和控制傳染病的蔓延,我們就必須了解一般傳染病毒的傳播機制,定量地研究傳染病的傳播規律,找出其傳播的數學模型。

本文結合實際影響因素建立了一個最接近現實情況的傳播模型,并通過建立的預測模型對病毒傳播未來一段時間發展態勢進行預測。

1 傳播模型的建立

1.1 傳染病傳播機制

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

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

α1為感染者I 的單位時間傳播率,即單位時間內一個感染者會導致α1S(t)/N 個健康人帶上病毒。其中帶病毒的人群有k1的比例抵抗力較差從而會直接患病成為感染者I,有k2比例因為自身有較強的抵抗力,沒有立即患病,但是在之后也有很大患病幾率,成為潛伏者E。

α2為潛伏者E 的單位時間傳播率,即單位時間內一個潛伏者會導致α2S(t)/N 個健康人帶上病毒,成為潛伏者E,因為潛伏者的傳染率比感染者的傳染率更弱,大部分潛伏者傳染別人使別人變成潛伏者E,而不會使得其直接變成感染者I。

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

β1為單位時間潛伏者E 變為感染者I 的比例,β2為單位時間潛伏者E 變為健康者S 的比例,變成健康者S 后仍然有被感染的風險。

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

圖1 傳染病傳播圖

1.2 病毒傳播模型微分方程

根據傳播機制建立關于s(t),e(t),i(t)和r(t)的方程:

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

1.2.2 等式左邊表示單位時間內健康者S 的變化,等式右邊從左向右第一部分是正,表示單位時間內潛伏者會有β2比例變成健康者。第二部分是負,表達單位時間未被隔離的感染者(比例為1-γ1)會感染健康者使得健康者變為感染者或潛伏者,傳播率為α1。第三部分是負,表達單位時間未隔離的潛伏者(比例為1-γ2)會感染健康者使得健康者變為潛伏者,傳播率為α2。

1.2.3 等式左邊表示單位時間內潛伏者E 的變化,等式右邊從左向右第一部分是正,表示單位時間未被隔離的感染者(比例為1-γ1)會感染健康者使得健康者變為潛伏者,比例為k2,傳播率為α1。第二部分是正,表達單位時間未被隔離的潛伏者(比例為1-γ2)會感染健康者使得健康者變為潛伏者,傳播率為α2。第三部分為負,表示單位時間內潛伏者會有β1比例變成感染者。第四部分為負,表示單位時間內潛伏者會有β2比例變成健康者。

1.2.4 等式左邊表示單位時間內感染者I 的變化,等式右邊從左向右第一部分是正,表示單位時間未隔離的感染者(比例為1-γ1)會感染健康者使得健康者變為感染者,比例為k1,感染率為α1。第二部分是正,表達單位時間內潛伏者會有β1比例變成感染者。第三部分為負,表示單位時間內有l1比例的感染者會被治愈成為排除者R,第四部分為負,表示單位時間內有l2比例的感染者會死亡成為排除者R。

1.2.5 等式左邊表單位時間內排除者R 的變化,等式右邊第一部分為正,表示單位時間內有l1比例的感染者會被治愈成為排除者R,第二部分為負,表示單位時間內有l2比例的感染者會死亡成為排除者R

1.2.6 單位時間內感染者會感染健康人群使得比例為k1的健康者變成感染者,比例為k2的健康者變成潛伏者,它們之和為1。

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

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

1.3 系數優化模型的建立

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

2 算例分析

2.1 傳播模型系數的估計

人口總N:經過網站[3]查詢某國2020 人口總數為143,964,709。

2.1.1 確定系數α1α2

α1:為未被隔離的感染者的傳播率,是一個重要參數,對其進行優化得出其值,α2為未被隔離的潛伏者的傳播率,由相關文獻[4]稱,4 成左右的潛伏者有感染能力,在此設α2=0.4×α1,α1范圍為:[0.01-1]。

2.1.2 確定系數γ1γ2

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

2.1.3 確定系數β1β2

β1為單位時間潛伏者變成感染者的比例,β2為單位時間潛伏者變為健康者的比例由查閱相關資料可得,最終潛伏者中有45%的比例變為感染者,那么55%的比例會變成健康者,基于目前的流行病學調查,潛伏期1-14 天,多為3-7 天,范圍為:[0.01-0.45]。

2.1.4 確定系數l1l2

根據真實數據用Matlab2018a 可繪出:

圖2 單位時間死亡率和治愈率曲線

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

圖3 單位時間治愈率擬合曲線圖

l1單位時間治愈率=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 確定系數k1k2

k1為被感染者感染的健康者變為感染者的比例,k2為被感染者感染的健康者變為潛伏者的概率。潛伏期為1～14 天,在此定義潛伏期2 天以下的就視為感染者,三天及以上視為潛伏者。

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

2.2 最佳系數優化解

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

α1=0.2923 β1=0.40

擬合值和真實值效果圖如圖4:

圖4 擬合值真實值比較圖

2.3 R0 的最終確定

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

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

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

根據相關文獻[6]可得:

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

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

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

TL表示感染者潛伏期的平均長度,SI 表示一個感染者被感染的時間和他感染的下一個人被感染的時間的間隔,TI=SI-TL,這三個數據都來自對現實病人情況的觀察,通過查詢文獻[7]可以估計SI 為8.4 天,TL為6 天。

最后剩一個未知數據r 為指數增長率,本文用擬合的無干預(即隔離率為零)數據來擬合r。

所以r=Ln1.1787=0.1644

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

3 病毒傳播數據的預測

傳播模型對病毒傳播數據的預測:

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

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

圖6 病毒傳播情況預測圖

列出相應表格可得出:

表1 病毒傳播情況預測表

同時可以看出在5.14-5.16,將會迎來病毒傳播的拐點,感染者人數會下降,新增確診感染者數也會下降。

4 模型結果分析

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

表2 誤差比較表

根據5.1、5.2 日預測和真實值比較可以得出,誤差相對小,證明模型合理,結果合理。