感染症流行の数理(SIRモデル)のグラフ化

雑記
スポンサーリンク

お久しぶりです。ブログ更新サボっていました(^_^;)

3月頃からの新型コロナウイルスのパンデミックにより、世界中が混乱に陥っていますね。私も日々ニュースをみて不安になったり、投資の含み損を見て嫌な気持ちになったりしています。 仕事もフルリモートワークとなって、オーバワーク気味で疲れが溜まっておりました。しかし、ゴールデンウィークの休みで元気を取り戻してきました。

今回の記事は気分転換とリハビリも兼ねて、雑記回です。コロナショックで不安を抱えている方への頭の体操や、プログラミングの入門となれば幸いです。

さて、今回は感染症流行の古典的モデルである「SIRモデル」をPythonでプログラミングしてみました。

SIRモデルは、すでに色々なサイトやYoutubeで解説されておりますが、感染症の流行を表すのに一番シンプルな数学的なモデルだそうです。100年近く前からあるモデルですが、未だに使えるそうです。SIRモデルの詳しい数理はWikipediaや後述の参考サイトに解説があるので、興味のある方は調べてみると面白いと思います。

スポンサーリンク

基本再生算数2.5のときのグラフ

SIRモデルにおいて、感染率βを0.25、回復率γを0.1(10日間で隔離、回復または死亡する)、つまり1人の感染者が感染を移す人数R0(基本再生算数)を2.5とした場合の流行は以下のようなグラフとなりました。

R0=2.5のときのSIRモデル

R0=2.5という数字は、専門家会議で提示されていたものです。実態はまだわからないそうですが、2~3くらいというのが見立てだそうです。

感染者が発生してから流行のピークまで約50日ほどかかり、20%以上の方が感染し、流行収束までに約90%が感染するとなりました。

日本では、重症者だけ受け入れるとしても、医療キャパシティはせいぜい数千から数万人程度でしょうから、 ほとんどの人が軽症とはいえ、放っておけば医療崩壊待ったなしですね。 再生算数を低く抑えながらみんなの免疫がつくまで長く付き合っていかざるをえなさそうですね。

SIRモデルのソースコード(Python)

グラフを表示するのに使ったPythonのソースコードです。参考サイトのものをほとんどコピーしたものですが、一部見やすく修正しています。差分方程式を解いたり、グラフ化したりするのにこんな短いコードでできるとは、素晴らしいですね。

Pythonはほとんど使ったことがないのですが、使いやすいですね。世界中で流行しているのもうなずけます。

import numpy as np
from scipy.integrate import odeint
from scipy.optimize import minimize
import matplotlib.pyplot as plt

# the effective contact rate of the disease
BETA = 0.25
# recovery rate
GAMMA = 0.1
R_nought = BETA / GAMMA

# max time(days)
MAX_DAYS = 200
# delta time
DT = 0.01

# Number of population
NUM_POPULATION = 1000
# Initial infected population
I_INIT = 1
# Initial susceptible population
S_INIT = NUM_POPULATION - I_INIT
# Initial removed population (either by death or recovery)
R_INIT = 0

# SIR Differential Equation
def diff_eq_func(y, time, NUM_POPULATION, BETA, GAMMA):
    susceptible = y[0]
    infectious = y[1]
    ds_dt = -BETA * susceptible * infectious / NUM_POPULATION
    di_dt = BETA * susceptible * infectious / NUM_POPULATION - GAMMA * infectious
    dr_dt = GAMMA * infectious

    return [ds_dt, di_dt, dr_dt]

# Solver SIR model
y_init = [S_INIT, I_INIT, R_INIT]
times = np.arange(0, MAX_DAYS, DT)
result = odeint(diff_eq_func, y_init, times, args=(NUM_POPULATION, BETA, GAMMA))

# plot
plt.plot(times, result)
plt.title('SIR model. R0: ' + str(R_nought) + '(beta: ' + str(BETA) + ', gamma: ' + str(GAMMA) + ')')
plt.legend(['Susceptible', 'Infected', 'Recovered\n(isolated/death/recovery)'])
plt.xlabel('Time(days)')
plt.ylabel('Population')
plt.show()

参考サイト

Wikipedia

The SIR epidemic model

タイトルとURLをコピーしました