傳染病的數學模型是數學建模中的典型問題,常見的傳染病模型有 SI、SIR、SIRS、SEIR 模型。html
考慮存在易感者、暴露者、患病者和康復者四類人羣,適用於具備潛伏期、治癒後得到終身免疫的傳染病。python
本文詳細給出了 SEIR 模型微分方程的建模、例程、結果和分析,讓小白都能懂。算法
『Python小白的數學建模課 @ Youcans』帶你從數模小白成爲國賽達人。編程
創建傳染病的數學模型來描述傳染病的傳播過程,要根據傳染病的發病機理和傳播規律, 結合疫情數據進行擬合分析,能夠認識傳染病的發展趨勢,預測疫情持續時間和規模,分析和模擬各類防控措施對疫情發展的影響程度, 爲傳染病防控工做提供決策指導,具備重要的理論意義和現實意義。數組
SI 模型是最簡單的傳染病傳播模型,把人羣分爲易感者(S 類)和患病者(I 類)兩類,經過 SI 模型能夠預測傳染病高潮的到來;提升衛生水平、強化防控手段,下降病人的日接觸率,能夠推遲傳染病高潮的到來。在 SI 模型基礎上發展的 SIS 模型考慮患病者能夠治癒而變成易感者,SIS 模型表面傳染期接觸數 \(\sigma\) 是傳染病傳播和防控的關鍵指標,決定了疫情終將清零或演變爲地方病長期存在。在 SI 模型基礎上考慮病癒免疫的康復者(R 類)就獲得 SIR 模型,經過 SIR 模型也揭示傳染期接觸數 \(\sigma\) 是傳染病傳播的閾值,知足 \(s_0>1/\sigma\) 纔會發生傳染病蔓延,由此能夠分析各類防控措施,如:提升衛生水平來下降日接觸率\(\lambda\)、提升醫療水平來提升日治癒率 \(\mu\),經過預防接種達到羣體免疫來下降 \(s_0\) 等。ide
傳染病大多具備潛伏期(incubation period),也叫隱蔽期,是指從被病原體侵入肌體到最先臨牀症狀出現的一段時間。在潛伏期的後期通常具備傳染性。不一樣的傳染病的潛伏期長短不一樣,從短至數小時到長達數年,但同一種傳染病有固定的(平均)潛伏期。例如,流感的潛伏期爲 1~3天,冠狀病毒感染的潛伏期爲4~7天,新型冠狀病毒肺炎傳染病(COVID-19)的潛伏期爲1-14天(* 來自:新型冠狀病毒肺炎診療方案試行第八版,潛伏時間 1~14天,多爲3~7天,在潛伏期具備傳染性),肺結核的潛伏期從數週到數十年。函數
SEIR 模型考慮存在易感者(Susceptible)、暴露者(Exposed)、患病者(Infectious)和康復者(Recovered)四類人羣,適用於具備潛伏期、治癒後得到終身免疫的傳染病。易感者(S 類)被感染後成爲潛伏者(E類),隨後發病成爲患病者(I 類),治癒後成爲康復者(R類)。這種狀況更爲複雜,也更爲接近實際狀況。工具
SEIR 模型的倉室結構示意圖以下:spa
考察地區的總人數 N 不變,即不考慮生死或遷移;3d
人羣分爲易感者(S 類)、暴露者(E 類)、患病者(I 類)和康復者(R 類)四類;
易感者(S 類)與患病者(I 類)有效接觸即變爲暴露者(E 類),暴露者(E 類)通過平均潛伏期後成爲患病者(I 類);患病者(I 類)可被治癒,治癒後變爲康復者(R 類);康復者(R類)得到終身免疫再也不易感;
將第 t 天時 S 類、E 類、I 類、R 類人羣的佔比記爲 \(s(t)\)、\(e(t)\)、\(i(t)\)、\(r(t)\),數量分別爲 \(S(t)\)、\(E(t)\)、\(I(t)\)、\(R(t)\);初始日期 \(t=0\) 時,各種人羣佔比的初值爲 \(s_0\)、\(e_0\)、\(i_0\)、\(r_0\);
日接觸數 \(\lambda\),每一個患病者天天有效接觸的易感者的平均人數;
日發病率 \(\delta\),天天發病成爲患病者的暴露者佔暴露者總數的比例;
日治癒率 \(\mu\),天天被治癒的患病者人數佔患病者總數的比例,即平均治癒天數爲 \(1/\mu\);
傳染期接觸數 \(\sigma = \lambda / \mu\),即每一個患病者在整個傳染期內有效接觸的易感者人數。
由
得:
SEIR 模型不能求出解析解,能夠經過數值計算方法求解。
SIS 模型是常微分方程初值問題,可使用 Scipy 工具包的 scipy.integrate.odeint() 函數求數值解。
scipy.integrate.odeint(func, y0, t, args=())
**scipy.integrate.odeint() **是求解微分方程的具體方法,經過數值積分來求解常微分方程組。
odeint() 的主要參數:
odeint() 的返回值:
常微分方程的導數定義(SIS模型)
def dySIS(y, t, lamda, mu): # SIS 模型,導數函數 dy_dt = lamda*y*(1-y) - mu*y # di/dt = lamda*i*(1-i)-mu*i return dy_dt
常微分方程組的導數定義(SEIR模型)
def dySEIR(y, t, lamda, delta, mu): # SEIR 模型,導數函數 s, e, i = y ds_dt = - lamda*s*i # ds/dt = -lamda*s*i de_dt = lamda*s*i - delta*e # de/dt = lamda*s*i - delta*e di_dt = delta*e - mu*i # di/dt = delta*e - mu*i return np.array([ds_dt,de_dt,di_dt])
Python 能夠直接對向量、向量函數進行定義和賦值,使程序更爲簡潔。但考慮讀者主要是 Python 小白,又涉及到看着就心煩的微分方程組,因此咱們寧願把程序寫得累贅一些,便於讀者將程序與前面的微分方程組逐項對應。
SEIR 模型是三元常微分方程,返回值 y 是 len(t)*3 的二維數組。
# modelCovid4_v1.py # Demo01 of mathematical modeling for Covid2019 # SIR model for epidemic diseases. # Copyright 2021 Youcans, XUPT # Crated:2021-06-13 # Python小白的數學建模課 @ Youcans # 1. SEIR 模型,常微分方程組 from scipy.integrate import odeint # 導入 scipy.integrate 模塊 import numpy as np # 導入 numpy包 import matplotlib.pyplot as plt # 導入 matplotlib包 def dySIS(y, t, lamda, mu): # SI/SIS 模型,導數函數 dy_dt = lamda*y*(1-y) - mu*y # di/dt = lamda*i*(1-i)-mu*i return dy_dt def dySIR(y, t, lamda, mu): # SIR 模型,導數函數 s, i = y # youcans ds_dt = -lamda*s*i # ds/dt = -lamda*s*i di_dt = lamda*s*i - mu*i # di/dt = lamda*s*i-mu*i return np.array([ds_dt,di_dt]) def dySEIR(y, t, lamda, delta, mu): # SEIR 模型,導數函數 s, e, i = y # youcans ds_dt = -lamda*s*i # ds/dt = -lamda*s*i de_dt = lamda*s*i - delta*e # de/dt = lamda*s*i - delta*e di_dt = delta*e - mu*i # di/dt = delta*e - mu*i return np.array([ds_dt,de_dt,di_dt]) # 設置模型參數 number = 1e5 # 總人數 lamda = 0.3 # 日接觸率, 患病者天天有效接觸的易感者的平均人數 delta = 0.03 # 日發病率,天天發病成爲患病者的潛伏者佔潛伏者總數的比例 mu = 0.06 # 日治癒率, 天天治癒的患病者人數佔患病者總數的比例 sigma = lamda / mu # 傳染期接觸數 fsig = 1-1/sigma tEnd = 300 # 預測日期長度 t = np.arange(0.0,tEnd,1) # (start,stop,step) i0 = 1e-3 # 患病者比例的初值 e0 = 1e-3 # 潛伏者比例的初值 s0 = 1-i0 # 易感者比例的初值 Y0 = (s0, e0, i0) # 微分方程組的初值 # odeint 數值解,求解微分方程初值問題 ySI = odeint(dySIS, i0, t, args=(lamda,0)) # SI 模型 ySIS = odeint(dySIS, i0, t, args=(lamda,mu)) # SIS 模型 ySIR = odeint(dySIR, (s0,i0), t, args=(lamda,mu)) # SIR 模型 ySEIR = odeint(dySEIR, Y0, t, args=(lamda,delta,mu)) # SEIR 模型 # 輸出繪圖 print("lamda={}\tmu={}\tsigma={}\t(1-1/sig)={}".format(lamda,mu,sigma,fsig)) plt.title("Comparison among SI, SIS, SIR and SEIR models") plt.xlabel('t-youcans') plt.axis([0, tEnd, -0.1, 1.1]) plt.plot(t, ySI, 'cadetblue', label='i(t)-SI') plt.plot(t, ySIS, 'steelblue', label='i(t)-SIS') plt.plot(t, ySIR[:,1], 'cornflowerblue', label='i(t)-SIR') # plt.plot(t, 1-ySIR[:,0]-ySIR[:,1], 'cornflowerblue', label='r(t)-SIR') plt.plot(t, ySEIR[:,0], '--', color='darkviolet', label='s(t)-SIR') plt.plot(t, ySEIR[:,1], '-.', color='orchid', label='e(t)-SIR') plt.plot(t, ySEIR[:,2], '-', color='m', label='i(t)-SIR') plt.plot(t, 1-ySEIR[:,0]-ySEIR[:,1]-ySEIR[:,2], ':', color='palevioletred', label='r(t)-SIR') plt.legend(loc='right') # youcans plt.show()
例程 2.3 的參數和初值爲:\(\lambda=0.3,\delta=0.03,\mu=0.06,(s_0,e_0,i_0)=(0.001,0.001,0.998)\),上圖爲例程的運行結果。
曲線 i(t)-SI 是 SI 模型的結果,患病者比例急劇增加到 1.0,全部人都被傳染而變成患病者。
曲線 i(t)-SIS 是 SIS 模型的結果,患病者比例快速增加並收斂到某個常數,即穩態特徵值 \(i_\infty=1-\mu/\lambda = 0.8\),代表疫情穩定,並將長期保持必定的患病率,稱爲地方病平衡點。
曲線 i(t)-SIR 是 SIR 模型的結果,患病者比例 i(t) 先上升達到峯值,而後再逐漸減少趨近於常數。
曲線 s(t)-SEIR、e(t)-SEIR、i(t)-SEIR、r(t)-SEIR 分別表示 SEIR模型中易感者(S類)、潛伏者(E類)、患病者(I類)和康復者(R 類)人羣的佔比。
圖中易感者比例 s(t) 單調遞減並收斂到接近於 0 的穩定值。潛伏者比例 e(t) 曲線存在波峯,先逐漸上升而達到峯值,而後再逐漸減少,最終趨於 0。患病者比例 i(t) 曲線與潛伏者比例曲線相似,上升達到峯值後逐漸減少,最終趨於 0;但患病者比例曲線發展、達峯的時間比潛伏者曲線要晚一些,峯值強度也較低。康復者比例 r(i) 單調遞增並收斂到非零的穩態值。以上分析只是對本圖進行的討論,並不是廣泛結論,取決於具體參數條件。
比較相同參數條件下 SIR 和 SEIR 模型的結果,SIR 模型中患病者比例 i(t) 的波形起點、峯值和終點到來的時間都顯著早於 SEIR 模型,峯值強度也高於 SEIR 模型。這代表具備潛伏期的傳染病,疫情發生和峯值的到來要晚於沒有潛伏期的傳染病,並且持續時間更長。
SEIR 模型中有日接觸率 \(\lambda\) 、日發病率 \(\delta\) 和日治癒率 \(\mu\) 三個參數,還有 \(i_0、e_0、s_0\) 等初始條件,咱們先用單因素分析的方法來觀察參數條件對於疫情傳播的影響。
SEIR 模型中有 \(i_0、e_0、s_0\) 等 3個初始條件,組合衆多沒法窮盡。考慮實際狀況中,疫情初始階段尚無康復者,而潛伏者比例每每高於確診的發病者,咱們假定 \(e_0/i_0=2,r_0=0\),考察不一樣 \(i_0\) 時的疫情傳播狀況。
經過對於該參數下不一樣的患病者、潛伏者比例初值條件的模擬,能夠看到患病者、潛伏者比例的初值條件對疫情發生、達峯、結束的時間遲早具備直接影響,但對疫情曲線的形態和特徵影響不大。不一樣初值條件下的疫情曲線,幾乎是沿着時間指標平移的。
這說明若是不進行治療防控等人爲干預,疫情傳播過程與初始患病者、潛伏者比例關係並不大,該來的總會來。
圖中患病率達到高峯後逐步下降,直至趨近於 0;易感率在疫情爆發後迅速降低,直至趨近於 0。但這一現象是基於具體的參數條件的觀察,僅由該圖並不能肯定其是否廣泛規律。
首先考察日接觸率 \(\lambda\) 的影響。
保持參數 \(\delta =0.1,\mu=0.06, (i_0=0.001,e_0=0.002,s_0=0.997)\) 不變,$\lambda = [0.12, 0.25, 0.5, 1.0, 2.0] $ 時 \(i(t), s(t)\) 的變化曲線以下圖所示。
經過對於該條件下日接觸率的單因素分析,能夠看到隨着日接觸率 \(\lambda\) 的增大,患病率比例 \(i(t)\) 出現的峯值更早、更強,而易感者比例 \(s(t)\) 逐漸下降,但最終都趨於穩定。
下面考察日發病率 \(\delta\) 的影響。保持參數 \(\lambda =0.25,\mu=0.06, (i_0=0.001,e_0=0.002,s_0=0.997)\) 不變,$\delta = [0.02, 0.05, 0.1, 0.2, 0.4] $ 時 \(i(t), s(t)\) 的變化曲線以下圖所示。
經過對於該條件下日接觸率的單因素分析,能夠看到隨着日接觸率 \(\lambda\) 的增大,患病率比例 \(i(t)\) 出現的峯值更早、更強,而易感者比例 \(s(t)\) 逐漸下降,但最終都趨於穩定。
下面考察日治癒率 \(\mu\) 的影響。保持參數 \(\lambda =0.25,\delta=0.1, (i_0=0.001,e_0=0.002,s_0=0.997)\) 不變,$\mu = [0.02, 0.05, 0.1, 0.2, 0.4] $ 時 \(i(t), s(t)\) 的變化曲線以下圖所示。
經過對該條件下日治癒率的單因素分析,能夠看到在日治癒率 \(\mu=0.4\) 時,患病者比例始終很是低並趨於 0,易感者比例幾乎不變,代表疫情不會傳播,這是由於患病者治癒的速度快於易感者被感染的速度。
在日治癒率 \(\mu=0.2\) 時,患病者比例也始終很是低(接近 0),易感者比例緩慢下降並趨於穩定值,代表疫情的傳播十分緩慢微弱,這是由於患病者治癒的速度較快,在易感者比例逐漸下降到某一水平後治癒人數與感染人數將達到平衡。
在日治癒率較低時 (\(\mu<0.2\) ),患病者比例曲線存在波峯,而後再逐漸減少,最終趨於 0。隨着日治癒率 \(\mu\) 的減少,患病率比例 \(i(t)\) 曲線的峯值更強、也稍早。易感者比例 \(s(t)\) 隨着患病者比例上升而迅速下降,最終趨於某個穩定值,也達到治癒與感染的平衡。
經過對不一樣參數的單因素實驗和結果分析,能夠發現雖然各個模型參數和初始條件都會影響疫情曲線,但日治癒率 與日接觸率的關係更爲重要。進一步的分析和模擬能夠發現,傳染期接觸數 \(\sigma = \lambda / \mu\)是傳染病蔓延的閾值,知足 \(s_0>1/\sigma\) 纔會發生傳染病蔓延。
這說明具備決定性影響的特徵參數,每每不是模型中的某個參數,而是多個參數特定關係的組合,僅從單因素實驗很難充分反映模型中的內在特徵。
# modelCovid4_v1.py # Demo01 of mathematical modeling for Covid2019 # SIR model for epidemic diseases. # Copyright 2021 Youcans, XUPT # Crated:2021-06-15 # Python小白的數學建模課 @ Youcans # 7. SEIR 模型,常微分方程組 相空間分析: e(t)~i(t) from scipy.integrate import odeint # 導入 scipy.integrate 模塊 import numpy as np # 導入 numpy包 import matplotlib.pyplot as plt # 導入 matplotlib包 def dySEIR(y, t, lamda, delta, mu): # SEIR 模型,導數函數 s, e, i = y ds_dt = -lamda*s*i # ds/dt = -lamda*s*i de_dt = lamda*s*i - delta*e # de/dt = lamda*s*i - delta*e di_dt = delta*e - mu*i # di/dt = delta*e - mu*i return np.array([ds_dt,de_dt,di_dt]) # 設置模型參數 number = 1e5 # 總人數 lamda = 0.3 # 日接觸率, 患病者天天有效接觸的易感者的平均人數 delta = 0.1 # 日發病率,天天發病成爲患病者的潛伏者佔潛伏者總數的比例 mu = 0.1 # 日治癒率, 天天治癒的患病者人數佔患病者總數的比例 sigma = lamda / mu # 傳染期接觸數 tEnd = 500 # 預測日期長度 t = np.arange(0.0, tEnd, 1) # (start,stop,step)# e0List = np.arange(0.01,0.4,0.05) # (start,stop,step) e0List = np.arange(0.01,0.4,0.05) # (start,stop,step) for e0 in e0List: # odeint 數值解,求解微分方程初值問題 i0 = 0 # 潛伏者比例的初值 s0 = 1 - i0 - e0 # 易感者比例的初值 ySEIR = odeint(dySEIR, (s0,e0,i0), t, args=(lamda,delta,mu)) # SEIR 模型 plt.plot(ySEIR[:,1], ySEIR[:,2]) # (e(t),i(t)) print("lamda={}\tdelta={}\mu={}\tsigma={}\ti0={}\te0={}".format(lamda,delta,mu,lamda/mu,i0,e0)) # 輸出繪圖 plt.title("Phase trajectory of SEIR models: e(t)~i(t)") plt.axis([0, 0.4, 0, 0.4]) plt.plot([0,0.4],[0,0.35],'y--') #[x1,x2][y1,y2] plt.plot([0,0.4],[0,0.18],'y--') #[x1,x2][y1,y2] plt.text(0.02,0.36,r"$\lambda=0.3, \delta=0.1, \mu=0.1$",color='black') plt.xlabel('e(t)') plt.ylabel('i(t)') plt.show()
上圖爲例程 4.2 的運行結果,是 SEIR 模型的相軌跡圖。
圖中每一條 e-s 曲線,從直線 i(t)+s(t)=1 上的某一初值點出發,最終收斂於 s軸上的某一點,對應着某一個初值條件下的患病者與易感者比例隨時間的變化關係。
SEIR 模型的相軌跡圖比較複雜,難以在本文展開進行討論。可是,相軌跡圖中的虛線仍是反映出了時間變化曲線圖中難以表達的一些重要特徵,以此爲線索能夠進行更深刻的研究。
最後,咱們簡單總結一下 SEIR 模型的特色:
SEIR 模型是一個單向模型,易感者(S)不斷轉變爲潛伏者(E),潛伏者(E)發病後成爲患病者(I),患病者(I)不斷轉變爲康復者(R),所以易感者比例 s(t) 單調遞減,康復者比例 r(i) 單調遞增。
SEIR 模型假設潛伏期無傳染性,所以潛伏期延遲了傳染期的開始,疫情發生和峯值的到來要晚於沒有潛伏期的 SIR 模型。但潛伏期不會改變感染人羣的累計數量,並且持續時間更長。
\(1/\sigma\) 是傳染病蔓延的閾值,知足 \(s_0>1/\sigma\) 纔會發生傳染病蔓延。所以,爲了控制傳染病的蔓延:一方面要提升閾值 \(1/\sigma\),這能夠經過提升衛生水平來下降日接觸率\(\lambda\)、提升醫療水平來提升日治癒率 \(\mu\);另外一方面要下降 \(s_0\),這能夠經過預防接種達到羣體免疫來實現。
在 SEIR 模型的基礎上,能夠根據不一樣傳染病病理特徵及疫情傳播特色,對模型進行進一步的改進,使模型與實際狀況更加吻合,以便更準確地預測疫情發展趨勢。
在 SEIR 模型的基礎上,能夠結合實際的疫情數據來擬合和估計模型參數,進而用來模擬和分析不一樣治療方案和防控措施對疫情發展的影響,爲新冠疫情的防控工做提供決策指導。
【本節完】
版權聲明:
歡迎關注『Python小白的數學建模課 @ Youcans』 原創做品
原創做品,轉載必須標註原文連接:。
Copyright 2021 Youcans, XUPT
Crated:2021-06-15
歡迎關注 『Python小白的數學建模課 @ Youcans』,每週更新數模筆記
Python小白的數學建模課-01.新手必讀
Python小白的數學建模課-02.數據導入
Python小白的數學建模課-03.線性規劃
Python小白的數學建模課-04.整數規劃
Python小白的數學建模課-05.0-1規劃
Python小白的數學建模課-06.固定費用問題
Python小白的數學建模課-07.選址問題
Python小白的數學建模課-09.微分方程模型
Python小白的數學建模課-B2.新冠疫情 SI模型
Python小白的數學建模課-B3.新冠疫情 SIS模型
Python數模筆記-PuLP庫
Python數模筆記-StatsModels統計迴歸
Python數模筆記-Sklearn
Python數模筆記-NetworkX
Python數模筆記-模擬退火算法