SIRモデルによる感染シミュレーション(一度低下した感染率が再上昇する場合)

SIRモデルによる感染シミュレーションで作成したプログラムを使って、ロックダウンや外出制限などで一度感染率が低下、その後ロックダウンが解除され感染率が上昇した場合のSIRモデルによる感染者数シミュレーションを行います。

SIRモデル（エスアイアールモデル）は、感染症の短期的な流行過程を決定論的に記述する古典的なモデル方程式である。(ウィキペディア　 SIRモデルより)

Susceptible 感染症への免疫がない人々
Infected 感染症に現在かかっている人々
Recovered 感染症から回復して感染症への免疫が生じた人々

の3項目について時間変化を計算するモデルで
時間が t から t+1 へ経過した際のそれぞれの増減を以下式で計算します。
S(t+1) - S(t) = -b × S(t)/S(0) × I(t)
R(t+1) - R(t) = I(t) × c
I(t+1) - I(t) = b × S(t)/S(0) - I(t) × c
ここでbは感染率(単位時間に1人の感染者がb人に伝染す。ただし免疫のある回復者には伝染らない)、cは回復率(単位時間に感染症から回復する割合) となります。

環境

windows10 home
Anaconda 3/ jupyter notebook 5.6.0
Python 3.7.0

コード

感染症への免疫がない人々の初期値は1億、感染症に現在かかっている人々の初期値を1、初期の感染率bを1.0、回復率cを0.1として ①bが変化しない場合 ②time=25でbが1から0.2に低下した場合(「接触八割減」のモデル) ③time=25でbが1から0.2に低下し、time=60でbが0.2から0.5に上昇する場合(「活動再開」のモデル) について計算しています。

# グラフをjupyternotebook内に表示するコマンド
%matplotlib inline  

import matplotlib.pyplot as plt

#S,I,Rを計算する関数
def sir(b0, c0, time0, susceptible0, infected0, recovered0, b1, c1, time1, b2, c2, time2):
    b = b0
    c = c0
    time = [time0]
    susceptible = [susceptible0]
    infected = [infected0]
    recovered = [recovered0]
    for i in range(1, 200):
        time.append(i)
        # 設定したタイミングでb,cを変更する場合
        if time[i] == time1:
            b = b1
            c = c1
        if time[i] == time2:
            b = b2
            c = c2
        recovered_delta = infected[i-1] * c
        susceptible_delta = -b * susceptible[i-1] / susceptible[0] * infected[i-1]
        infected_delta = -susceptible_delta - recovered_delta
        susceptible.append(susceptible[i-1] + susceptible_delta)
        infected.append(infected[i-1] + infected_delta)
        recovered.append(recovered[i-1] + recovered_delta)
    return time, susceptible, infected, recovered

# 異なる3つの設定で関数sirを呼び出し、計算結果を得る
time1, susceptible1, infected1, recovered1 = sir (1, 0.1, 0, 100000000, 1, 0, 1, 0.1, 25, 1, 0.1, 25)
time2, susceptible2, infected2, recovered2 = sir (1, 0.1, 0, 100000000, 1, 0, 0.2, 0.1, 25, 0.2, 0.1, 60,)
time3, susceptible3, infected3, recovered3 = sir (1, 0.1, 0, 100000000, 1, 0, 0.2, 0.1, 25, 0.5, 0.1, 60,)

# グラフ描画オブジェクトの生成とサイズ設定
fig = plt.figure(figsize = (6, 8))

# 3つのグラフをfigに追加(グラフを縦に3つ並べる)
ax1 = fig.add_subplot(3, 1, 1)
ax2 = fig.add_subplot(3, 1, 2)
ax3 = fig.add_subplot(3, 1, 3)

# グラフ1の設定
ax1.plot(time1, susceptible1, label='susceptible')
ax1.plot(time1, infected1, label='infected')
ax1.plot(time1, recovered1, label='recovered')
ax1.set_title('b=1 c=0.1')
ax1.set_xlabel('Time')
ax1.set_ylabel('Populstion')
ax1.legend()

# グラフ2の設定
ax2.plot(time2, susceptible2, label='susceptible')
ax2.plot(time2, infected2, label='infected')
ax2.plot(time2, recovered2, label='recovered')
ax2.set_title('b=1 c=0.1 --> b=0.2 c=0.1 @25')
ax2.set_xlabel('Time')
ax2.set_ylabel('Populstion')
ax2.legend()

# グラフ3の設定
ax3.plot(time3, susceptible3, label='susceptible')
ax3.plot(time3, infected3, label='infected')
ax3.plot(time3, recovered3, label='recovered')
ax3.set_title('b=1 c=0.1 --> b=0.2 c=0.1 @25 --> b=0.5 c=0.1 @60')
ax3.set_xlabel('Time')
ax3.set_ylabel('Populstion')
ax3.legend()

# グラフ並べ方の設定
fig.tight_layout()