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

2021 年 7 月 4 日
筆記
Python 學習

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

SIS 模型型將人群分為 S 類和 I 類，考慮患病者可以治癒而變成易感者，但不考慮免疫期。

本文詳細給出了 SIS 模型的建模、常式、運行結果和模型分析，讓小白都能懂。

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

1. 疫情傳播 SIS 模型

傳染病動力學是對傳染病進行定量研究的重要方法。它依據種群繁衍遷移的特性、傳染病在種群內產生及傳播的機制、醫療與防控條件等外部因素，建立可以描述傳染病動力學行為的數學模型，通過對模型進行定性、定量分析和數值計算，模擬傳染病的傳播過程，預測傳染病的發展趨勢，研究防控策略的作用。

1.1 SI 模型

SI 模型把人群分為易感者（S類）和患病者（I類）兩類，易感者（S類）與患病者（I類）有效接觸即被感染，變為患病者，無潛伏期、無治癒情況、無免疫力。

SI 模型適用於只有易感者和患病者兩類人群，且無法治癒的疾病。

按照 SI 模型，最終所有人都會被傳染而變成病人，這是因為模型中沒有考慮病人可以治癒。因此只能是健康人患病，而患病者不能恢復健康（甚至也不會死亡，而是不斷傳播疫情），所以終將全部被傳染。

1.2 SIS 模型

SIS 模型將人群分為 S 類和 I 類，考慮患病者（I 類）可以治癒而變成易感者（S 類），但不考慮免疫期，因此患病者（I 類）治癒變成易感者以後還可以被感染而變成患病者。

SIS 模型適用於只有易感者和患病者兩類人群，可以治癒，但會反覆發作的疾病，例如腦炎、細菌性痢疾等治癒後也不具有免疫力的傳染病。

SIS 模型假設：

考察地區的總人數 N 不變，即不考慮生死或人口流動；
人群分為易感者（S類）和患病者（I類）兩類；
易感者（S類）與患病者（I類）有效接觸即被感染，變為患病者；患病者（I類）可被治癒而變為易感者，無潛伏期、無免疫力；
每個患病者每天有效接觸的易感者的平均人數（日接觸數）是 $\lambda$，稱為日接觸率；
每天被治癒的患病者人數占患病者總數的比例為 $\mu$ ，即日治癒率；
將第 t 天時 S類、I 類人群的佔比記為 $s(t)$、$i(t)$，數量為 $S(t)$、$I(t)$；初始日期 $t=0$ 時， S類、I 類人群佔比的初值為 $s_0$、$i_0$。

需要說明的是，不考慮生死或人口流動，通常是由於考慮一個封閉環境而且假定疫情隨時間的變化比生死、遷移隨時間的變化顯著得多, 因此後者可以忽略不計。

SIS 模型的微分方程：

由

\[\begin{align*}
N\frac{di}{dt} = N \lambda s i – N \mu i
\end{align*}
\]

得：

\[\begin{align*}
\frac{di}{dt} = \lambda i (1-i) – \mu i，\ i(0) = i_0
\end{align*}
\]

由日治癒率 $\mu$ 可知平均治癒天數為 $1/\mu$，也稱平均傳染期。定義 $\sigma = \lambda / \mu$，其含義是每個病人在傳染期內所傳染的平均人數，稱為傳染期接觸數。例如，平均傳染期 $1/\mu = 5$，日接觸率 $\lambda = 2$（每天傳染 2人），則傳染期接觸數 $\sigma = 10$。

SIS 模型的解析解為：

\[\begin{cases}
\begin{aligned}
& i(t)=\frac{i_0}{1+\lambda t i_0}&,\lambda = \mu\\
& i(t)=[\frac{\lambda}{\lambda-\mu} + (\frac{1}{i_0}-\frac{\lambda}{\lambda-\mu})*e^{-(\lambda – \mu) t}]^{-1} &,\lambda \neq \mu\\
\end{aligned}
\end{cases}\\
\]

注意：網上有些博文中解析解的公式誤寫成 $exp((\lambda-\mu)t)$ ，漏掉了一個負號。

2. SIS 模型的 Python 編程

2.1 Scipy 工具包求解 SIS 模型

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$，對於常微分方程組 $y_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 值。

odeint() 的編程步驟：

導入 scipy、numpy、matplotlib 包；
定義導數函數 $f(i,t)=\lambda i (1-i)- \mu i$ ；
定義初值 $y_0$ 和 $y$ 的定義區間 $[t_0,\ t]$；
調用 odeint() 求 $y$ 在定義區間 $[t_0,\ t]$ 的數值解。

2.2 Python常式：SIS 模型的解析解與數值解

# 1. SIS 模型，常微分方程，解析解與數值解的比較
from scipy.integrate import odeint  # 導入 scipy.integrate 模組
import numpy as np  # 導入 numpy包
import matplotlib.pyplot as plt  # 導入 matplotlib包

def dy_dt(y, t, lamda, mu):  # SIS 模型，導數函數
    dy_dt = lamda*y*(1-y) - mu*y  # di/dt = lamda*i*(1-i)-mu*i
    return dy_dt

# 設置模型參數
number = 1e5  # 總人數
lamda = 1.2  # 日接觸率, 患病者每天有效接觸的易感者的平均人數
sigma = 2.5  # 傳染期接觸數
mu = lamda/sigma  # 日治癒率, 每天被治癒的患病者人數占患病者總數的比例
fsig = 1-1/sigma
y0 = i0 = 1e-5  # 患病者比例的初值
tEnd = 50  # 預測日期長度
t = np.arange(0.0,tEnd,1)  # (start,stop,step)
print("lamda={}\tmu={}\tsigma={}\t(1-1/sig)={}".format(lamda,mu,sigma,fsig))

# 解析解
if lamda == mu:
    yAnaly = 1.0/(lamda*t +1.0/i0)
else:
    yAnaly= 1.0/((lamda/(lamda-mu)) + ((1/i0)-(lamda/(lamda-mu))) * np.exp(-(lamda-mu)*t))
# odeint 數值解，求解微分方程初值問題
ySI = odeint(dy_dt, y0, t, args=(lamda,0))  # SI 模型
ySIS = odeint(dy_dt, y0, t, args=(lamda,mu))  # SIS 模型

# 繪圖
plt.plot(t, yAnaly, '-ob', label='analytic')
plt.plot(t, ySIS, ':.r', label='ySIS')
plt.plot(t, ySI, '-g', label='ySI')

plt.title("Comparison between analytic and numerical solutions")
plt.axhline(y=fsig,ls="--",c='c')  # 添加水平直線
plt.legend(loc='best')  # youcans
plt.axis([0, 50, -0.1, 1.1])
plt.show()

2.3 SIS 模型解析解與數值解的比較

本圖為常式 2.2 的運行結果，圖中對解析解（藍色）與使用 odeint() 得到的數值解（紅色）進行比較。在該例中，無法觀察到解析解與數值解的差異，表明數值解的誤差很小。

本圖也比較了對相同日接觸率和患病者初值下 SI模型與 SIS模型進行了比較。SI 模型更早進入爆發期，最終收斂到 100%；SIS 模型下進入爆發期較晚，患病者的比例最終收斂到某個常數（與模型參數有關）。

考察 SI 模型與 SIS模型的關係，顯然 SI 模型是 SIS 模型在 $\mu = 0$ 時的特殊情況。

3. SIS 模型參數的影響

對於 SIS 模型，需要考慮日接觸率 $\lambda$ 與日治癒率 $\mu$ 的關係、患病者比例的初值 $i_0$ 的影響，總人數 N 沒有影響。

3.1 日接觸率 $\lambda$ 與日治癒率 $\mu$ 關係的影響

直觀地考慮，如果每天治癒的人數高於感染的人數，則疫情逐漸好轉，否則疫情逐漸嚴重。因此日接觸率 $\lambda$ 與日治癒率 $\mu$ 的關係非常關鍵，這就是傳染期接觸數 $\sigma = \lambda / \mu$ 的意義。

（1） $\sigma \leq 1$

當 $\sigma<1$ 時，傳染期接觸數小於 1，日接觸率小於日治癒率，患病率單調下降，最終清零，與患病率初值無關。 $\sigma$ 越小，疫情清零速度越快； $\sigma$ 越接近於 1，疫情清零越慢，但最終仍將清零。

分析其實際意義，傳染期接觸數小於 1，表明在傳染期內經過接觸而使易感者變成患病者的數量，小於在傳染期內治癒的患病者的數量，因此患病者數量、比例都會逐漸降低，所以最終可以清零，稱為無病平衡點。

當 $\sigma=1$ 時，不論患病率初值如何，患病率也是單調下降，最終趨近於 0。雖然在數學上患病率只能趨近於 0 而不等於 0，但考慮到總人數 N 是有限的，而患病者和易感者人數需要取整，因此 $\sigma=1$ 時最終也會清零。

（2） $\sigma > 1$

當 $\sigma>1$ 時，傳染期接觸數大於 1，日接觸率大於日治癒率，患病率的升降有兩種情況：

當患病率很低時，患病者人數少而易感者人數多，患病率上升；但隨著患病率增大，患病者越來越多而易感者越來越少，患病率雖然仍然上升但上升速度趨緩，最終趨於定值。
當患病率很高時，患病者人數多而易感者人數少，患病率下降；但隨著患病率減小，患病者越來越少而易感者越來越多，患病率雖然仍然下降但下降速度趨緩，最終也趨於相同的定值。
患病率最終都會收斂到穩態特徵值 $i_\infty=1-1/\sigma$。當 $i_0>i_\infty$ 即患病率初值大於穩態特徵值時，疫情曲線單調上升收斂；當 $i_0<i_\infty$ 即患病率初值小於穩態特徵值時，疫情曲線單調下降收斂；當 $i_0 = i_\infty$ 時，患病率始終大於穩態特徵值，疫情曲線為水平直線。

這表明，當 $\sigma>1$ 時疫情終將穩定但不會清零，而是長期保持一定的患病率，稱為地方病平衡點。

當 $\sigma=1$ 時，不論患病率初值如何，患病率都單調下降並最終趨於 0。

3.2 傳染期接觸數 $\sigma$ 與 $ di/dt$ 的關係

患病率的一階導數 $di/dt$ 的變化曲線，表明不論傳染期接觸數和初值如何，患病率的變化率都將收斂到 0，因此疫情終將穩定。當 $\sigma<1$ 時， $di/dt$ 始終是負值，單調上升趨近於 0；當 $\sigma>1$ 時， $di/dt$ 始終是正值，先上升達到峰值後再逐漸減小趨近於 0。

本圖為患病率 $i(t)$ 與一階導數 $di/dt$ 在不同傳染期接觸數下的關係曲線（相空間圖）。當 $\sigma\leq 1$ 時，曲線收斂到原點 $(0,0)$，即存在無病平衡點；當 $\sigma>1$ 時，曲線收斂到 $(1-1/\sigma,0)$，即存在地方病平衡點。

3.3 Python常式：傳染期接觸數 $\sigma$ 與 $ di/dt$ 的關係

# 4. SIS 模型，模型參數對 di/dt的影響
from scipy.integrate import odeint  # 導入 scipy.integrate 模組
import numpy as np  # 導入 numpy包
import matplotlib.pyplot as plt  # 導入 matplotlib包

def dy_dt(y, t, lamda, mu):  # SIS 模型，導數函數
    dy_dt = lamda*y*(1-y) - mu*y  # di/dt = lamda*i*(1-i)-mu*i
    return dy_dt

# 設置模型參數
number = 1e5  # 總人數
lamda = 1.2  # 日接觸率, 患病者每天有效接觸的易感者的平均人數
# sigma = np.array((0.1, 0.5, 0.8, 0.95, 1.0))  # 傳染期接觸數
sigma = np.array((0.5, 0.8, 1.0, 1.5, 2.0, 3.0))  # 傳染期接觸數
y0 = i0 = 0.05  # 患病者比例的初值
tEnd = 100  # 預測日期長度
t = np.arange(0.0,tEnd,0.1)  # (start,stop,step)

for p in sigma:
    ySIS = odeint(dy_dt, y0, t, args=(lamda,lamda/p))  # SIS 模型
    yDeriv = lamda*ySIS*(1-ySIS) - ySIS*lamda/p
    # plt.plot(t, yDeriv, '-', label=r"$\sigma$ = {}".format(p))
    plt.plot(ySIS, yDeriv, '-', label=r"$\sigma$ = {}".format(p)) #label='di/dt~i'
    print("lamda={}\tmu={}\tsigma={}\t(1-1/sig)={}".format(lamda,lamda/p,p,(1-1/p)))

# 繪圖
plt.axhline(y=0,ls="--",c='c')  # 添加水平直線
plt.title("i(t)~di/dt in SIS model") # youcans-xupt
plt.legend(loc='best')
plt.show()

4. SIS 模型結果討論

SIS 模型表明：

若 $\sigma > 1$，則 $\lim\limits_{t \to \infty} i(t) = 1-1/\sigma$，表明患病者始終存在，成為地方病。
若 $\sigma \leq 1$，則 $\lim\limits_{t \to \infty} i(t) = 0, (\sigma\leq 1)$ ，表明患病者人數不斷減少，最終可以清零。
SIS 模型說明，對於傳染病，需要對患病者進行隔離以減少有效接觸，通過減少日接觸率 $\lambda$ 來減小接觸數 $\sigma$ ，打破傳播鏈，最終控制疫情。

需要指出的是，本文討論的 SIS模型是把考察地區視為一個疫情均勻分布的整體進行研究。實際上，在考察區域的疫情分布必然是不均衡的，可能在局部區域發生疫情爆發導致該區域患病人數激增，是否會影響 SIS 模型的演化過程和穩定性呢？相關研究表明，擴散速度的不同可能導致種群空間分布的差異，在低風險區域將達到無病平衡點，在高風險區域仍將達到地方病平衡點。

【本節完】

歡迎關注『Python小白的數學建模課 @ Youcans』原創作品

原創作品，轉載必須標註原文鏈接：(//www.cnblogs.com/youcans/p/14968504.html)。

Crated：2021-06-09

歡迎關注『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數模筆記-模擬退火演算法

Tags: Python 學習