感染症流行過程の数理モデルをPythonで解く

注意！この記事はなんとなくやってみた以上の情報を含みません。

感染症数理モデル

昨今のCOVID19事情を受けて色々感染症の数理モデルの話は出ていると思うわけですが、それをPythonのodeintで解いて、そのままmatplotlibで可視化してみよう、という事です。
くどいようですがこの記事のなんかを適当に実装した結果からCOVID19に関して何かのお気持ちを表明しては駄目だぞ。

で、解くのは一番簡単なSIRモデルとしまして、

\begin{align}
\frac{\mathrm{d} S}{\mathrm{d} t} & = - \beta I S\\
\frac{\mathrm{d} I}{\mathrm{d} t} & = \beta I S - \gamma I\\
\frac{\mathrm{d} R}{\mathrm{d} t} & = \gamma I
\end{align}

としておきます。
$\beta$は感染率(/(人x時間))で、病気にかかった人間ひとりが単位時間あたりに感染す事のできるまだ伝染していない人の人数、$\gamma$は治癒率(/時間)で、治癒までにかかる時間の逆数です。
$S$はまだ感染した事がなく、かつ感染能力も無い人間の数、$I$は感染しておりかつ他人に感染させる能力をもつ人間の数、$R$は隔離・治癒・死亡などで感染しているあるいはいたが他人に感染させる能力を失った人間の数です。

Python code

面倒くさいので特に解説とかは無しのコードベタ貼り。

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider, Button

def solve(value, t, foi, rr):
	# fio = force of infection (beta)
	# rr = recovery rate
	S, I, R = value
	dSdt = - foi * I * S
	dIdt = foi * I * S - rr * I
	dRdt = rr * I
	return [dSdt, dIdt, dRdt]


span = 100
t = np.linspace(0, span, 256)
initial_value = [0.9, 0.1, 0] # initial values for S, I, R
solution = odeint(solve, initial_value, t, args = (0.5, 0.5))

fig, ax = plt.subplots(figsize=(12, 6))
plt.xlim(0, span)
plt.ylim(0, 1)
plt.subplots_adjust(bottom = 0.2)

data_S, = plt.plot(t, solution[:, 0], "blue", label = "Susceptible")
data_I, = plt.plot(t, solution[:, 1], "red", label = "Infectious")
data_R, = plt.plot(t, solution[:, 2], "orange", label = "Removed")

plt.legend(loc = "best")

rr_slider = Slider(plt.axes([0.2, 0.02, 0.65, 0.02]), "recovery rate", 0, 1, valinit = 0.5, valstep = 0.1)
foi_slider = Slider(plt.axes([0.2, 0.05, 0.65, 0.02]), "force of inf.", 0, 1, valinit = 0.5, valstep = 0.1)

def update(val):
	rr = rr_slider.val
	foi = foi_slider.val
	solution = odeint(solve, initial_value, t, args=(foi, rr))
	data_S.set_ydata(solution[:, 0])
	data_I.set_ydata(solution[:, 1])
	data_R.set_ydata(solution[:, 2])
	fig.canvas.draw_idle()

rr_slider.on_changed(update)
foi_slider.on_changed(update)

plt.show()

実行結果

こんなのが出てくるんじゃないかなと思います。 ちなみに感染率と治癒率はどの程度の値が正しいのかよくわからなかったので、適当につまみを作って値をいじれるようにしておいたので、遊んでみてください。

で

まあ遊んでみてください。
あなたもパラメーターをうねうねさせながら微分方程式を解きたくなったらぜひパクってください。