モチベーション
今回、メルケルさんやジョンソンさんが言及している内容や厚生労働省が発表している内容を正しく理解するひとつの端緒となるべく、感染が広がる様子をシミュレーションをしてみました。
なぜ Qiita?
今の時期に(2020年4月1日)こんなシミュレーションをするのはけしからんという向きもあるでしょうが、Qiita を読むような人なら正しくその趣旨や内容を理解すると思ったので。
時期が来たら消すかもしれません。
内容は？
ソースは github で公開しています。(https://github.com/ryos36/sim-covid-19) 。医学的根拠や数学的根拠(とくに統計)はないことに注意してください。また、結果はあくまでシミュレーションです。パラメタを変えれば結果は自分の求めるものに合わせられてしまいます。
プログラムの中身
プログラムは Python で書いてあります。かなり RAPIDに(やっつけで)つくったのと十分なテストがされてないため、バグはあるでしょう。そこそこの結果が出ているので現時点では個人的に良しとしています。パラメタは変えられます。一応、github の Param.md (https://github.com/ryos36/sim-covid-19/blob/master/Param.md) に説明をつけてあります。たぶんそれだけじゃわかりません。興味のある方はソースを読んでください。ソースの説明をするつもりは現在ではありません。
ソースは汚い
クラス化されておらず(int の方がスピード上がりそうだったから。検証した訳ではない)、その上、無駄にループがあります。
前半のループ(各所の更新)と後半のループ(統計情報を集める)と絵を書くループ、移動をシミュレーションするループは恐らく１つのループ内で達成できます（今考えたらちゃんと整理すればシミュレーション速度が4倍になるじゃないか）。
シミュレーションが遅い
高速化をねらってないのでシミュレーションには時間がかかります。将来的には都市間シミュレーションもしたいところですが、そうするなら、クラス化してちゃんと整理した方がよいでしょう。
CUDA とかで書いたらもっとはやくなりそう。
UI がないに等しい
グラフ化もしてないし、パラメータもフレキシブルに与えにくくなっています。だれか JavaScript とかと連携させるなりして Web アプリにしてください。
シミュレーションをする。
シミュレーションの前にモチベーションともなったよさそうな情報(正しそうな情報源)を列挙しておきます。
- メルケル独首相「全人口の6～7割感染も」(https://www.nikkei.com/article/DGXMZO56691720S0A310C2000000/)
- 英政府、独自の新型コロナ「集団免疫」戦略を修正へ (https://www.technologyreview.jp/nl/the-uk-is-scrambling-to-correct-its-coronavirus-strategy/)
- Impact of non-pharmaceutical interventions (NPIs) to reduce COVID-19 mortality and healthcare demand (https://www.imperial.ac.uk/media/imperial-college/medicine/sph/ide/gida-fellowships/Imperial-College-COVID19-NPI-modelling-16-03-2020.pdf)
COVID-19 の感染は恐らく２つのシナリオがあって
- インフルエンザのようにいずれ収束する
- ６～７割感染して（するまで）収束する
で、今回、後者の集団免疫を獲得するまで感染が広がるという立場をとりました。
分子モデルのようなものでシミュレーションしている人もいる。(元ソースはワシントンポストらしい）
- https://www.yomiuri.co.jp/topics/covid19/wp/
簡単にシミュレーションしてみる
一応、フレキシブルになっているので巷で流れている”噂”を数値化してシミュレーションしてみました。このシミュレーションで重要なパラメタは次の通り。
- R0=2.8
- 潜伏期間 5日（中央値） <- プログラム上確率的にはなっていない
- 治るまで 7日
- 重症化率 5%
- 重症化した人がお亡くなりになる率 5%
- 重症化している日数 20日
R0 は医学的な用語に似せてあるがそれを指しているわけではないことに注意。モデルは中学生レベルの数学。治るまで7日かかり、そのまわりに8か所(これはプログラム的な方向。現実を表していないところに注意)あり、その状況で 2.8 人に移すなら、単純に 2.8/(7*8) を確率とみなして計算します。
移動も考慮しています
近距離移動と遠距離移動をシミュレーション可能なようにしてあります。
- jump_distance = 30
- jump_distance_long = int(width * 0.5) -> 1920 * 0.5 = 960
結果の一例を示す。
グラフにすると次のようになりました。
繰り返しになるが、これはあくまでシミュレーションです。そして、こうなるであろうという予測のもとにパラメタを調整しています。
これの意味するところは何か？
このモデルの意味するところは単純に R0 値とシミュレーション対象の人口の2つのパラメータから、最終的にどれだけの人が感染するか？を表しているに過ぎません。シミュレーションしなくても実は計算でわかります。人口が 1920*1080 = 2073600 で、感染率が 80%、重症化率が 5%、重症化の人がなくなる率が 5% (= 致死率 0.25 % に相当) なら 4147.2 人が犠牲になります。
インフルエンザの入院率と致死率はもっとひくいらしい。( https://www.hosp.med.osaka-u.ac.jp/home/hp-infect/file/ictmon/ictmon162.pdf) 同じ人口に感染率 16%、致死率 0.001% を割り当ててみると、3.46 人。このシミュレーションとは別にある確率で犠牲者が予測されるのはどの感染症も同じです。
このシミュレーションが示唆する第一の点は犠牲者の数ではなくその期間です。収束までに 300 日を超す時間が必要です。示唆という言葉を使いましたが、シミュレーションがそう設計されているので示唆ではなくそうシミュレーションしたというのが正しいかもしれません。
犠牲者の数はどう見積もるか？
ここからはセンシティブな話になります。同じ事象でも、社会が受け入れなければならない犠牲と個人が受け入れることが出来ない犠牲は別でしょう。感情の問題も多々含まれます。ここではそれを横に置いておいてシミュレーションのみをします。
シミュレーションでは病床数も考慮に入れました。重症化する人とあらかじめ設定されたベット数(人口の 1/1000 を想定)を比較し、ベッドの数が足りなくなった場合、確率的になくなる率を高くしています。一人一人の症状をおってシミュレーションをしているわけではありません。
ベッド数が十分にある場合とない場合を比較シミュレーションしてみます。
左が別途数不足、右が潤沢にある場合。グラフのオレンジ色が犠牲者。ベッド数が少ないため治るというより無言の帰宅をすることで重症者数が減っていくというかなり衝撃的なシミュレーション結果になってしまいました。シミュレーションの上だけの話になってほしいところです。
長距離移動を抑える
このシミュレーションでは長距離移動を抑えると感染の広がる速度はなだらかになります。このなだらかになるが曲者です。もともと 300 日を超えるシミュレーションなので、なだらかになるが意味することはシミュレーション期間が数年を要するという事を意味します。
上のグラフは重症者数の推移。400 日に及ぶ長距離移動の抑制(赤で色付けされた期間)により常にベッド数を確保可能な状況をつくることができました。完全収束はには 700 日を超す日数を要すという結果に。
断続的な長距離移動抑制(ロックダウン)は有効か？
長距離移動を抑えることが有効そうな結果は（シミュレーション上で）でましたが、なるべく生活に差し支えないようにしたいところです。「断続的に長距離移動を抑えることでそれは実現可能であろう。」という予測のもとにシミュレーションしてみます。
上の図は感染者数のグラフ。7,14日間の抑制は微々たる影響しか与えません。28日でグラフに明らかに違いは認められます。400日連続で行えばほぼ横ばいにすることができました。感染拡大後の抑制は効果が薄いこともわかりはじめます。もう少し早い段階で抑制を開始すべきでしょう。
Imperial-College-COVID19-NPI-modelling-16-03-2020.pdf との比較
Imperial-College-COVID19-NPI-modelling-16-03-2020.pdf との比較をしてみます。次に pdf から図を拝借します。
この pdf では最初は 6週間、2週間空けて、4週間の suppression(抑制) を提案しています。同じようなことをこのシミュレーションでも実施してみます。
結果の評価は非常に難しいものとなりました。青が重症者の数、オレンジが亡くなられた方の数、灰色が抑制期間です。３つのグラフで一番上は抑制が早すぎでした。一番下は遅いタイミングでしたが効果がありました。うまいタイミング で抑制を 適切な期間 行えばよいようで、ここのシミュレーション環境では人口の0.3% 程度の重症者数が出た前後をカバーするように抑制を行えばよいようでした。現実の世界ではパラメタが多すぎて適切なタイミングを計るのはかなり難しいでしょう。
シミュレーション上の課題
シミュレーション結果はプログラム上の構造に依存しています。例えば、始点が真ん中の場合と左端の場合で違った結果を算出します。地理的要因が考慮されているともとれますが、やはりプログラム上の課題でしょう。実際の都市は地続きであり、また人口密度も一様ではないので、それらが考慮されてシミュレーションされるべきでしょう。
図のように始点が真ん中だと感染の広がりが早く被害は（シミュレーション上）大きくなります。
あくまでシミュレーションであるという事
最後になるがこれはあくまでシミュレーションであるという事。何が言えるか？何が検証できるか？実際のデータが正確に探し出せてないのでなんとも難しい。
恐らく、”長期戦”であると言い始めている人がいることからも、どこかで正確なデータとより数理的、解析的な方法でシミュレーションしていると思います。それらのプログラムや数式が公開されて複数の研究者やプログラマで検証されれば、すくなくともシミュレーションという意味では理解が進むことでしょう。
直感的には人間は移動はするが引っ越しは頻繁にしないのでそういった人間の行動を数理的に完全に組み入れるのは難しい気がします。そういう意味でこれらのプログラマ的なシミュレーションは意味を持つのではないでしょうか？
1年後には単なるシミュレーションはシミュレーションだったね。とハッピーな気持ちで見返せればよいと思います。
追補
緊急事態宣言直後(2020/4/10)
新型コロナウイルスに関するQ&A（一般の方向け）から図を引用します。新型コロナウイルス感染症対策の基本方針（令和２年２月25日）ではピークに関して図のように説明していました。
おそらく、集団免疫を想定していたと思われます。（山の面積が違うのが気になりますが）
緊急事態宣言後ちょっと情報が変わり始めました。「新型コロナウイルス感染症対策専門家会議」からの情報を掲げます。
「人との接触を最低７割、極力８割まで減らすことを目指しています。」そうすれば、グラフにあるように患者数が減るということらしいです。これはどうもシミュレーションをしたわけではなく、R0 の数字に 0.2 (8割減) を掛けて、潜伏期間＋治癒までの期間すぎればこうなるという数字を示したように見えます。
北海道の実績を見ると（北海道庁から引用）横ばいになるのが精いっぱいのようです。
今後の推移を見守りたいと思います。