Python小白的數學建模課-B5. 新冠疫情 SEIR模型

youcans發表於2021-07-10
Python模型

傳染病的數學模型是數學建模中的典型問題，常見的傳染病模型有 SI、SIR、SIRS、SEIR 模型。

考慮存在易感者、暴露者、患病者和康復者四類人群，適用於具有潛伏期、治癒後獲得終身免疫的傳染病。

本文詳細給出了 SEIR 模型微分方程的建模、例程、結果和分析，讓小白都能懂。

『Python小白的數學建模課 @ Youcans』帶你從數模小白成為國賽達人。


1. SEIR 模型

1.1 SEIR 模型的提出

建立傳染病的數學模型來描述傳染病的傳播過程，要根據傳染病的發病機理和傳播規律， 結合疫情資料進行擬合分析，可以認識傳染病的發展趨勢，預測疫情持續時間和規模，分析和模擬各種防控措施對疫情發展的影響程度， 為傳染病防控工作提供決策指導，具有重要的理論意義和現實意義。

SI 模型是最簡單的傳染病傳播模型，把人群分為易感者（S 類）和患病者（I 類）兩類，通過 SI 模型可以預測傳染病高潮的到來；提高衛生水平、強化防控手段，降低病人的日接觸率，可以推遲傳染病高潮的到來。在 SI 模型基礎上發展的 SIS 模型考慮患病者可以治癒而變成易感者，SIS 模型表面傳染期接觸數 \(\sigma\) 是傳染病傳播和防控的關鍵指標，決定了疫情終將清零或演變為地方病長期存在。在 SI 模型基礎上考慮病癒免疫的康復者（R 類）就得到 SIR 模型，通過 SIR 模型也揭示傳染期接觸數 \(\sigma\) 是傳染病傳播的閾值，滿足 \(s_0>1/\sigma\) 才會發生傳染病蔓延，由此可以分析各種防控措施，如：提高衛生水平來降低日接觸率\(\lambda\)、提高醫療水平來提高日治癒率 \(\mu\)，通過預防接種達到群體免疫來降低 \(s_0\) 等。

傳染病大多具有潛伏期（incubation period），也叫隱蔽期，是指從被病原體侵入肌體到最早臨床症狀出現的一段時間。在潛伏期的後期一般具有傳染性。不同的傳染病的潛伏期長短不同，從短至數小時到長達數年，但同一種傳染病有固定的（平均）潛伏期。例如，流感的潛伏期為 1～3天，冠狀病毒感染的潛伏期為4~7天，新型冠狀病毒肺炎傳染病（COVID-19）的潛伏期為1-14天（* 來自：新型冠狀病毒肺炎診療方案試行第八版，潛伏時間 1～14天，多為3～7天，在潛伏期具有傳染性），肺結核的潛伏期從數週到數十年。

SEIR 模型考慮存在易感者（Susceptible）、暴露者（Exposed）、患病者（Infectious）和康復者（Recovered）四類人群，適用於具有潛伏期、治癒後獲得終身免疫的傳染病。易感者（S 類）被感染後成為潛伏者（E類），隨後發病成為患病者（I 類），治癒後成為康復者（R類）。這種情況更為複雜，也更為接近實際情況。

SEIR 模型的倉室結構示意圖如下：


1.2 SEIR 模型假設

  1. 考察地區的總人數 N 不變，即不考慮生死或遷移；

  2. 人群分為易感者（S 類）、暴露者（E 類）、患病者（I 類）和康復者（R 類）四類；

  3. 易感者（S 類）與患病者（I 類）有效接觸即變為暴露者（E 類），暴露者（E 類）經過平均潛伏期後成為患病者（I 類）；患病者（I 類）可被治癒，治癒後變為康復者（R 類）；康復者（R類）獲得終身免疫不再易感；

  4. 將第 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\)；

  5. 日接觸數 \(\lambda\)，每個患病者每天有效接觸的易感者的平均人數；

  6. 日發病率 \(\delta\)，每天發病成為患病者的暴露者佔暴露者總數的比例；

  7. 日治癒率 \(\mu\)，每天被治癒的患病者人數佔患病者總數的比例，即平均治癒天數為 \(1/\mu\)；

  8. 傳染期接觸數 \(\sigma = \lambda / \mu\)，即每個患病者在整個傳染期內有效接觸的易感者人數。


1.3 SEIR 模型的微分方程

由

\[\begin{align} & N \frac{ds}{dt} = - N \lambda s i\\ & N \frac{de}{dt} = N \lambda s i - N \delta e\\ & N \frac{di}{dt} = N \delta e - N \mu i\\ & N \frac{dr}{dt} = N \mu i\\ \end{align} \]

得：

\[\begin{cases} \begin{align*} & \frac{ds}{dt} = -\lambda s i, &s(0)=s_0\\ & \frac{de}{dt} = \lambda s i - \delta e, &e(0)=e_0\\ & \frac{di}{dt} = \delta e - \mu i, &i(0)=i_0 \end{align*} \end{cases} \]

SEIR 模型不能求出解析解，可以通過數值計算方法求解。


2. SEIR 模型的 Python 程式設計

2.1 Scipy 工具包求解微分方程組

SIS 模型是常微分方程初值問題，可以使用 Scipy 工具包的 scipy.integrate.odeint() 函式求數值解。

scipy.integrate.odeint(func, y0, t, args=())

**scipy.integrate.odeint() **是求解微分方程的具體方法，通過數值積分來求解常微分方程組。

odeint() 的主要引數：

  • func: callable(y, t, …) 　　導數函式 \(f(y,t)\) ，即 y 在 t 處的導數，以函式的形式表示
  • y0: array：　　初始條件 \(y_0\)，注意 SEIR模型是二元常微分方程組， 初始條件為陣列向量 \(y_0=[i_0, s_0]\)
  • t: array：　　求解函式值對應的時間點的序列。序列的第一個元素是與初始條件 \(y_0\) 對應的初始時間 \(t_0\)；時間序列必須是單調遞增或單調遞減的，允許重複值。
  • args: 嚮導數函式 func 傳遞引數。當導數函式 \(f(y,t,p1,p2,..)\) 包括可變引數 p1,p2.. 時，通過 args =(p1,p2,..) 可以將引數p1,p2.. 傳遞給導數函式 func。

odeint() 的返回值：

  • y: array 　　陣列，形狀為 (len(t),len(y0)，給出時間序列 t 中每個時刻的 y 值。

2.2 odeint() 求解 SEIR 模型的程式設計步驟

  1. 匯入 scipy、numpy、matplotlib 包。
  2. 定義導數函式 \(f(y,t)\)。注意對於常微分方程（例如 SI模型）和常微分方程組（SEIR模型），y 分別表示標量和向量，函式定義略有不同，以下給出兩種情況的例程以供對比。

常微分方程的導數定義（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 小白，又涉及到看著就心煩的微分方程組，所以我們寧願把程式寫得累贅一些，便於讀者將程式與前面的微分方程組逐項對應。

  1. 定義初值 \(y_0\) 和 \(y\) 的定義區間 \([t_0,\ t]\)，注意初值為陣列向量 \(y_0=[s_0,e_0,i_0]\)。
  2. 呼叫 odeint() 求 \(y\) 在定義區間 \([t_0,\ t]\) 的數值解。

SEIR 模型是三元常微分方程，返回值 y 是 len(t)*3 的二維陣列。


2.3 Python例程：SEIR 模型

# 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.4 SI /SIS/SIR 模型與SEIR 模型的比較

例程 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 模型。這表明具有潛伏期的傳染病，疫情發生和峰值的到來要晚於沒有潛伏期的傳染病，而且持續時間更長。


3. SEIR 模型引數的影響

SEIR 模型中有日接觸率 \(\lambda\) 、日發病率 \(\delta\) 和日治癒率 \(\mu\) 三個引數，還有 \(i_0、e_0、s_0\) 等初始條件，我們先用單因素分析的方法來觀察引數條件對於疫情傳播的影響。

3.1 初值條件 \(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。但這一現象是基於具體的引數條件的觀察，僅由該圖並不能確定其是否普遍規律。


3.2 日接觸率 \(\lambda\) 的影響

首先考察日接觸率 \(\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)\) 逐漸降低，但最終都趨於穩定。


3.3 日發病率 \(\delta\) 的影響

下面考察日發病率 \(\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)\) 逐漸降低，但最終都趨於穩定。


3.4 日治癒率 \(\mu\) 的影響

下面考察日治癒率 \(\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\) 才會發生傳染病蔓延。

這說明具有決定性影響的特徵引數，往往不是模型中的某個引數，而是多個引數特定關係的組合，僅從單因素實驗很難充分反映模型中的內在特徵。


4. SEIR 模型的相空間分析

4.1 Python例程：SEIR 模型的相軌跡

# 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 模型的相軌跡分析

上圖為例程 4.2 的執行結果，是 SEIR 模型的相軌跡圖。

圖中每一條 e-s 曲線，從直線 i(t)+s(t)=1 上的某一初值點出發，最終收斂於 s軸上的某一點，對應著某一個初值條件下的患病者與易感者比例隨時間的變化關係。

SEIR 模型的相軌跡圖比較複雜，難以在本文展開進行討論。但是，相軌跡圖中的虛線還是反映出了時間變化曲線圖中難以表達的一些重要特徵，以此為線索可以進行更深入的研究。


5. SEIR 模型結果討論

最後，我們簡單總結一下 SEIR 模型的特點：

  1. SEIR 模型是一個單向模型，易感者（S）不斷轉變為潛伏者（E），潛伏者（E）發病後成為患病者（I），患病者（I）不斷轉變為康復者（R），因此易感者比例 s(t) 單調遞減，康復者比例 r(i) 單調遞增。

  2. SEIR 模型假設潛伏期無傳染性，因此潛伏期延遲了傳染期的開始，疫情發生和峰值的到來要晚於沒有潛伏期的 SIR 模型。但潛伏期不會改變感染人群的累計數量，而且持續時間更長。

  3. \(1/\sigma\) 是傳染病蔓延的閾值，滿足 \(s_0>1/\sigma\) 才會發生傳染病蔓延。因此，為了控制傳染病的蔓延：一方面要提高閾值 \(1/\sigma\)，這可以通過提高衛生水平來降低日接觸率\(\lambda\)、提高醫療水平來提高日治癒率 \(\mu\)；另一方面要降低 \(s_0\)，這可以通過預防接種達到群體免疫來實現。

  4. 在 SEIR 模型的基礎上，可以根據不同傳染病病理特徵及疫情傳播特點，對模型進行進一步的改進，使模型與實際情況更加吻合，以便更準確地預測疫情發展趨勢。

  5. 在 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數模筆記-模擬退火演算法

相關文章