SIRモデルからはじめる微分方程式と離散時間確率過程（前編）

今年はデング熱やエボラで騒がれました。そのような感染症の伝播によって感染人数がどのように変化するかを表すモデルはいくつかありますが、最もシンプルなものはSIRモデルというものです。Wikipediaの記事はこちら。

総人口をNとして、Sが感受性人口（まだ感染してないけど感染する可能性がある人）、Iが発症感染者人口、Rが除外人口（Iから治ってもう感染しない人）を表します。現時点のS,I,Rが与えられた時に、○○時間経った後のS,I,Rはどのようになるでしょうか？○○時間が短い極限をとって連続時間で変化のルールを記述したものを微分方程式と呼びます。SIRモデルでは次のような変化のルールに従うと仮定します。

上記のルールは化学反応式とほとんど同じです。SとIが接触したらある反応速度係数に従ってSもIになる。接触頻度はS*Iに比例するはず。I自体も単体で、ある反応速度係数に従ってRになる。これらを式に書き下して連立微分方程式にすると以下のようになります。

1行目の読み方は「Sがある短い時間に変化する量は-betaSI」です。その他も同様です。微分方程式の挙動を知るにはいくつか方法がありまして、一番よいのは紙と鉛筆で厳密解を求めること。しかし、微分方程式がちょっとややこしくなると厳密解を求めるのはすぐに不可能になり、コンピュータで初期値を適当に与えて数値計算するしかなくなります。それでも紙と鉛筆は活躍します。例えば連立微分方程式が2変数の場合、その2変数を横軸縦軸にしてある初期値から時間とともにどんな感じで動いていくかが知りたいことになります。その平面を相平面（phase plane; 相空間 phase spaceという単語もあります）と呼びます。相平面での挙動を探るには、微分方程式の時間変化の途中（軌跡やtrajectoryと呼ばれます）はさておいて、最終的な落ち着き先（アトラクター;attractorと呼ばれます）をまずは知りたいです。それは時間変化がなくなるところ、すなわち連立微分方程式の右辺が全部ゼロになるところがアトラクターになる可能性を持っていますので求めます。連立微分方程式の各式の右辺がゼロと置いて出てきた境界線みたいのをヌルクライン（nullcline）と呼びます。式が2つならヌルクラインは2つ以上になり、その交点がアトラクターになる可能性を持っています。その交点を求めたあとは、その交点近傍で線形安定性を求めることになります。線形化して行列の固有値を求めて、その符号を見て、全部マイナスなら安定な不動点でアトラクターになります。その近傍の点は時間とともにそこへ吸い込まれます。その次にヌルクラインの間がどんな流れになるのか式の符号を求めたり、数値計算で初期値をいくつか振って見てみます。そしてvector fieldなんか描いたりします。ここまででだいたい微分方程式の挙動を知ることができます。

さらに式が少し複雑になるとアトラクターが点ではなく周期的に振動したりします。ヌルクラインの交点だけでは分かりません。そのようなアトラクターはリミットサイクル（limit cycle）と呼ばれます。またアトラクターが複数になって、初期値の場所によって行き着き先が変わることもあります。そのようなあるアトラクターに吸い込まれる初期値の領域をそのアトラクターのベイシン（basin）と呼びます。英語で洗面器の意味です。アトラクターからゆらぎで少しはみ出ても、ベイシンの中ならすぐに元のアトラクターに吸い込まれます。また、線形安定性を見た時に、固有値が正負混じっていると鞍点になります。ある方向からは吸い寄せるけど、ある方向へは遠ざかっていくような点になります。そのような点や安定点がたくさんできるようなモデルだと、かなり複雑な挙動になります。また、変数の数が増えて次元数が上がると、なかなかヌルクライン描いて…というわけにもいかなくなり、初期値をたくさん振るのも現実的な時間では不可能になってきます。

それでも微分方程式は解の挙動を定性的に把握できるメリットがあります。現象の安定性をアトラクターとベイシンで簡潔に記述できることにメリットがあります。時間変化のルールに既知の知識を組み込めるメリットがあります。ゲーム理論なんかでは「みんなが同じその最適戦略をとったらどうなるのか？」という問題がよく出ますが、そのような際には微分方程式でモデル化することで、あるアトラクターに落ち着くのか、発散するのかというようなことが分かります。また、最近の有名になった例ですと、ヒット現象の数理モデル（slideshare）なんかにも使われています。

前置きが長くなりましたが、今回のSIRモデルではRの{deSolve}パッケージを使って数値的に解きます。また、SIRモデルでは人数の挙動なので本来は離散的な値ですが、ここではNが十分に大きいとしてS,I,Rを各々Nで割って、化学反応の「濃度」のような連続的な値を考えます。するとN*betaをbetaに置き換え直すだけで元の式の形が変わりません。ただし人数が少ない時の効果がクリティカルに効いてくる場合は濃度のように置き換え直すのがNGになり、微分方程式はダメで、人数の確率過程で正確に記述しないといけません。また定性的というより定量的な把握が重要な場合も確率過程を使うことになります。確率過程バージョンは次回の記事で触れる予定です。

さて今からbeta=1.0、gamma=0.3の場合の挙動を見ていきます。Rの値はあまり関係ないのでSとIの相平面に注目します。元の意味からSとIが0以上だけを見ることにします。ヌルクラインは、Sの時間発展の式からS=0 or I=0、Iの時間発展の式からはS=gamma/beta or I=0、Rの時間発展の式からはI=0です。それらの交点は今回は線となってしまい、線形安定性解析の例題としては向いていませんね…。流れのvector fieldを描くと以下のようになりました。

黒線はヌルクラインです。ヌルクラインをまたぐと微分方程式の右辺の符号が変わることを意味するので、流れの向きが変わります。これを見てもS<gamma/betaの範囲のI=0は安定なアトラクターだと分かります。ベイシンはSとIが0以上の全域です。どんな初期値でも最後は上記のアトラクターに落ち着くと分かります。パンデミックの視点からは、ワクチンの意図は、十分にSを減らしておくことでIが0からちょっと増えたとしても、すぐにI=0へ吸い込めるようにしておくことだと考えることができます。今回は初期値がS=0.999, I=0.001, R=0の場合を考えます。vector fieldの図では右下のところから出発です。はじめは遅く、途中から速度をあげて真ん中のS=gamma/betaあたりに来ます。そして、ヌルクラインをまたいだ後は速度を落としながら左下のI=0に向かっていくことが分かります。その挙動は横軸に時間・縦軸に値をとると以下のような図になります。

さて実際には病院に来る感染者数（I）だけが観測できると思いますので、そこに今回は簡単のため正規分布に従う観測誤差を加えてデータを作ってみます。元々の意味を考えるとマイナスになっちゃいけないのですが、アホな観測機器だと思って許してください。詳しい話は次回の記事で扱います。

作られたI_obsのデータ（SIRmodel_data.txt）は以下のようになっています。

time	S	I	R	I_obs
0.0	0.9990	0.0010	0.0000	-0.0111
0.2	0.9988	0.0012	0.0001	0.0039
0.4	0.9985	0.0013	0.0001	0.0122
…	…	…	…	…
49.6	0.0408	0.0000	0.9591	0.0027
49.8	0.0408	0.0000	0.9591	0.0027
50.0	0.0408	0.0000	0.9591	0.0044

このデータをもとに、RからStanを呼び出すパッケージである{rstan}を使って、Stan2.5以降の新機能であるintegrate_ode()で微分方程式を解きつつ、観測誤差の大きさを推定します。Stanコードは以下になりました。

functions {
   real[] sir(
      real t,
      real[] y,
      real[] theta,
      real[] x_r,
      int[] x_i
   ) {
      real dydt[3];
      dydt[1] <- - theta[1]*y[1]*y[2];
      dydt[2] <- theta[1]*y[1]*y[2] - theta[2]*y[2];
      dydt[3] <- theta[2]*y[2];
      return dydt;
   }
}

data {
   int<lower=1> T;
   real I_obs[T];
   real T0;
   real TS[T];
   real Y0[3];
}

transformed data {
   real x_r[0];
   int x_i[0];
}

parameters {
   real<lower=0, upper=10> beta;
   real<lower=0, upper=10> gamma;
   real<lower=0> sigma;
}

transformed parameters {
   real y_hat[T,3];
   real theta[2];
   theta[1] <- beta;
   theta[2] <- gamma;
   y_hat <- integrate_ode(sir, Y0, T0, TS, theta, x_r, x_i);
}

model {
   sigma ~ cauchy(0, 2.5);
   for (t in 1:T)
      I_obs[t] ~ normal(y_hat[t,2], sigma);
}

• 37,41行目: integrate_ode()についてはStanのマニュアルを読んで下さい。y_hatはdeterministic（決定論的）に決まるので値が後から使いたい場合はtransformed parametersブロックに書くことになります。
• 22行目: 今回は初期値をdataとして与えていますが、マニュアルの例のようにparameterとして推定することも可能です。ですがモデルによっては非常に時間がかかるようになるようで今後に期待です。

キックするRコードは以下になりました。計算時間は1chainあたり15秒程度でした。

library(rstan)

df <- read.delim("input/SIRmodel_data.txt", sep="\t")
T <- length(df$time)-1
data <- list(
   T=T,
   I_obs=df$I_obs[-1],
   Y0=c(0.999, 0.001, 0),
   T0=df$time[1],
   TS=df$time[-1]
)

init_fun <- function() {
   list(beta=runif(1,0.5,1), gamma=runif(1,0,0.5), sigma=runif(1,0,0.01))
}

model_fit <- stan(file='model/model.stan', chains=0)
fit <- stan(
   fit=model_fit, data=data, init=init_fun, pars=c('y_hat', 'beta', 'gamma', 'sigma'),
   iter=1000, warmup=400, thin=1, chains=3,
   seed=1234
)

save.image("output/result.RData")

・4～11行目: integrate_ode()は現状では初期値の時間はTSに含めてはいけないようなので注意です。

結果は以下になります。 Iの推定結果は以下の図になります。黒線はデータであるI_obs、青線は真のI、赤線は推定されたIのMCMCサンプルの中央値、帯は95%信用区間ですが、細すぎてよく見えないです。

beta, gamma, sigmaの事後分布は以下のとおりです。黒点は真の値、カラーの点はMCMCサンプルの中央値、線は95%信用区間ですが、やっぱり見えない…。violin plotも見えない…。バッチリ求まりすぎです。

---

最後に上記の他に使用したRコードを載せます。 初期値をS=0.999, I=0.001, R=0にしてSIRモデルを{deSolve}で解いた例です。{deSolve}は速いし使い勝手がいいです。

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)
set.seed(1234)
ode_out <- transform(ode_out, I_obs=rnorm(length(times), mean=I, sd=0.01))
write.table(ode_out, file='input_ode.txt', sep='\t', quote=F, row.names=F)

• 19行目: Iに観測誤差を付け加えてI_obsを作っています。

vector fieldを描いたRコードは以下になります。

library(doParallel)
library(foreach)
library(ggplot2)
library(deSolve)
library(scales) # for muted

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))
   })
}

exec_ode_from_states <- function(df_state=df_state, times=times, func=func, parms=parms){
   cl <- makeCluster(detectCores()-1)
   registerDoParallel(cl)
   on.exit(stopCluster(cl))
   res_df <- foreach(i=1:nrow(df_state), .packages='deSolve', .combine=rbind) %dopar% {
      ode_out <- ode(y=c(unlist(df_state[i,])), times=times, func=func, parms=parms)
      data.frame(i.state=i, ode_out)
   }
   return(res_df)
}

times <- seq(0, 0.5, by=0.02)
beta <- 1.0
gamma <- 0.3
parameters <- c(beta=beta, gamma=gamma)
df_state <- expand.grid(S=seq(.03, 1, 0.05), I=seq(.03, 1, 0.05))
df_state$R <- 0
res <- exec_ode_from_states(df_state=df_state, times=times, func=SIRmodel, parms=parameters)

p <- ggplot(data=res, aes(x=S, y=I, group=i.state, color=time))
p <- p + geom_line(size=2) + geom_point()
p <- p + geom_vline(xintercept=0) + geom_vline(xintercept=gamma/beta) + geom_hline(yintercept=0)
p <- p + scale_colour_gradient2(low="white", high=muted("red"))
p <- p + coord_cartesian(xlim=c(-0.01,1), ylim=c(-0.01,1))
ggsave(plot=p, file="output/vector_field.png", dpi=300, width=6, height=5)

• 31-33行目: 初期値をexpand.grid()で作ってその各々を初期値としてSIRモデルを{deSolve}で少しの時間だけ解いています。
• 19行目: こちらの記事参照。@kos59125さん @hoxo_mさんThanks!

---

微分方程式に関する個人的にオススメの書籍を挙げておきます。

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

Hirsch・Smale・Devaney 力学系入門 ―微分方程式からカオスまで―

非線形力学の展望I カオスとゆらぎ

非線形力学の展望II 複雑さと構造

「微分方程式で数学モデルを作ろう」は易しくて読みやすいのですが、途中から出てくるモデルが厳密に解けるモデルばかりでやや実践的ではありません。Hirsch・Smaleの本は例も証明もしっかりしていて読みやすいです。残り二冊は証明はあまりなくて、あってもざっくりしていますが、網羅的で図がいっぱいで楽しいです。難度はやや高めです。いずれにせよ、自分でモデル式をたてて解いたりするのが必須で、それで理解がすすむと思います。