SIR及SEIR建模的简单示例

领域模型 专栏收录该内容
1 篇文章 0 订阅

概述

看了一些2019-nCoV相关的文章,摘录并总结了一些关于SIR和SEIR模型的定义。并通过代码进了简单模型的实现,使用R语言作为编程工具。

1.一些定义

1.1 一些名词

序号词汇解释
1Pneumonia肺炎
2Coronavirus冠状病毒
3Incubation潜伏
4Quarantined隔离
5Susceptible易感的,这里用于指易感人群,在CDC的公告中,所有人都属于易感人群。
6zoonotic动物传播的
7Reproductive再生,复制。这里指感染者再传播并产生新的感染者。
8novel新的
92019-nCoV2019新型冠状病毒

1.2 一些符号

序号符号含义
1 R 0 R_0 R0基础再传播人数
2 τ \tau τ传染概率,一名易感者和一名感染者接触时被感染的概率
3 c ˉ \bar{c} cˉ单位时间内和感染者接触的易感人员的平均比例值
4 d d d感染暴露时长
5 S S S易感人群
6 I I I感染人群
7 R R R移除人群,是指,被隔离或治愈而消除影响的感染人群
8E暴露人群, 处于潜伏期

1.3 一些定义

R 0 = τ ⋅ c ˉ ⋅ d R_0 = \tau \cdot \bar{c} \cdot d R0=τcˉd

用1小时表示1单位时长,通俗的理解:
再 传 播 人 数 = 传 染 率 ⋅ 接 触 易 感 人 数 的 平 均 值 ⋅ 感 染 暴 露 时 长 再传播人数=传染率\cdot 接触易感人数的平均值 \cdot 感染暴露时长 =
假定 τ = 0.2 , c ˉ = 5 , d = 4 \tau=0.2,\bar{c}=5,d=4 τ=0.2,cˉ=5,d=4,则 R 0 = 4 R_0=4 R0=4

2.方法论

2.1 SIR

SIR(Susceptible-Infected-Removed)模型,用于传染病的传播人数进行建模。

关于模型的一些定义和假设:

  1. N,封闭系统内的总人口。
  2. 传染率、排除率等为常量。
  3. 不考虑出生和自然死亡。
  4. 人群是均匀混合的。任何感染者可以以概率接触任何一名易感者,这里的概率可以用总体均值来替代。

通过模型的假设我们可以看到,SIR实际上是对传染病流行早期传播行为的建模。模型形式如下:
{ d s d t = − β s i d i d t = β s i − γ i d r d t = γ i \left \{ \begin{aligned} & \frac{ds}{dt} = -\beta s i \\ & \frac{di}{dt} = \beta s i -\gamma i \\ & \frac{dr}{dt} = \gamma i \end{aligned} \right. dtds=βsidtdi=βsiγidtdr=γi
其中 β = τ c ˉ \beta = \tau\bar{c} β=τcˉ表示有效接触率, v v v表示消除率, $\gamma 表 示 移 除 率 , 都 是 常 量 。 因 而 传 染 暴 露 时 长 实 际 上 是 表示移除率,都是常量。因而传染暴露时长实际上是 v 的 倒 数 , 有 的倒数,有 d=\gamma^{-1}$。

当传染病爆发时,感染者人数随着时间上升,因此有 d i / d t ≥ 0 di/dt\geq0 di/dt0,从而有:
β s i − γ i > β s i γ > i \beta s i-\gamma i\gt \\ \frac{\beta s i}{\gamma}>i βsiγi>γβsi>i
在一场疫情的爆发之初,每个人都可视为易感染群,因此这里 s ≈ 1 s\approx 1 s1,代入上式有:
β γ = τ c ˉ d = R 0 > 1 \frac{\beta}{\gamma}=\tau\bar{c}d=R_0>1 γβ=τcˉd=R0>1
因而当 R 0 > 1 R_0>1 R0>1时疫情是处于传染阶段的。

对上述微分方程组求解得到:
I = ( S 0 + I 0 ) − S + 1 R 0 ln ⁡ S S 0 I=(S_0+I_0)-S+\frac{1}{R_0}\ln\frac{S}{S_0} I=(S0+I0)S+R01lnS0S
其中 S 0 , I 0 S_0,I_0 S0,I0表示初值。当 S 0 < 1 R 0 S_0<\frac{1}{R_0} S0<R01时传染降低,为了达到这个目的可采取以下措施:

  1. 降低 S 0 S_0 S0,即减少易感者数量,在这里也就是戴口罩、不外出等。
  2. 提高 1 R 0 \frac{1}{R_0} R01,即降低$\beta $(接触率),在这里也就是减少人员流动。
  3. 提高 γ \gamma γ(移除率),提高隔离率或治愈率。

2.2 SEIR

SEIR(Susceptible-Exposed-Infected-Removed),类似于SIR,但是增加了对潜伏期的定义,因此更适用于具有一定潜伏期的传染病。状态之间的转化如下所示:
在这里插入图片描述

其中, λ \lambda λ表示易感人群的输入(人口增加), μ \mu μ表示死亡率, k k k表示从暴露人群到确诊感染者的比率, γ \gamma γ是感染者的移除率。模型由如下四个等式组成:
S ˙ = − β S I + λ − μ S E ˙ = β S I − ( μ + k ) E I ˙ = k E − ( γ + μ ) I R ˙ = γ I − μ R \begin{aligned} \dot{S} &=-\beta S I+\lambda-\mu S \\ \dot{E} &=\beta S I-(\mu+k) E \\ \dot{I} &=k E-(\gamma+\mu) I \\ \dot{R} &=\gamma I-\mu R \end{aligned} S˙E˙I˙R˙=βSI+λμS=βSI(μ+k)E=kE(γ+μ)I=γIμR

2.3 代际传播

2.3.1 传播矩阵

对于同一种传染病而言,每一个患者的感染途径可能是不同的,比如:蝙蝠传人、男人传女人、狗传人等。定义代际传播矩阵 G \mathbf{G} G,其中元素 g i , j g_{i,j} gi,j表示下一代中由一个 j j j类病患导致的 i i i类病患的数量,由此可知 G \mathbf{G} G是方阵。同时 R 0 R_0 R0对应 G G G的谱半径(绝对值最大的特征值)。从 G \mathbf{G} G的数学性质来看,它是非奇异的,同时具有一个正的特征值并且严格大于其它特征值,实际上这个特征值就是 R 0 R_0 R0

对于只有两种状态的传染病而言有:
G = [ a b c d ] \mathbf{G}=\left[ \begin{array}{ll}{a} & {b} \\ {c} & {d} \end{array}\right] G=[acbd]
其特征值为: λ ± = a + d 2 ± ( ( a + d ) / 2 ) 2 − ( a d − b c ) \lambda_{\pm}=\frac{a+d}{2}\pm\sqrt{((a+d)/2)^2-(ad-bc)} λ±=2a+d±((a+d)/2)2(adbc)

此外, G \mathbf{G} G也可以写成入下形式:
G = F V − 1 \mathbf{G}=FV^{-1} G=FV1
其中 F = [ ∂ F i ( x 0 ) ∂ x j ] F=\left[\frac{\partial F_i(x_0)}{\partial x_j}\right] F=[xjFi(x0)]表示新增感染者, V = [ ∂ V i ( x 0 ) ∂ x j ] V=\left[\frac{\partial V_i(x_0)}{\partial x_j}\right] V=[xjVi(x0)]表示感染者的区间传播。

3 模型实现

编码采用R语言,对于SIR系列的建模,R语言中有现成的软件包:SimInf。为了展示数据的迭代过程,这里并不打算直接采用。

3.1 参数设定

结合当前实际情况,有如下设定:

  • N,以武汉为中心的辐射人口:19 000 000。
  • μ \mu μ, 患病死亡率:0.02。
  • I 0 I_0 I0,期初患病人数:1。
  • β γ \frac{\beta}{\gamma} γβ, 2.68。

3.2 SIR

library(deSolve) 
library(ggplot2)

(1)模型

基于如下微分方程组构建模型
d s d t = − β s i d i d t = β s i − γ i d r d t = γ i \begin{aligned} & \frac{ds}{dt} = -\beta s i \\ & \frac{di}{dt} = \beta s i -\gamma i \\ & \frac{dr}{dt} = \gamma i \end{aligned} dtds=βsidtdi=βsiγidtdr=γi

sir <- function(time, state, pars) {
  with(as.list(c(state, pars)), {
    dS <- -beta * S * I/N
    dI <- beta * S * I/N - gamma * I
    dR <- gamma * I
    return(list(c(dS, dI, dR)))
  })
}

(2)参数

N <- 1.9e8 # 总人口
I0 <- 1 # 初始感染者数量
RM0 <- 0 # 初始移除人员数量
S0 <- N - I0 - RM0 # 初始易感人群数量
init <- c(S = S0, I = I0, R = RM0) # 初始值
# 以下参数在模型假定下是常量
pars <- c(
  beta = 0.55, # 有效接触率
  gamma = 0.2, # 移除率
  N = N # 人口
  ) 
# 迭代次数,以天计
times <- seq(0, 150, by = 1) 

(3)计算

res <- as.data.frame(ode(y = init, times = times, func = sir, parms = pars))

(4)绘图

ggplot(res) +
  geom_line(aes(x = time, y = S, col = '易感'))+
  geom_line(aes(x = time, y = I, col = '感染'))+
  geom_line(aes(x = time, y = R, col = '移除'))+
  theme_light(base_family = 'Kai') +
  scale_colour_manual("",
  values=c("易感" = "cornflowerblue", "感染" = "darkred", "移除" = "forestgreen")
  ) +
  scale_y_continuous('')

在这里插入图片描述

3.3 SEIR

(1)模型

基于如下微分方程组构建模型
d s d t = − β S I N d e d t = β S I N − k E d I d t = k E − ( γ + μ ) I R d t = γ I \begin{aligned} \frac{ds}{dt} &= -\beta S \frac{I}{N} \\ \frac{de}{dt} &= \beta S \frac{I}{N}-k E \\ \frac{dI}{dt} &=k E-(\gamma+\mu) I \\ \frac{R}{dt} &=\gamma I \end{aligned} dtdsdtdedtdIdtR=βSNI=βSNIkE=kE(γ+μ)I=γI

seir<-function(time, state, pars){ 
  with(as.list(c(state, pars)),{ 
    dS <-- S * beta * I/N 
    dE <- S * beta * I/N - E * k 
    dI <- E * k - I * (mu + gamma) 
    dR <- I * gamma
    dN <- dS + dE + dI + dR 
    
    list(c(dS,dE,dI,dR,dN)) 
  }) 
} 

(2)参数

N <- 1.9E8 # 总人口
I0 <- 89 # 期初感染数
E0 <- 0 # 期初潜伏数
RM0 <- 0 # 期初移除数
S0 = N - I0 - RM0 # 期初易感人数
init<-c(S = S0, E = E0, I = I0, R = RM0, N = N)	
time <- seq(0, 150, 1) 
pars<-c( 
  beta = 0.55,	#有效接触率
  k = 1,	#潜伏到感染的转化率 
  gamma = 0.2,	#RECOVERY 
  mu=0.02	#感染期死亡率 
) 

(3)计算

res.seir<-as.data.frame(lsoda(y = init, times = time, func = seir, parms = pars)) 

(4)绘图

ggplot(res.seir) +
  geom_line(aes(x = time, y = S, col = '2 易感'))+
  geom_line(aes(x = time, y = E, col = '3 潜伏'))+
  geom_line(aes(x = time, y = I, col = '4 感染'))+
  geom_line(aes(x = time, y = R, col = '5 移除'))+
  geom_line(aes(x = time, y = N, col = '1 人口'))+
  theme_light(base_family = 'Kai') +
  scale_colour_manual("",
  values=c(
    "2 易感" = "cornflowerblue", "3 潜伏" = "orange",
    "4 感染" = "darkred", "5 移除" = "forestgreen", 
    "1 人口" = "black"
    )
  ) +
  scale_y_continuous('')

在这里插入图片描述

3.4 小结

对比SIR和SEIR模型的结果,可以看到,相同条件下具有潜伏期的疾病其感染人数峰值的到来要晚于没有潜伏期的疾病,并且持续时间更长。

参考文献

  1. Joseph T Wu*, Kathy Leung*, Gabriel M Leung. Nowcasting and forecasting the potential domestic and
    international spread of the 2019-nCoV outbreak originating
    in Wuhan, China: a modelling study.Lancet,2020.
  2. Gerardo Chowell.Fitting dynamic models to epidemic outbreaks with quantified uncertainty: A primer for parameter uncertainty, identifiability, and forecasts.Infectious Disease Modelling,2017.
  3. James Holland Jones.Notes On R0.2007.
©️2022 CSDN 皮肤主题:点我我会动 设计师:我叫白小胖 返回首页

打赏作者

禾雨辰子

你的鼓励将是我创作的最大动力

¥2 ¥4 ¥6 ¥10 ¥20
输入1-500的整数
余额支付 (余额:-- )
扫码支付
扫码支付:¥2
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值