素人がパンデミックをなるべく簡単なモデルでシミュレーションしてみる

Python+numpyでパンデミックをなるべく単純にシミュレーション

1週間ほど前、遂に我が家の近所の大型商業施設で、新型コロナウィルスの感染者が出てしまいました。
うちの地元では誰でもショッピングに行くところで、今までどこか他人事なところ(自分には感染しないよ)もあったのですが、急にリアリティを感じる様になってきました。
報道や自治体の発表を見る限り、同じ施設利用者で追加で感染が確認されたのはこの1週間で1名だけですが、毎日多数の人が利用する商業施設でそんなに感染少ないわけないよな..と。
例えば1人の感染者が毎日別の1人に感染させると、感染者数は倍々ゲームで増えて行って、単純計算だと33日目には世界人口超えてしまうわけですが、世界の状況を見るにそこまでまだ酷い状況では無さそうです。
とは言え、WHOは既にパンデミックと言い始めました。

自分は感染症の専門家でも無いし、数理シミュレーション屋さんでも無いですが、以下の様なモデルで極簡単にシミュレーションをしてみました。
あくまでも素人の思考実験です。グラフの形だけは割とそれっぽくなりましたが、現実はもっと複雑なんだろうと思います。

1. 感染者1名が1日に平均何名に感染させてしまうかをパラメーター化する
   
   1日0.2人に感染させるとする。(これは何度もシミュレーションして最初の感染者が出た日(Day0)から3ヶ月くらいの感染者発生数が割と現実のオーダーに近い(あくまでも桁が合う程度)値になる様に調整しました)
   
2. 感染者の平均感染期間を設定し、その期間のみ上記のパラメーターに基づき、感染者が増える
   
   取り敢えず平均を10日に設定。(要するに1.と合わせて一人の感染者が10日間で二人に感染させてしまうという単純な計算です)
   
3. 感染者のうち、20%が重症化、同時に感染者として検出されカウントされる
   
4. 重症化するのは感染してから平均7日後とする
   
5. 感染者は一度感染したら、2回以上感染しない事とする(新型コロナウィルス感染でどの程度免疫が獲得できるか明らかではありませんが、取り敢えずこのように仮定します)
   
   この考えで、(1-トータル感染者数/全世界人口(ざっくり77億人))を最初のパラメーターに毎日乗じて行くことで、感染者数がある程度増えると増加率が抑制される様にする。
   
6. 総感染者数(感染検出者ではない、症状が出ないあるいは軽症で風邪と区別付かない人も含む)が世界人口の8割に達したらシミュレーション停止
   
7. 季節、気温によるウィルスの感染力変化は考慮しない(データが無いので出来ない)

ちなみに重症者は隔離されて、新規感染者を発生させにくくなるとかもあるはずですが、それを入れると複雑化しそうなのと、結局は1.のレートの調整に帰結すると思うので考慮していません。

上記の前提で計算した最初の感染者が発生したDay0〜Day30までのシミュレーションをグラフ化したものです。

day0〜day30まで

上段左側がトータル感染者数、上段右側が日々の新規感染者数、
下段左側がトータルの感染検出者数(重症化した人)、下段右側が日々の新規感染検出者数です。
最初に誰かが感染してから、実際に感染者が確認される確率が50%以上になるまで24日くらい掛かります。

次にday100までのシミュレーションです。

day0〜day100まで

3/14のWHOの日報によると、
https://www.who.int/docs/default-source/coronaviruse/situation-reports/20200314-sitrep-54-covid-19.pdf
Globally
142 539 confirmed (9769 new)
だそうです。これに近い数字(1日の新規感染者確認数が9769人もしくはトータルの感染確認者数が14万人)をシミュレーション結果から探すとday94〜99あたりがオーダー的には近いものになります。
現実のコロナの最初の感染者が見つかったのが11/7という説がある様で、既にそこから4ヶ月経過しています。下記だと3ヶ月なので1か月くらいずれていますが、あくまでも桁が違わないというレベルで合わせこむとこんな感じになります。

day = 94, new_infected = 125675.626608, new_detected = 9394.323170,total_infected = 958165.299817,total_detected = 71617.911461
day = 99, new_infected = 253795.579841, new_detected = 18974.472213,total_infected = 1935202.998709,total_detected = 144659.256943

さて、このシミュレーションを200日くらい続けると、以下の様になります。

day0〜day187

187日目総感染者数が世界人口の8割を超えてシミュレーションが止まりました。

day = 187, new_infected = 8630412.553553, new_detected = 4937914.125984,total_infected = 6164594617.184158,total_detected = 1213025199.994220
Total infected count exceeds World population. x 0.800000

まさにパンデミックです。ただしこのシミュレーションは世界中の人々がほぼ均一な密度で分布し触れ合う様な仮定のもとに成り立っており、国境や海による人の流れの抑制や長距離移動に伴う遅延とか全く考慮していないので、現実はもう少し緩やかであろうことを期待したいです。
とりあえずこのシミュレーション結果だと、ピーク時の新規感染検出者(重傷者)/日が5千万人くらいになってしまい、こうなる遥か手前で実際には医療崩壊しそうです。

でも、これは、あくまでも素人が作った単純なモデルによるシミュレーション結果です。
まあ、小松左京先生の「復活の日」とかのフィクションに使うにはいい感じのグラフだと思います。

これだけだと救いが無いので、例えば適当なタイミングで各国政府が日々の感染率抑制対策を講じて、これが上手く行った場合のシミュレーションを入れてみました。
具体的にはday90から日々の感染者発生レートをそれ以前の4割に抑える。(感染者一人あたり0.2人/日の割合で新規感染者を発生させていたのを、0.2x0.4=0.08人まで抑制する)
この意味は感染者が感染期間の10日間で他の人に感染させる人数を平均1人未満に抑えるというものです。

day90から感染者発生率を抑制

現実がこういう感じで収束すると良いのですが。（ただ、これだとday90以降、収束に更に3ヶ月かかる)

シミュレーションソースコード: バグってるとか計算が変だとかのコメントは歓迎です。あまり検算、検証していません。

# -*- Coding: utf-8 -*-

# numpy
import numpy as np

# Matplotlib
import matplotlib as mpl
mpl.use('TkAgg')  # or whatever other backend that you want
import matplotlib.pyplot as plt

# シミュレーション日数
SIM_DAYS = 180
# 感染後、陰性になるまでの平均日数
AVG_RECOVERY_DAYS = 10
# 重症となる割合
SEVERE_PATIENT_RATE = 0.20
# 重症化するまでの平均日数
AVG_SEVERE_DAYS = 7
# 一人の陽性者が1日あたりに感染させる平均人数
AVG_INFECTION_RATE_PER_PERSON = 0.2
# 世界人口
WORLD_POPULATION = 7700000000
# シミュレーション停止レシオ(対世界人口)
SIM_STOP_RATIO=0.8

# 移動制限等の対策パラメータ
# 対策有無
CONTROL_MEASURE=False
# 対策を開始する日
CONTROL_START_DAY=60
# 対策による感染率低減
CONTROL_FACTOR = 0.4

def controlFactor(day):
    global CONTROL_MEASURE
    global CONTROL_START_DAY
    global CONTROL_FACTOR
    retvalue = 1.0
    if CONTROL_MEASURE:
        if day >= CONTROL_START_DAY:
            retvalue = CONTROL_FACTOR
    return retvalue

# 日別の新規感染者
new_infected = np.zeros(SIM_DAYS+1)
new_infected[0] = 1
# 日別の累積感染者
total_infected = np.zeros(SIM_DAYS+1)
total_infected[0] = 1
# 重症化して陽性判定される患者数
new_detected = np.zeros(SIM_DAYS+1)
# 重症化して陽性判定される患者数の累積
total_detected = np.zeros(SIM_DAYS+1)


sim_stop_day = SIM_DAYS
# シミュレーション
for day in range(1,SIM_DAYS + 1):
    # その日新規に感染する人数
    day_idx = (day - AVG_RECOVERY_DAYS) if day - AVG_RECOVERY_DAYS >= 0 else 0
    ctrl_factor = controlFactor(day)
    infected_cnt = np.sum(new_infected[day_idx:day] * AVG_INFECTION_RATE_PER_PERSON * ctrl_factor * (WORLD_POPULATION - total_infected[day - 1])/WORLD_POPULATION)
    new_infected[day] = infected_cnt
    total_infected[day] = total_infected[day - 1] + infected_cnt

    day_idx2 = day -  AVG_SEVERE_DAYS
    if day_idx2 >= 0:
        new_detected[day] = new_infected[day_idx2] * SEVERE_PATIENT_RATE
    total_detected[day] = total_detected[day - 1] + new_detected[day]

    print('day = %d, new_infected = %f, new_detected = %f,total_infected = %f,total_detected = %f' % (day,infected_cnt,new_detected[day],total_infected[day],total_detected[day] ))

    if total_infected[day] >= WORLD_POPULATION * SIM_STOP_RATIO: # 感染者が世界人口x停止レシオ超えたら停止
        print('Total infected count exceeds World population. x %f' % SIM_STOP_RATIO)
        sim_stop_day = day
        break



# グラフ描画
days = np.arange(sim_stop_day + 1)
fig = plt.figure()
ax1 = fig.add_subplot(221)
ax2 = fig.add_subplot(222)
ax3 = fig.add_subplot(223)
ax4 = fig.add_subplot(224)

ax1.plot(days[0:sim_stop_day + 1],total_infected[0:sim_stop_day + 1],label='Total Infected')
ax2.bar(days[0:sim_stop_day + 1],new_infected[0:sim_stop_day + 1],label='Newly Infected')
ax3.plot(days[0:sim_stop_day + 1],total_detected[0:sim_stop_day + 1],label='Total Detected')
ax4.bar(days[0:sim_stop_day + 1],new_detected[0:sim_stop_day + 1],label='Newly Detected')

ax1.legend()
ax2.legend()
ax3.legend()
ax4.legend()

plt.show()