SIRモデルとロジスティック方程式の類似点　どれだけ似ていて何が違うのか

　日本では新規感染者数も収まりつつある新型コロナウイルス感染症(COVID-19)。累積感染者数のグラフはいよいよS字曲線になってきた。

f:id:argos_m1111:20200608042500j:plain — 6月7日までの日本の累積感染者数の推移(ジョン・ホプキンス大のCOVID-19 Dashboardより)

　前の記事でも触れたが、ウイルスという生物的なふるまいをする物体の増殖曲線でS字といえば、大学レベルの生物学を学んだことがあれば真っ先に思い出すのが「ロジスティック方程式」(ロジスティック曲線)だ。

$\displaystyle \frac{dN}{dt} = r(1-\frac{N}{K})N$ ……ロジスティック方程式

( $0 \lt N \lt K$ , $r \gt 0$ )

　 $N$ はある時点 $t$ における個体数だ。ロジスティック方程式は分かりやすい特徴を持っている。第一に、環境収容力( $K$ )と呼ばれる定数で個体数の上限が決まること。第二に、増殖速度が緩やかになり始める変曲点が $N=K/2$ であること。つまり理論的には、変曲点が分かれば環境収容力が大体どの辺なのか推測できる。

f:id:argos_m1111:20200608042748j:plain — ロジスティック方程式(曲線)の例

　しかし、実際の生物の増殖曲線がロジスティック方程式に従うのは限定的といわれている。理由はいろいろあるが、ざっくり言えば、現実に増殖速度に影響を及ぼす要因はもっと複雑だからだ。

　そして、ウイルス由来のものを含めた感染症の増殖曲線については、ロジスティック方程式ではなくSIRモデルといわれる感染症数理モデルをベースとするのが一般的だ。

　SIRモデルはロジスティック方程式のように1つの微分方程式では表さず、3つの式で表される。

$\displaystyle \frac{dS}{dt} = -βIS$

$\displaystyle \frac{dI}{dt} = βIS - γI$

$\displaystyle \frac{dR}{dt} = γI$

( $0 \le S \lt K$ , $I \ge 0$ , $0 \le R \lt K$ , $S + I + R = K$ )

　 $S$ は感受性人口(感染症にかかりうる個体数)、 $I$ は感染者数(感染症にかかっていて、他人に感染させうる個体数)、 $R$ は回復数(感染から回復、あるいは隔離され、もう感染しない個体数)を表す。ロジスティック方程式の $N$ と同じように、どれもある時点 $t$ における個体数だ。注目している全人口 $K$ の各人がいずれかに属するので、全て足し合わせると $K$ になるのはすぐ分かる。

f:id:argos_m1111:20200608043616j:plain — SIRモデルの各曲線。Sはまだ感染していない人の数、Iは感染していて他人に感染させることができる人の数、Rは隔離・回復してもう感染しない人の数

　ロジスティック方程式もSIRモデルも微分方程式という形を取っており、個体数と増殖速度の関係式ともいえる。つまり、個体数と時間の関係式ではないため、S字曲線自体を表す式にはなっていない。

　ロジスティック方程式ではこれを解き、個体数と時間の関係式に直すことでS字曲線を表す関数を得られるのだが、SIRモデルでそれができるかというと自明ではない。少し調べた感じだと、SIRモデルは厳密には解けず、近似的に描かれることが一般的なようだった。

　しかし、SIRモデルの近似解による曲線を見ると、やはりどうもロジスティック曲線と似ているように見える。この2つのモデルはどこまで似ていて、どこが違うのだろうか。私はこれが気になった。

　この記事では、私が調べた限りでの、SIRモデルとロジスティック方程式の類似点や違いを、数式をベースに紹介したい。これを調べ始めたのは4月下旬のことで、GWの時間のほとんどを数式やグラフの比較に費やしてしまった。正直今も考えている最中ではあるのだが、キリがないので一つの結論が見えたところで記事としてまとめることにした。

　ただ、調べたと言っても文献にはほとんど当たっていない。あくまで、ずっと自分で式とグラフをいじくり回していただけだ。もしかしたらこの内容は感染症の教科書には載っていることなのかもしれない(少なくとも車輪の再発明ではあるだろう)が、自分の一連の思考をまとめるという意味の記事だと思ってほしい。

　レベル感としては学部の数理生物学のレポート1本分くらいだろうか。一応、説明や式変形はできるだけ丁寧に行うつもりだ。理系の高校生くらいなら分かってくれると思う。

　長い記事になるので、端的な結論だけまず述べる。

「SIRモデルのある係数を0と置いたとき、

SIRモデルはロジスティック方程式と一致する」

　分かる人はもう分かったかもしれない。分からない人はもやもやしながら記事を読み進んですっきりしてほしい。

　では、まず解説の順番を鳥瞰しよう。

1. ロジスティック方程式の解法を知る

1-1. 微分方程式からロジスティック方程式の特徴を得る

1-2. 式を正規化する

1-3. 微分方程式を解き、解(関数)を得る

2. SIRモデルにロジスティック方程式の手法がどこまで使えるのか検討する

2-0. 表計算ソフトで近似曲線を書いてみる

2-1. 微分方程式を整理してSとIの関係式を得る

2-2. 変数変換と正規化でロジスティック方程式と比べやすくする

2-3. SIRモデルは本当に解けないのか?

2-4. 微分方程式の形でロジスティック方程式と比較してみて分かること

　ロジスティック方程式の解法や特徴の調べ方を知り、それらの手法をSIRモデルに適用できるのかどうか検討していく。このような順番で説明していこう。ロジスティック方程式を熟知している人は前半は読み飛ばして構わないが、ちょっと不安、もしくはよく知らないという人は目を通しておいた方が全体を通して分かりやすくなると思う。

1. ロジスティック方程式の解法を知る

1-1. 微分方程式からロジスティック方程式の特徴を得る

　まず、ロジスティック方程式を初めより少し簡単な形で表す。

$\displaystyle \frac{dN}{dt} = mN$

ただし、 $\displaystyle m = r(1 - \frac{N}{K})$

( $0 \lt N \lt K$ , $r \gt 0$ )

　いきなりこの式を理解せよと言われても面食らうので、少し整理してみよう。

　 $N$ は個体数で、時間 $t$ の関数だ。左辺 $\dfrac{dN}{dt}$ は $N$ の1階導関数、つまり増殖曲線の傾き、つまり個体数の増殖速度を表している(進んだ距離を時間で割ると速度になるのと同じ考え方)。

　右辺は個体数 $N$ にある係数 $m$ を掛けている。つまり、個体数の増殖速度は、その時点の個体数 $N$ に $m$ を掛けたものという意味だ。

　 $m$ は本来正負を問わない実数でよいが、ここでは想像を簡単にするために正の実数としておく。

　さて、この式の形は微分方程式と呼ばれる。元の関数とそれを微分したものの関係式をうまく解けば、元の関数の形を得られる。例えば、元の関数 $y$ とそれを微分したもの $y'$ が等しいなら、元の関数が $±e^ x$ という指数関数になることは、高校の微分積分を知っていればすぐに分かるだろう。

　もしこの $m$ が定数なら、この微分方程式は単純な指数関数になる。しかしそうは問屋が卸さない。

　ロジスティック方程式では、「個体数の増加に比例して、増殖率が減る」という考え方を導入する。つまり、横軸に個体数、縦軸に増殖率を取ったグラフを描いたとき、直線で増殖率が減っていき、個体数の上限であるKで増殖率がゼロになるような線を描く方程式が必要だ。

　増殖率の初期値を $r$ とし、増殖率を $m$ とすると、この直線は

$\displaystyle m = - \frac{r}{K} N + r$

　となる。( $m$ を $y$ 、 $N$ を $x$ と読み換えれば、親しみのある式の形だと分かるはずだ)

　これを整理すると、

$\displaystyle m = r(1 - \frac{N}{K})$

　となるので、先程の式に代入すると、最初に提示した

$\displaystyle \frac{dN}{dt} = r(1-\frac{N}{K})N$

　が得られる。

　つまり、 $N = f(t)$ の具体的な $f(t)$ が分からなくても、この微分方程式を解釈するだけで $N = K$ で増殖速度はゼロになる(厳密には極限を取るとだが)ことは見えてくる。

　初めは指数関数的に個体数が増え、最終的にはKに漸近するのだから、どこかに変曲点があるはずだ。どこにあるのだろうか。

　関数の具体的な形が分かっていれば、2階導関数まで求めて増減表を書いてみてもいいが(学生時代にこの方法で求めた記憶がある)、今の説明ではまだ具体的な関数は求めていない。

　ロジスティック方程式の場合、変曲点は「増殖速度(曲線の傾き)が最大の場所」と考えられる。そしてこの微分方程式は傾きに関する関数とも取れる。ということはこの微分方程式を横軸に個体数、縦軸に傾きを取るグラフ上に当てはめて、最大値を求められるならそのときの個体数が変曲点の座標になりそうだ。

　右辺を式変形してみると、

$\displaystyle -\frac{r}{K} ((N - \frac{K}{2})^ 2 - \frac{K}{4})$

　という2次関数であることが分かるから、 $\displaystyle N = \frac{K}{2}$ で傾きは最大となる。これで冒頭で説明したロジスティック方程式の特徴が分かった。

1-2. 式を正規化する

　さて、もし実際に生物の増殖を調べるなら、 $K$ という個体数の上限を示す定数があるのはありがたいことだが、ここでは $0 \lt N \lt K$ の範囲で動く変数 $N$ を $0 \lt n \lt 1$ と変数変換することで正規化したい。そうすると解を求めるのが少し楽になるのと、SIRモデルとの比較もしやすくなるからだ。

　 $\displaystyle 0 \lt \frac{N}{K} \lt 1$ だから、 $\displaystyle \frac{N}{K} = n$ とおける。両辺を $t$ について微分すると $\displaystyle \frac{1}{K} \frac{dN}{dt} = \frac{dn}{dt}$ だから、これらを与式に代入すると、

$\displaystyle \frac{dn}{dt} = r(1-n)n$

　が得られる。最初より少し簡単な式になったような気がする。

1-3. 微分方程式を解き、解(関数)を得る

　では、この式を解いていこう。この微分方程式には「変数分離法」と呼ばれる解法を用いる。具体的には、両辺を $(1-n)n$ で割ることで左辺に $n$ を集め、両辺を $t$ について積分する。割るに当たっては、 $0 \lt (1-n)n \lt 1$ なので0割りにはならないことに一応注意しておきたい。

$\displaystyle \int \frac{1}{(1-n)n} \frac{dn}{dt} dt = \int r dt$

　すると、右辺は定数の積分なので

$\displaystyle \int r dt = rt + C_1$ $C_1$ は積分定数

　となる。一方の左辺はどうだろう。まず、置換積分の考え方から $dt$ を消去して、 $n$ で $\dfrac{1}{(1-n)n}$ を積分すると考えられる。

$\require{cancel} \displaystyle \int \frac{1}{(1-n)n} \frac{dn}{\cancel{dt}} \cancel{dt} = \int \frac{1}{(1-n)n} dn$

　そして、 $\dfrac{1}{(1-n)n}$ はなんとうまいことに部分分数分解により、 $\dfrac{1}{n} + \dfrac{1}{1-n}$ と変形できる。各項に対してそれぞれ積分を考えたらいいので、

$\displaystyle \int \dfrac{1}{x} dx = \log{|x|} + C$

　という分数関数の積分の公式から、

$\displaystyle \int (\dfrac{1}{n} + \dfrac{1}{1-n}) dn = \log{|n|} - \log{|1-n|} + C_2$

　と積分できる。

　 $n$ も $1-n$ も正なので絶対値は外せる。logはまとめてしまおう。後処理の関係で $1-n$ を分子に持ってくる。さきほど積分した右辺との式は、

$\displaystyle -\log{\frac{1-n}{n}} = rt + C$

　両辺に $-1$ をかけ、指数化すると、

$\displaystyle \frac{1-n}{n} = e^{-rt-C} = e^{-C} e^{-rt}$

　となる。

　ここで、 $e^{-C}$ は定数なので、簡単のために $A$ とおく。

　式を整理すると、以下のようになる。

$\displaystyle \frac{1}{n} = 1 + A e^{-rt}$

$\displaystyle ∴ n = \frac{1}{1 + A e^{-rt}}$

　ここで止めてもいいが、 $A$ が若干正体不明なので、個体数の初期値 $\dfrac{N_0}{K}$ と $t=0$ を式に代入して、

$\displaystyle \frac{N_0}{K} = \frac{1}{1 + A e^{0}}$

$\displaystyle ∴ A = \frac{K}{N_0} - 1$

　を得る。先程の式に戻して、最終的に

$\displaystyle n = \frac{1}{1 + (\frac{K}{N_0} - 1) e^{-rt}}$

　という解が得られた。これがS字曲線を表す「ロジスティック曲線」の関数だ。

　試しに、スプレッドシートでグラフを書いてみよう。

f:id:argos_m1111:20200608041934j:plain — ロジスティック曲線の各係数に適当な値を入れ、スプレッドシートで描画した。傾きはk項の前後k±1項を引いて近似している。K/2で確かに傾きが最大になっていることが分かる

2.SIRモデルにロジスティック方程式の手法がどこまで使えるのか検討する

　ここまでロジスティック方程式の解法を見てきた。ここからが本題だ。

2-0.表計算ソフトで近似曲線を書いてみる

　いきなり微分方程式の整理を始めてもいいが、その前にSIRモデルがどんな曲線を描くのか、近似することで確かめてみたい。

　だが、微分方程式を解いて曲線の関数を求めるのが難しそうなのにどうやって近似するのか。ここで、微分の定義を振り返ってみよう。

　微分して求めた導関数 $f'(x)$ に数 $a$ を代入した $f'(a)$ は、 $x = a$ における曲線 $y$ の傾きを表している。

　傾きは、 $x, y$ の増分を $⊿x, ⊿y$ と表すと、 $\dfrac{⊿y}{⊿x}$ だ。つまり、 $b \gt a$ として、 $\dfrac{f(b)-f(a)}{b-a}$ と書ける。

　この $b$ を限りなく $a$ に近づけたときに得られるのが $x = a$ における傾きであって、この操作こそが微分だった。

$\displaystyle \lim_{b\to a}\frac{f(b)-f(a)}{b-a} = f'(a)$

　 $b$ が $a$ より少し大きい数 $a+h$ と考えると、以下のようにも書ける。

$\displaystyle \lim_{h\to 0}\frac{f(a+h)-f(a)}{h} = f'(a)$

　では、あえて極限を取らないとどうなるか。例えば、 $x=t$ と $x=t+1$ の間で傾きを考えてみると、

$\displaystyle \frac{f(t+1)-f(t)}{(t+1) - t} = f(t+1)-f(t)$

　となる。 $1$ ではなく、もっと小さな隙間 $δ$ で考えると、

$\displaystyle \frac{f(t+δ)-f(t)}{(t+δ) - t} = \frac{f(t+δ)-f(t)}{δ}$

　 $f(t+1)$ や $f(t+δ)$ を $F_{t+1}$ 、 $f(t)$ を $F_{t}$ と見なすと、微分は数列 $F$ の漸化式に近似できる。この近似手法はオイラー法ともいわれる。

　例えば、 $\dfrac{dy}{dx} = y$ という微分方程式に適用すると、

$\displaystyle \frac{y_{n+1} - y_{n}}{δ} = y_{n}$

∴ $\displaystyle y_{n+1} = y_{n} + δy_{n} = (1 + δ)y_{n}$ という漸化式が得られる。

　数列の2項間漸化式が得られたのだから、あとは初期値の $F_0$ を与えるだけでコンピュータが次々に各項の近似値をはじき出してくれるという寸法だ。SIRモデルの3つの式を漸化式にして、近似的な曲線を書かせたのが以下だ。

　現実の感染症で現れるS字の曲線の代表的なものは累積感染者数なので、 $I$ と $R$ を足し合わせた上図の緑の線と考えられる。あとで変数変換するが、上記でも $I+R=N$ とおいた曲線も書いた。この微分方程式の右辺は、Sの右辺の符号を反転させたものと同じだ。

　SIRモデルの式を再掲する。

$\displaystyle \frac{dS}{dt} = -βIS$

$\displaystyle \frac{dI}{dt} = βIS - γI$

$\displaystyle \frac{dR}{dt} = γI$

( $0 \le S \lt K$ , $I \ge 0$ , $0 \le R \lt K$ , $S + I + R = K$ )

　 $\dfrac{dI}{dt}$ と $\dfrac{dR}{dt}$ を足すと、右辺は $βIS$ となり、 $S$ を反転させたものだと分かる。

2-1.微分方程式を整理してSとIの関係式を得る

　先ほど $I$ と $R$ を足してみたが、式をよく見れば分かるように $I$ は $R$ を内包している。ということは、 $S$ と $I$ の関係が重要そうだ。ここを調べてみたい。

$\displaystyle \frac{dI}{dt} = βIS - γI$

　は、以下のように整理できる。

$\displaystyle \frac{dI}{dt} = βIS(1 - \frac{γ}{βS})$

　 $\dfrac{dS}{dt}=-βIS$ だから、代入すると、

$\displaystyle \frac{dI}{dt} = \dfrac{dS}{dt}(\frac{γ}{βS} - 1)$

　両辺をtについて積分すると、

$\require{cancel} \displaystyle \int \frac{dI}{\cancel{dt}} \cancel{dt} = \int (\frac{γ}{β}\frac{1}{S} - 1)\dfrac{dS}{\cancel{dt}} \cancel{dt}$

∴ $\displaystyle I = \frac{γ}{β}\log{S} - S + C$

　を得る。

　ここで、積分定数をはっきりさせたいので、初期値 $S_0$ と $I_0$ を代入する。 $S_0 + I_0 = K$ であることに注意すると、

$\displaystyle C = K - \frac{γ}{β}\log{S_0}$

　を得る。

　元の式に代入すると、

$\displaystyle I = \frac{γ}{β}\log{S} - S + K - \frac{γ}{β}\log{S_0}$

∴ $\displaystyle I = \frac{γ}{β}\log{\frac{S}{S_0}} + K - S$

　という、 $S$ と $I$ の関係式を得られた。これを $\displaystyle \frac{dS}{dt} = -βIS$ に代入すると、

$\displaystyle \frac{dS}{dt} = - \{γ\log{\frac{S}{S_0}} + β(K - S)\}S$

　という $S$ の微分方程式が得られた。

2-2.変数変換と正規化でロジスティック方程式と比べやすくする

　このまま $S$ についての整理を進めていきたいところでもあるのだが、近似曲線でみたように $S$ は逆S字曲線だから、ロジスティック方程式と比べるには上下反転させたい。そして $S$ の曲線を上下反転させたのが $I+R$ の曲線だった。 $S+I+R=K$ なので、 $I+R=N$ , $S=K-N$ として変数を変換すると、 $S$ と $I$ の関係式は、

$\displaystyle I = \frac{γ}{β}\log{\frac{K - N}{S_0}} + N$

　という $N$ と $I$ の関係式になる。 $\dfrac{dS}{dt}= -\dfrac{dN}{dt}$ であること( $S=K-N$ より、両辺をtで微分すれば得られる)に気を付けながら $S$ の微分方程式にこれを代入すると、

$\displaystyle -\frac{dN}{dt} = -\{γ\log{\frac{K-N}{S_0}} + βN\}(K - N)$

∴ $\displaystyle \frac{dN}{dt} = \{γ\log{\frac{K-N}{S_0}} + βN\}(K - N)$

　を得る。式自体はやや複雑にはなったが、ロジスティック方程式と同じように $N$ と $t$ しかない微分方程式を得られた。

　ここで、ロジスティック方程式と同様に正規化を行う。 $\dfrac{N}{K}=n$ とするので、 $\dfrac{1}{K}\dfrac{dN}{dt} = \dfrac{dn}{dt}$ に気を付けながら変換していくが、ここでひとつごまかしを入れる。 $S$ (これから感染しうる人口数)の初期値 $S_0$ は $K$ (全人口)とほぼ等しいと考えられるので、 $S_0 = K$ としてみる。すると、

$\displaystyle \frac{dN}{dt} = \{γ\log{(1 - \frac{N}{K})} + βN\}(K - N)$

$\displaystyle K \frac{dn}{dt} = \{γ\log{(1 - n)} + βKn\}K(1 - n)$

∴ $\displaystyle \frac{dn}{dt} = \{γ\log{(1 - n)} + βKn\}(1 - n)$

　となる。この式をロジスティック方程式と比べようという考えだ。

2-3.SIRモデルは本当に解けないのか?

　実際に比べる前に、この微分方程式は解けないのか、ということを考えたい。実は、これについては私は検討しきれていない。

　ロジスティック方程式と同様に変数分離法で解きたいところなのだが、部分分数分解がうまくいかず、きれいに項を分けられない。

$\displaystyle \int \frac{1}{\{\log{(1 - n)} + \frac{βK}{γ}n\}(1 - n)} dn = \int γ dt$ としても、左辺をうまく積分する方法が見つからない……

　すると複雑な分数関数の積分になり、計算が困難になる。試しにwolfram alphaにも入れてみたが、有効そうな答えは出てこなかった。

　私は微分方程式についてはほとんど知らないため、今はマセマの常微分方程式の本を読んで手がかりを探している。

スバラシク実力がつくと評判の常微分方程式キャンパス・ゼミ

作者:敬之, 馬場
発売日: 2019/08/01
メディア: 単行本

2-4.微分方程式の形でロジスティック方程式と比較してみて分かること

　さて、あらためて正規化したロジスティック方程式とSIRモデルの方程式を並べてみよう。

$\displaystyle \frac{dn}{dt} = rn(1-n)$ ……ロジスティック方程式

$\displaystyle \frac{dn}{dt} = \{γ\log{(1 - n)} + βKn\}(1 - n)$ ……SIRモデル(累積感染者数)

　かなり似ているといえるのではないか。違うのは、SIRモデルには $\log$ の項があるということだ。この項には回復率 $γ$ が掛けられている。普通に考えて、感染症患者は回復するかもしくは死亡することで感染者数 $I$ からは除外されるため、 $γ$ が0になることはないが、もしも $γ=0$ になった場合、

$\displaystyle \frac{dn}{dt} = βKn(1 - n)$

　の式が残る。すると、ロジスティック方程式とSIRモデルで異なるのは係数 $r$ (※SIRモデルの $γ$ ではない)と $βK$ だけだ。

　 $r$ はロジスティック方程式の増加率。 $β$ はSIRモデルの感染率で、 $βK$ は感染者1人当たりが生み出す新規感染者数を示しているから、やはり増加率だ。

　つまり、γ=0のとき、SIRモデルの累積感染者数はロジスティック曲線となるといえる。

　それは具体的にどんな感染症なのか? 一度感染すると治ることはなく、死ぬこともなく、ずっと感染者を生み出し続ける、ゾンビのような病気かもしれない。

　そんな妄想は置いておいて、 $γ=0$ でロジスティック方程式と一致するのは実はもっと前の段階でも確かめられる(私は気付かなかったが)。

　 $\displaystyle \frac{dI}{dt} = βIS - γI$ で $γ=0$ として、 $S$ との関係を導くと $\dfrac{dS}{dt} = - \dfrac{dI}{dt}$ より $S=-I+C$ (Cは積分定数)だから、 $I=C-S$ 。積分定数 $C$ は $C = S_0 + I_0$ より $C = K$ 。