はじめに
要旨
- インフルエンザ
- コロナウイルス
について数理モデルを用いて分析してみました。
背景
今回は感染症モデルを数値計算することを試みましたが、筆者は感染症モデル等の数理モデルのプロフェッショナルではありませんのでご理解ください。大学・大学院では物理専攻でしたので少し計算機科学を触れたことがある程度です。本投稿の背景としては、今ホットなトピックとして挙がっているCoronavirusについてなにか定量的に分析できないか興味が湧いたからです。以前の感染症事例でSARSやMERSがありますが、感染症の収束時期の予想などのためにどのような数理モデルがあるのか調べてみるのは価値があるのではと思った次第です。
感染症モデル
SIRモデル
感染症の流行過程を記述するモデル方程式としてSIRモデルというものがあります。変数としてはSuscceptible(感染の疑いがある者)、Infected(感染者)、Recovered(免疫保持者)の3変数です。1927年の論文(“A Contribution to the Mathematical Theory of Epidemics")で発表されたモデルです。特に社会物理をはじめとした身近な現象を扱う分野でよく用いられることが多いです。
S+I \rightarrow 2I \\
I \rightarrow R \\
素過程は上の2過程となります。上は感染過程、下は回復過程です。化学反応式みたいなものだと思えば良いと思います。実際に微分方程式に落とし込むときには化学における反応速度論と同じ計算ロジックを使えば良いです。
数式
SIRモデルを定式化すると以下のようになります。
\frac{dS(t)}{dt} = -\beta S(t)I(t) \\
\frac{dI(t)}{dt} = \beta S(t)I(t)-\gamma I(t)\\
\frac{dR(t)}{dt} = \gamma I(t) \\
ある時刻$t$に対して、$S(t)$は感染の疑いがある人数、$I(t)$は感染者数、$R(t)$は免疫保持者数です。
- 上段の式➤感染に疑いがある人数は感染者との接触回数に比例して減ります
- 中段の式➤感染者数は接触回数に比例しますが治癒する人が増加すればその分減ります
- 下段の式➤回復する人数が増えれば感染者数が減ります
ここで$\beta$は感染率、$\gamma$は回復率や隔離率を表しています。また
-\beta I(t) \\
は感染力に相当します。
特になぜ上記の微分方程式のように記述できるのかは
http://www2.meijo-u.ac.jp/~tnagata/education/react/2019/react_03_slides.pdf
を参考にすると良いと思います。高校化学の反応速度の式を出すときと同じロジックです。
ここで上記の連立微分方程式を解くのですが私はこれを見て大学時代に学んだローレンツ方程式を思い出しました。ローレンツ方程式は物理学におけるカオスの領域で語られる非線形常微分方程式(Ordinary Differential Equation:通称ODE)です。かなりマニアックだと思うので詳細は語りませんが式だけ書いておきます。
\frac{dx}{dt} = -px+py \\
\frac{dy}{dt} = -xz+rx-y\\
\frac{dz}{dt} = xy-bz \\
昔大学の課題でこの数値シミュレーションをしたことがあったのでその時と同じようにシミュレーションすればいいはずと考えました。その時はルンゲクッタで解いた記憶がありますが今回は勉強のためにTensorFlowを使ってみようと思います。
TensorFlowとは
TensorFlowはGoogleが開発したオープンソースのMachine Learning向けのソフトウェアライブラリです。一言で言えば高速数値解析用のPythonライブラリですね。今回はMachine Learningはしませんが数値計算用として使用してみます。テンソルは聞きなれないと思いますが、行列を次元拡張したものだと私は理解しております。本投稿の趣旨とは外れるのでどこかで投稿しようと思います。TensorFlowの特徴としてはソースコードが簡単に書けることです。その反面ほとんどの計算処理がブラックボックス的なので計算処理がクリアではないです。
実際に数値シミュレーションする(1)
ソースコード
パラメーターなどは結果の部分で解説します。
import numpy as np
import tensorflow.compat.v1 as tf
tf.disable_v2_behavior()
#TensorFlowのバージョンによってうまくいかない場合がある(互換性が原因と思われる)
#import tensorflow as tf
import matplotlib.pyplot as plt
import time
from scipy.integrate import solve_ivp
%matplotlib inline
def main():
x=tf.placeholder(dtype=tf.float64,shape=(3),name='x')
class SIR(object):
def __init__(self, beta=0.00218, gamma=0.452, **kwargs):
self.beta = float(beta)
self.gamma = float(gamma)
def __call__(self, t, x):
dx_dt = -self.beta * x[0]*x[1]
dy_dt = self.beta*x[0] *x[1] - self.gamma*x[1]
dz_dt = self.gamma*x[1]
dX_dt = tf.stack([dx_dt, dy_dt, dz_dt])
return dX_dt
f_sir = SIR()
fx =f_sir(None, x) # ODE
x0 = np.array([762, 1, 0], dtype=np.float64) # initial value
sess = tf.Session()
sess.run(tf.initializers.global_variables())
print('###value of f(x)###')
print(sess.run(fx, feed_dict={x:x0}))
dt=0.1
tstart=0.0
tend=30.0
ts=np.arange(tstart, tend+dt, dt) # 0.1ステップで30日分を出力させる
def f_sir_tf(t,xt):
return sess.run(fx, feed_dict={x: xt})
start_time =time.time() # measurement of time
sol_sir = solve_ivp(fun=f_sir_tf,
t_span=[tstart, tend], y0=x0, t_eval=ts)
integration_time_tf = time.time() - start_time
t_lo = sol_sir['t'] # 各ステップの時刻を取得
x_lo = sol_sir['y'] # 各ステップのx(t)の値を取得
#処理時間を計算する
print("processing time (tf) %.5f" % integration_time_tf)
fig=plt.figure(figsize=(14,8))
ax=fig.add_subplot(111)
ax.plot(t_lo,x_lo[0],'*',color='blue',alpha=0.2,label="Susceptible",markersize=10)
ax.plot(t_lo,x_lo[1],'^',color='red',alpha=0.2,label="Infected",markersize=10)
ax.plot(t_lo,x_lo[2],'.',color='green',alpha=0.2,label="Recovered",markersize=10)
ax.grid()
ax.legend()
ax.set_xlabel('day')
ax.set_ylabel('number')
if __name__ == "__main__":
main()
シミュレーション結果
今回のシミュレーションは以下のスライドのp14を参考にしました。
http://www.actuaries.jp/lib/meeting/reikai20-7-siryo.pdf
1978年に調査された、ある寄宿学校のインフルエンザの事例です。合計763人中、$S＝762$人で$I＝1$人を初期条件とします。また感染率$\beta$、回復率$\gamma$については
$$\beta=0.00218,\gamma=0.452$$
としてシミュレーションしております。パラメーターについてもう少し詳しく言うと、感染率$\beta$は一日あたりに感染する率で回復率$\gamma$は非感染になるにはどれだけの時間を要するのかに相当します。
ちなみに、この時の$S$と$I$の関係をグラフにすると以下のようになります。右下から左側に行くイメージで原点に収束致します。
解釈
緑の回復した人数$R$はここではあまり重要ではないです。一番重要なのは赤の感染者数$I$です。今回の場合初期状態$I=1$の状態から、感染者と接触した$S$が$I$になっていき急激に増え1週間後にピークを迎え、大体14~15日後に0に収束していってる様子が伺えます。実際に参照スライドp14の白丸プロットのようにうまくフィッティングされていることが分かります。
実際に数値シミュレーションする(2)
インフルエンザについては数値計算し終わったので次は「コロナウイルス(Coronavirus)」です。インフルエンザのところでは説明しませんでしたが、感染症を題材として扱う上で重要なパラメーターである「基本再生産数$R_0$」から解説します。
基本再生産数
定義としては、「全感染期間において1人の感染者が生み出す2次感染者数の期待値のこと」です。例えば$R_0=2$であれば一人が平均で2人に感染させるということです。詳しくは語りませんが$$R_0<1$$であれば感染流行は収束することが一般に知られております。SIRモデルの2番目の式で$右辺＝0$として計算すれば出せます。
今回の$R_0$について以下のURLを引用します。
https://wired.jp/2020/01/26/scientists-predict-wuhans-outbreak-will-get-much-worse/
大きな不確定要素に、2019-nCoVの感染力がある。実際、どの程度の感染力があるのだろうか？
リードのモデルでは、1人の感染者が感染させる可能性のある人数（ウイルスの基本再生産数）は3.6〜4.0と推定されている。比較として、SARS（重症急性呼吸器症候群）は2〜5であり、ヒトへの感染力が最も高い麻疹（はしか）は、なんと12〜18だ。
今回のコロナウイルスを基本再生産数$R_0=3.4$で計算してみます。計算は省きますが感染率$\beta$は以下の式で算出できます。
$$\beta=\frac{\gamma R_0}{S(0)}$$
ここで初期条件について考えます。コロナウイルスが蔓延している武漢市は1108万人の中国の大都市ですので
$$(S,I,R)=(11079999,1,0)$$
として計算することにします。
すると$\beta=2.19\times10^{-8}$が得られます。ここで$\gamma=0.071428$(1/14で2週間に相当)とおきました。$S(0)$は11,079,999としております。
結果
ソースコードはインフルエンザの場合と基本的には同じです。パラメーターだけ変えてあります。
解釈
初期条件では最初の感染者が発見された時点を想定しております。現実的には2019年12月1日だと思えば良いと思います。今ちょうど2月初めなので大体2ヶ月が経過したと思うと60日時点です。と思うとめちゃめちゃ恐ろしいです。。。今回はR=3.4というかなり大きな数字を使ったので最悪のケースを想定しております。実際WHO(世界保健機関)によると$R_0$=2.2~2.8だそうなので、もっと早く収束するはずです。第一感染者発見時をいつととるのかもキーポイントです。もし今時点が感染者のピーク(グラフで100日付近)だと仮定とすればここから感染者数が0に収束していくのに40日程度かかると見込まれます。お気づきの方も多々いると思いますが、このモデルでは死者は考えておりません。したがってS+I+Rの系が保存されます。より現実的なモデルを考える場合はSEIRモデルなどが挙げられます。
まとめ
本投稿
初めてQiitaに投稿しました。今回は「感染モデル」の中で最も基本的なものであるSIRモデルについて数値計算しました。実際に用いたパラメーターは過去の事例に習いましたが、パラメータ推定などを詳細にやる場合はもっと高度な数学が必要になります。時間があれば今後そういった記事も書いてみようと思います。個人としてはMachine Learningや競技プログラミングに興味があるのでそちらの備忘録等含めたリポジトリにしようとも思っております。
展望
- SEIRモデルの数値計算(SIRに潜伏期間を含めたモデル)
- ネットワークモデルにおけるSIR
- TensorFlowでいろいろ遊ぶ
- 本気で手洗いうがいをする