廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

感染症のモデル(SIRモデル)に入門した

導入

 記法はウィキペディアに合わせる.
 SIRモデル - Wikipedia

 時刻 t において,

  • S(t):感染可能者(Susceptible)の数. これから病気にかかるおそれのある人たち.
  • I(t):感染者(Infected)の数. いま病気にかかっていて人に病気を移す可能性がある人たち.
  • R(t):除外者(Recovered)の数. すでに病気にかかって隔離されている, または回復して免疫ができている, あるいは死んでいる人たち.

とする.

f:id:abrahamcow:20140713151357p:plain

 全体の人数は N(定数)とする.
 S(t)+I(t)+R(t)=N \tag{1}.

  • S と S が接触 → 感染者増えない
  • S と I が接触 → 感染者増える
  • S と R が接触 → 感染者増えない

なので, SSI に比例する速度で減少するとすると, 感染症の流行過程は下記の微分方程式で記述できる.

\displaystyle \frac{d}{dt}S(t) = -\beta S(t) I(t) \tag{2}
\displaystyle \frac{d}{dt}I(t) = \beta S(t) I(t) - \gamma I(t) \tag{3}
\displaystyle \frac{d}{dt}R(t) = \gamma I(t) \tag{4}

ここで定数βは伝染速度, γは除外速度と解釈できる.

(1) 式  S(t)+I(t)+R(t)=N は, (2) + (3) + (4) を積分した結果であることに留意する.

\displaystyle \frac{d}{dt}S(t) + \frac{d}{dt}I(t)+ \frac{d}{dt}R(t) =0.

I について解く.

(1), (2) より,

\displaystyle \frac{\frac{d}{dt}I(t)}{\frac{d}{dt}S(t)}=\frac{\beta S(t) I(t) - \gamma I(t)}{-\beta S(t) I(t)} = -1 + \gamma /\beta,

\displaystyle \frac{dI}{dS} = -1 + \gamma /\beta,
S積分すると,
\displaystyle I=-S+(\gamma /\beta) \log S +C,
C積分定数.

S(0)=N-I(0)(すなわちR(0) =0)とすると, 初期値は
\displaystyle C=N-(\gamma /\beta) \log S(0).

よって,

\displaystyle I(t)=N-S(t)+(\gamma /\beta)\log\{S(t)/S(0)\} \tag{5}.

S について解く.

(5) 式より,

\displaystyle N-S(t) -I(t)= -(\gamma /\beta)\log\{S(t)/S(0)\}
\displaystyle  R(t) = -(\gamma /\beta)\log\{S(t)/S(0)\}
\displaystyle  S(t) = S(0) \exp \left\{ - \frac{R(t)}{\gamma /\beta}\right\} \tag{6}.

R について解く.

 残念ながら R について解くことはできないらしい.

 しかし, 近似解は出せる.
 (1), (4) より
\displaystyle \frac{d}{dt}R(t) = \gamma(N(t)-R(t)-S(t)),
これに (6) を代入して,
\displaystyle \frac{d}{dt}R(t) = \gamma \left(N(t)-R(t)-S(0) \exp \left\{ - \frac{R(t)}{\gamma /\beta}\right\} \right).

マクローリン展開
 \exp \left\{ - \frac{R(t)}{\gamma /\beta}\right\} \approx 1-\frac{R(t)}{\gamma /\beta}+\frac{(\frac{R(t)}{\gamma /\beta})^2}{2}+ \cdots

により近似した方程式,
\displaystyle \frac{d}{dt}R(t) = \gamma \left( N(t)-R(t)-S(0)\left\{1-\frac{R(t)}{\gamma /\beta}+\frac{(\frac{R(t)}{\gamma /\beta})^2}{2}\right\} \right) .
は解,
\displaystyle R(t) = (\gamma /\beta)^2 [S(0)/ (\gamma /\beta) - 1 + \alpha \tanh (\alpha \gamma t/r -\phi)]/S(0) .
を持つ(らしい). ここで,
\displaystyle \alpha = \left\{\left(\frac{S(0)}{\gamma /\beta} -1\right)^2 + \frac{2S(0)(N-S(0))}{(\gamma /\beta)^2}\right\},
\displaystyle \phi=\tanh^{-1}\left\{\frac{S(0)/(\gamma/\beta)-1}{\alpha}\right\}.

(検算はしていません.)

参考文献・参考になりそうな文献

微分方程式で数学モデルを作ろう

微分方程式で数学モデルを作ろう

付録

今回図は Graphviz を使って書いてみた.

digraph SIR{
  graph [rankdir = LR];
  S -> I[label = "β"]
  I -> R[label = "γ"];
}

追記:数値計算

 R では deSolve パッケージを使うと微分方程式が数値的に解けるみたいです.

http://heartruptcy.blog.fc2.com/blog-entry-164.html

が, 使いこなせるかなあ...

 βγ には適当な数字を与えています.

library(deSolve)

parameters <- c(beta=1.0, gamma=0.3)
state <- c(S=0.999, I=0.001, R=0)

SIRmodel <- function(t, state, parameters) {
  with(as.list(c(state, parameters)),{
    dS <- - beta*S*I
    dI <- beta*S*I - gamma*I
    dR <- gamma*I
    list(c(dS, dI, dR))
  })
}

times <- seq(0, 50, by=0.2)

ode_out <- ode(y=state, times=times, func=SIRmodel, parms=parameters)

###ggplot2###
library(tidyr)
library(dplyr)
library(ggplot2)
library(scales)
ode_out2 <-as.data.frame(ode_out) %>% 
  gather(variable, value, -time)

theme_set(theme_bw(20))

ggplot(ode_out2)+ 
  geom_line(aes(x=time,y=value,colour=variable),size=1.5)+
  scale_colour_manual(values = c("royalblue","orange2","forestgreen" ),name="")

f:id:abrahamcow:20150121043711p:plain

###matplot###
matplot(ode_out[,1],ode_out[,2:4],type="l",lty=1,lwd=2,
        col=c("royalblue","orange2","forestgreen"),
        xlab="time",ylab="value")
legend(locator(1),legend=c("S","I","R"),lty=1,lwd=2,
       col=c("royalblue","orange2","forestgreen"))

f:id:abrahamcow:20150217113412p:plain