SIRモデルによる感染シミュレーション
SIRモデルを使った簡単な感染シミュレーションを行います(現実のコロナウイルスの状況やロックダウンの状況とはリンクしていません) 。
SIRモデル（エスアイアールモデル）は、感染症の短期的な流行過程を決定論的に記述する古典的なモデル方程式である。(ウィキペディア SIRモデル より)
の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
①感染率bを0.25、0.5、1と変えた場合
感染症への免疫がない人々の初期値は1億、感染症に現在かかっている人々の初期値を1、回復率cを0.1として計算しています。
コード
# グラフをjupyternotebook内に表示するコマンド %matplotlib inline import matplotlib.pyplot as plt #S,I,Rを計算する関数 def sir(b0, c0, time0, susceptible0, infected0, recovered0, b1, c1, time1): 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 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 (0.25, 0.1, 0, 100000000, 1, 0, 0.25, 0.1, 25) time2, susceptible2, infected2, recovered2 = sir (0.5, 0.1, 0, 100000000, 1, 0, 0.5, 0.1, 25) time3, susceptible3, infected3, recovered3 = sir (1, 0.1, 0, 100000000, 1, 0, 1, 0.1, 25) # グラフ描画オブジェクトの生成とサイズ設定 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=0.25 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=0.5 c=0.1') 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') ax3.set_xlabel('Time') ax3.set_ylabel('Populstion') ax3.legend() # グラフ並べ方の設定 fig.tight_layout()
実行結果
感染率を変えると感染者数(infected)の最大値とピークが現れるタイミングが大きく異なることが分かります。感染率を抑えることが感染爆発(オーバーシュート)の 抑制に重要であることが分かります。
②初期の感染率bを1としてtime=25において感染率bを1のまま、0.5に減らす、0.2に減らす場合。外出を大幅に制限することで感染率が下がる場合の想定です。
感染症への免疫がない人々の初期値は1億、感染症に現在かかっている人々の初期値を1、回復率cを0.1として計算しています。
コード
# グラフをjupyternotebook内に表示するコマンド %matplotlib inline import matplotlib.pyplot as plt #S,I,Rを計算する関数 def sir(b0, c0, time0, susceptible0, infected0, recovered0, b1, c1, time1): 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 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) time2, susceptible2, infected2, recovered2 = sir (1, 0.1, 0, 100000000, 1, 0, 0.5, 0.1, 25) time3, susceptible3, infected3, recovered3 = sir (1, 0.1, 0, 100000000, 1, 0, 0.2, 0.1, 25) # グラフ描画オブジェクトの生成とサイズ設定 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.5 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') ax3.set_xlabel('Time') ax3.set_ylabel('Populstion') ax3.legend() # グラフ並べ方の設定 fig.tight_layout()
実行結果
time=25で感染率を1から0.2に下げることで、感染者数の最大値を大幅に減らせることが分かります。外出自粛や「3つの密」を避けることが大事ということが分かります。
③初期の感染率bを1として0.2に減らす場合。減らすタイミングをtime=15,25,35と変えています。外出自粛やロックダウンのタイミングで感染者数がどのように変わるかを想定しています。
感染症への免疫がない人々の初期値は1億、感染症に現在かかっている人々の初期値を1、回復率cを0.1として計算しています。
コード
# グラフをjupyternotebook内に表示するコマンド %matplotlib inline import matplotlib.pyplot as plt #S,I,Rを計算する関数 def sir(b0, c0, time0, susceptible0, infected0, recovered0, b1, c1, time1): 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 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, 0.2, 0.1, 15) time2, susceptible2, infected2, recovered2 = sir (1, 0.1, 0, 100000000, 1, 0, 0.2, 0.1, 25) time3, susceptible3, infected3, recovered3 = sir (1, 0.1, 0, 100000000, 1, 0, 0.2, 0.1, 35) # グラフ描画オブジェクトの生成とサイズ設定 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 --> b=0.2 c=0.1 @15') 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 @35') ax3.set_xlabel('Time') ax3.set_ylabel('Populstion') ax3.legend() # グラフ並べ方の設定 fig.tight_layout()
実行結果
感染率を1から0.2に下げるタイミングで、感染者数の最大値を変わることが分かります。特にtime=35の場合最大値が著しく大きくなり、感染率を下げるタイミングは感染者数のピークよりも前に持ってくる必要があることが分かります。
④time=25で回復率が初期値0.1から変化する場合のモデルです。0.1のまま、0.3に上昇、0.9に上昇する場合の計算です。有効な薬が見つかり回復率が上昇した場合を想定しています。
感染症への免疫がない人々の初期値は1億、感染症に現在かかっている人々の初期値を1、感染率bを1として計算しています。
コード
# グラフをjupyternotebook内に表示するコマンド %matplotlib inline import matplotlib.pyplot as plt #S,I,Rを計算する関数 def sir(b0, c0, time0, susceptible0, infected0, recovered0, b1, c1, time1): 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 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) time2, susceptible2, infected2, recovered2 = sir (1, 0.1, 0, 100000000, 1, 0, 1, 0.3, 25) time3, susceptible3, infected3, recovered3 = sir (1, 0.1, 0, 100000000, 1, 0, 1, 0.9, 25) # グラフ描画オブジェクトの生成とサイズ設定 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 --> b=1 c=0.1 @25') 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=1 c=0.3 @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=1 c=0.9 @25') ax3.set_xlabel('Time') ax3.set_ylabel('Populstion') ax3.legend() # グラフ並べ方の設定 fig.tight_layout()
実行結果
time=25において回復率を上昇させることで劇的に感染者数の最大値が下がります。特効薬が見つかることが非常に重要であることが分かります。
