Susceptible-Infected-Recovered (SIR) モデル
Susceptible-Infected-Recovered (SIR) モデルは、 代表的な感染症モデルの１つで、社会に伝染病が持ち込まれた時に、それがどのように伝播し、やがて収束するかを記述します。SIRモデルでは人口は3つのグループに分かれると仮定します。1
- 感染可能者（Susceptible）：これから感染する可能性のある人
- 感染者（Infected）：現在感染している人
- 除外者（Recovered）：死亡・回復・隔離された人
それぞれの人口における割合を $x, y, z$ とすると、SIRモデルは次の微分方程式で表されます。
\begin{align}
\dot{x} &= - \beta x y \\
\dot{y} &= \beta x y - \gamma y \\
\dot{z} &= \gamma y
\end{align}
ここで、$\dot{x}$ は変数の時間変化（$\mathrm{d}x/\mathrm{d}t$）を表します。
$\beta$ は感染率を意味するパラメータで、SとIが接触した際に病気が移る割合と解釈できます。単位時間で新たに感染する人数は、SとIの接触頻度に依存することから、$x, y$ に比例すると仮定します。
$\gamma$ は、Iが回復・死亡・隔離などの理由によって、単位時間で感染者でなくなる割合です。病気が治りやすい場合 $\gamma$ は高くなる一方、非常に深刻で感染者がすぐに死に至る場合にもやはり $\gamma$ は高くなります。
人は S→I→R とグループを推移していき、Rになった人はもう感染しないし、他人を感染させもしない（免疫を保持）と仮定しています。ただし、全員が感染するとは限らず、Sのまま、Iの人数がゼロになり収束を迎えることもあり得ます。
なお、 $\dot{x} + \dot{y} + \dot{z} = 0$ であるため、社会全体の人口は一定です。これは比較的短期間の減少と考え人口動態を無視できるという仮定です。
この記事では、SIRモデルの基本的な挙動をシミュレーションと少しの数式で理解し、モデルをもとにコロナウィルス対策について考えます。
SIRモデルの挙動
数式を考える前に、まずはシミュレーションで感触を得てみます。可視化の便利さに鑑みてRを使います。
library(dplyr)
library(tidyr)
library(ggplot2)
library(gridExtra)
sir_simulate <- function(beta, gamma, maxt=150, y0=1e-6, z0=0) {
x0 <- 1 - y0 - z0
x <- numeric(1+maxt)
y <- numeric(1+maxt)
z <- numeric(1+maxt)
x[1] <- x0
y[1] <- y0
z[1] <- z0
for (i in 1:maxt) {
x[i+1] <- x[i] - beta*x[i]*y[i]
y[i+1] <- y[i] + beta*x[i]*y[i] - gamma*y[i]
z[i+1] <- z[i] + gamma*y[i]
}
data.frame(t=0:maxt, susceptible=x, infected=y, recovered=z)
}
# beta=0.3, gamma=0.1 でシミュレーション
simdata <- sir_simulate(0.3, 0.1)
pivot_longer(simdata, cols=-t,
names_to="type", values_to="fraction") %>%
ggplot(aes(t, fraction, color=type, linetype=type)) +
geom_path() +
theme_bw()
時間を通じて、S は減少、Rは増加を続けます。I（感染者）の人数は当初は増加していき、いずれ頭打ちになり減少に転じ、やがて収束します。
感染者を減らすためには
感染者が減少するための条件については、もとの微分方程式から考えるとわかりやすい。
$\dot{y} = \beta x y - \gamma y $ だから、
\dot{y} < 0
\Longleftrightarrow \beta x y - \gamma y < 0
\Longleftrightarrow x < \frac{\gamma}{\beta}, y \not= 0 \tag{1}
パラメータ $\beta, \gamma$ を固定すると、$x$（感染可能者）が十分に少なくなるしか、感染者を減らす方法はないことになります。つまり多くの人が感染すれば、新たな感染者が出にくくなり、いずれ収束に向かうという理屈です。これがいわゆる「集団免疫」という考え方です。
政策は、パラメータ $\beta, \gamma$ を外的に変化させる試みと解釈できます。
上の条件式 (1)から、感染者を減少させる方法が２つ考えられます。
- $\gamma$ を高くする：感染者の回復を早める、隔離を徹底する、薬ができる
- $\beta$ を低く抑える：人と人の接触を避けるよう外出規制などを課す、ワクチンができる
1つ目の案が可能なら理想的ですが、現実には難しいかもしれません。特効薬は現時点では存在せず、回復が治まるまで１〜２週間かかるのは避けられないのが現状のようです。またそもそも全ての感染者を確認できているわけではないので、隔離の強化にも限界があります。
現実にしばしば考えられているのは、$\beta$ を低く抑える2つ目のアプローチです。ワクチンの開発は理想的ですがいつ実現するかは不明です。そこで、人同士の接触を避けるため、外出の自粛を促したり、罰則をつけて取り締まったりする政策が講じられます。
これは、根本的というより一時的な緩和措置という側面が強いようです。
たとえば、現在感染者が増加している局面とします。つまり、$x > \frac{\gamma}{\beta}$ が成り立っている状況です。ここで一時的に $\beta$ を十分低くして感染者を減少させます。一方、$\beta$ が低いためこの間に新たな感染者があまり増えず、したがって $x$ はあまり変化しません。そのため、政策を止めると再び上の条件が復活する可能性が高いのです。このとき、感染者がわずかでも残っていると、そこから再び感染が蔓延してしまいます。
SIRモデルのもとでは、根本的な解決には集団免疫の獲得（$x$ が低くなり感染者が減る）か新薬・ワクチンの開発（$\beta, \gamma$ が恒常的に変化し、感染者が減る）が必要です。
この点は位相図で考えると直感的です。位相図とは、２つの変数について、２次元グラフ上の各点で変数がどのように推移していくかを可視化したものです。
sir_phase_diagram <- function(beta, gamma, ngrid=21, ...) {
expand.grid(x=seq(0,1,length=ngrid), y=seq(0,1,length=ngrid)) %>%
mutate(delta_x=-beta*x*y, delta_y=beta*x*y-gamma*y) %>%
ggplot(aes(x, y, xend=x+delta_x, yend=y+delta_y)) +
geom_segment(...) +
xlab("Susceptible") + ylab("Infected") +
coord_cartesian(xlim=c(0,1), ylim=c(0,1)) +
theme_bw()
}
sir_phase_diagram(0.3, 0.1, arrow=arrow(length = unit(0.1,"cm")))
上の図は、横軸にS（感染可能者）、縦軸にI（感染者）の割合をとり、各変数の推移を矢印で示したものです。当初右下から感染拡大していき、ある時点で接触規制を行い $\beta$ を外生的に下げます。すると、感染者は減少し、一方あまり感染可能者はあまり減りません。規制を終了すると、再び感染が始まります。
とはいえ、この図の規制は失敗例とは言えません。感染は再発するものの、規制の結果、ピーク時の感染者数を抑えることができています。いっときに感染者の数が拡大しすぎると医療機関が機能不全に陥るいわゆる「医療崩壊」の危険があるため、感染者の数を抑えながら集団免疫の獲得していくことに寄与していると見ることができます。
外出規制のタイミング
目下、外出禁止措置を行うか否かが議論になっています。これはモデルにおける $\beta$ を一時的に下げようという強い政策です。この政策はタイミングが難しく、早すぎる規制は集団免疫の獲得を遅らせるだけになってしまい、逆に遅すぎれば感染者数が拡大し医療崩壊の恐れがあります。規制のタイミングでその後の経過にどのような違いが出るか、シミュレーションで検証してみます。2
simulate_temporary_policy <- function(beta, gamma, t_policy, beta_policy,
maxt=150, y0=1e-6, z0=0) {
x0 <- 1 - y0 - z0
x <- numeric(1+maxt)
y <- numeric(1+maxt)
z <- numeric(1+maxt)
x[1] <- x0
y[1] <- y0
z[1] <- z0
for (i in 1:maxt) {
beta_ <- if (i %in% t_policy) beta_policy else beta
x[i+1] <- x[i] - beta_*x[i]*y[i]
y[i+1] <- y[i] + beta_*x[i]*y[i] - gamma*y[i]
z[i+1] <- z[i] + gamma*y[i]
}
data.frame(t=0:maxt, susceptible=x, infected=y, recovered=z)
}
simdata_policy <- NULL
for (pstart in seq(25, 80, by=5)) {
tmp <- simulate_temporary_policy(0.3, 0.1, pstart:(pstart+19), 0.15) %>%
mutate(policy_period=sprintf("%d-%d", pstart, pstart+19))
simdata_policy <- rbind(simdata_policy, tmp)
}
ggplot(simdata_policy, aes(t, infected)) +
facet_wrap(vars(policy_period), ncol=3) +
geom_line(linetype="dotted") +
geom_line(aes(t, infected), data=simdata) +
theme_bw()
上の図は、$\beta = 0.3, \gamma = 0.1$ におけるシミュレーションをベースラインとして(実線)、途中20日間だけ $\beta = 0.15$ に抑えた場合（点線）と比較しています。各パネルは異なる政策実施タイミングを示しています。
結果を見ると早すぎる介入はあまり役立たないことが分かります。開始時期が50期以前だと、感染者数は単に右へシフトしただけで、ピークの高さはあまり変わりません。時間稼ぎにはなるものの、医療崩壊のリスクを軽減させる効果は見られません。
開始時期55期目あたりから徐々にピークの感染者数を下げる効果が見られます。特に65〜75期目に開始すると、感染者のピークを抑えながら、減少に向かわせることに成功しています。
また、80期目のように介入が遅すぎるとほとんど効果がありません。
この結果からは、感染者数が急速に拡大したあたり、ちょうど変曲点を過ぎたあたりで介入すると、効果が大きいように見られます。変曲点は、増加している局面に置いては2階微分がゼロ、1階微分が極大となるので、感染者数の日次変化が下がり始めたタイミングということになります。
下記コードでは、感染者数の1階・2階微分（差分で近似したもの）と変曲点を可視化しています。シミュレーション結果における変曲点（1階微分が最大となる点）は、70期目でした。
tmp <- simdata %>%
mutate(d_infected=(lead(infected) - lag(infected))/2,
d2_infected=(lead(d_infected) - lag(d_infected))/2) %>%
select(t, infected, d_infected, d2_infected) %>%
pivot_longer(-t, names_to="variable", values_to="value") %>%
mutate(variable=factor(variable,
levels=c("infected", "d_infected", "d2_infected")))
inflection <- filter(tmp, variable=="d_infected") %>%
arrange(desc(value)) %>%
head(1) %>%
select(t) %>%
as.numeric()
ggplot(tmp, aes(t, value)) +
geom_line() +
geom_vline(xintercept=inflection, size=0.5, linetype="dotted") +
facet_wrap(vars(variable), ncol=1, scale="free_y") +
ylab(element_blank()) +
theme_bw()
日本・東京のいま
ここまでの検討を踏まえて日本の感染者数推移を見てみます。データは、有志がデータを集約してくれている、2019-ncov-japan というGitHubレポジトリを利用します。
jpn <- read.csv("https://raw.githubusercontent.com/swsoyee/2019-ncov-japan/master/Data/detailByRegion.csv",
stringsAsFactors=FALSE) %>%
setNames(c("date", "pref",
"confirmed", "hospitalized", "recovered", "dead")) %>%
mutate(date=as.Date(as.character(date), format="%Y%m%d"))
head(jpn)
データは、3月19日から日別・都道府県別に「感染確認者」「入院中」「退院者」「死亡者」の人数を捕捉しています。入院中以外は累積数です。
感染確認者から退院・死亡者を除いた人数をSIRモデルにおける「感染者」と考え、その時系列推移と1階・2階微分を調べてみます。
g1 <- mutate(jpn, infected=confirmed-recovered-dead) %>%
group_by(date) %>%
summarize(infected=sum(infected)) %>%
ungroup() %>%
mutate(pref="Total") %>%
select(date, pref, infected) %>%
mutate(d_infected=(lead(infected) - lag(infected))/2,
d2_infected=(lead(d_infected) - lag(d_infected))/2) %>%
pivot_longer(-c(pref, date), names_to="variable", values_to="value") %>%
mutate(variable=factor(variable, levels=c("infected", "d_infected", "d2_infected"))) %>%
ggplot(aes(date, value)) +
geom_point() +
geom_smooth(size=0.5, alpha=0.1) +
facet_grid(variable ~ pref, scales="free_y") +
ylab(element_blank()) +
theme_bw()
library(gridExtra)
g2 <- mutate(jpn, infected=confirmed-recovered-dead) %>%
filter(pref=="東京都") %>%
mutate(pref="Tokyo") %>%
select(date, pref, infected) %>%
mutate(d_infected=(lead(infected) - lag(infected))/2,
d2_infected=(lead(d_infected) - lag(d_infected))/2) %>%
pivot_longer(-c(pref, date), names_to="variable", values_to="value") %>%
mutate(variable=factor(variable, levels=c("infected", "d_infected", "d2_infected"))) %>%
ggplot(aes(date, value)) +
geom_point() +
geom_smooth(size=0.5, alpha=0.1) +
facet_grid(variable ~ pref, scales="free_y") +
ylab(element_blank()) +
theme_bw()
grid.arrange(g1, g2, nrow=1)
このグラフを見ると、3月末あたりで変曲点（2階微分がゼロ）を越えているように見えます。すると、そろそろ強い外出規制が高い効果を期待できるタイミングなのかもしれません。
ただし、これは3月27日頃の自粛要請の影響で感染者の増加が緩やかになった可能性もあるので、実際にはもう少し先なのかもしれません。
ただし、実際の規制実施のタイミングにおいては医療機関のキャパシティや経済効果をはじめたくさんの考慮要素があることは言うまでもありません。仮にもう少し後に実施したほうが高い政策効果が期待できるとしても、医療現場の限界が近いのであればそちらを優先して早期に規制を強化することも考えられます。
-
Burghes & Borrie (1981) Modelling with Differential Equations（垣田・大町訳 (1990)『微分方程式で数学モデルを作ろう』）の7.4節を参照 ↩
-
このシミュレーションは https://sites.google.com/site/aatoda111/misc/covid19 を参考に、簡素化したものです。 ↩