感染症の流行シミュレーションを試してみよう
精密工学科の学生の皆さん．慣れない自宅での履修苦労されているかと思います．いま社会で起こっている新型コロナウイルス感染拡大の最前線の現場は病院等の医療現場です．しかし，私たちエンジニアもけっして別の世界に居るわけではありません．たとえば，マシンダイナミクスで学習するバネ・マス・ダッシュポット系の運動は微分方程式で表すことができます．実は同じように感染症の流行も微分方程式で表現をすることができます．
学科開講科目のマシンダイナミクスは機械系の振動問題についての講義が主ですが，この機会に，その中で番外編として少しだけ紹介している感染症流行のシミュレーションについて再度紹介をしたいと思います．すでにマシンダイナミクスを受講した学生も受講していない学生もいることでしょう．受講済の学生は思い出してみてください．これから受講する学生はお楽しみに！
よいエンジニアになるためには，目の前で起こっていることを批判的に（ホントかどうか？），定量的に（何秒？何％？）評価できるようになることが大切です．ずっと家にいて窮屈な日々を過ごしているかと思いますが，エンジニアとしてできることは何か？，ニュースで流れている情報は本当か？，ノートと鉛筆そして，エンジニアだから電卓（パソコン）を使ってじっくり考えてみましょう．
１．意外と単純！？微分方程式を使ったモデル化
さて本題に入ります．古典的な感染症の流行モデルにSIRモデルというのがあります．Sは未感染者，Iは感染者，Rは免疫保持者のことです．この３つの集団を考えて，その振る舞いを数式で表現することにしましょう．即ち変数として，
S：未感染者数，I：感染者数，R：回復者数（免疫保持者数）
を考えます．S + I + R は考えている系（たとえば日本）の全体の人数になります．ここでは，一度感染して回復した人には免疫ができて二度と感染しなくなると仮定します．
まず， I（感染者数）がどのように増えるのかを考えてみます．感染するシチュエーションを想像してください．感染は感染者と未感染者が接触することで起こります．したがって，「感染者という集団に属する人数(I)」が多いほど，かつ「未感染者という集団に属する人数(S)」が多いほど接触の機会が多くなります．したがって，I × Sに比例して感染者が増えやすくなると考えます．また，長い時間が経過するほど接触の機会が増えますので，いま考えている時間経過の長さΔtにも比例すると考えます．ここまでの議論から，感染者数Iは S × I × Δt に比例して増えそうだ！と，とても単純化して考えます．でもちょっとまってください．接触したら確実に感染するわけではないですよね．そこで，感染率βを導入します．感染率βとは，単位時間あたりに一人の感染者が未感染者の中から新たに感染者を作り出す割合のことです．以上より，いま考えている時間Δtの間に新たに感染する人数は β × I × S × Δt となります．
どうでしょうか？「ちょっと待てよ」と思った人もいるでしょう．たしかに，新たに感染する人数は上記のβ × I × S × Δt ですが，時間が経過するにしたがって，回復して，感染者（I）の集団から離脱する人も居るはずです．その人たちも考慮して感染者数（I）の変化を考える必要があります．その回復者はどのように考えればよいでしょうか．そもそも感染している人が居なければ回復する人も居ないので，考えている時間Δtの間に回復する人の数は感染者の数(I)に比例すると考えることにします．したがって，Δtの間に新たに回復する人の数はI × Δtに比例すると考えます．ここでも，全員がその時間の間に回復するわけではありませんから，回復率γを導入します．回復率γは単位時間あたりに感染者の中から回復する割合です．これを掛けると，Δtの間に新たに回復する人の数はγ × I × Δtとなります．
以上では新たに感染する人数はβ × I × S × Δt であり，新たに回復する人数はγ × I × Δtであることを説明しました．このことから，いま考えている時間Δtの間の感染者の変化ΔIを求めると，単純に引き算をして，
ΔI = （ β × S × I × Δt ） -（ γ × I × Δt ）
となります．
これを，Δtは限りなく小さいと考えて極限をとると，
dI/dt = β × S × I - γ × I
が得られます．これが，感染者数（I）の振舞いを表す微分方程式です．また，上記の考察から落ち着いて考えると，未感染者数については，
dS/dt = - β × S × I
回復者（免疫保持者）については，
dR/dt = γ × I
となります．ここは，ゆっくり落ち着いて考えてみたら理解できると思います．
２．計算してみよう
さて，上記の３つの式が，S：未感染者数，I：感染者数，R：回復者数（免疫保持者数）の振る舞いを表す微分方程式です．よく見ると，左辺は各集団の人数が変化する速度（変化率）を示していて，右辺はその速度が，今見ている時刻におけるSやIの値だけで決まるんだよ．ということを言っています．あとは，βとγも影響しますが，それはその感染症の固有の性質などによる値で，時々刻々変化するわけではありません．そして，この式を差分化して組み込んだExcelのファイルが，
（クリックしてください）感染症の流行のシミュレーション.xlsx
これです．差分化については今は気にしないでください．微分方程式をExcelのような数値計算（加減乗除）のできるソフトに具体的に組み込む方法のことです．
このファイルでは，１万５０００人の全体の集団を想定しており，初期の感染者が１０名です．βとγは適当な値が入っています．t の単位は日で，0.01日刻みで60日分の計算を行っています．
＜検討方法例１＞
各設定値（パラメータ）を自由に変化させてみて，どうすると終息が早まるかを検証してみてください．たとえば，いま皆さんが実践しているようにソーシャルディスタンスをとったり接触機会を減らす行為は，感染率βを減らすことにつながっています．たとえばβを80％減らしたりしてみてください．感染者数の最大値や終息するまでの期間を評価すると良いでしょう．また，場合によっては全員が感染せずに終息する条件が見つかるかもしれません．
＜検討方法例２＞
また，インフルエンザのように予防接種がすでに行われている感染症の場合には，初期の免疫保持者数を大きくすることで，その状態を再現することができます．例えば，予防接種率50％の状態を想定して，S0を7490人，R0を7500人，I0を10人として計算をしてみてください．
＜検討方法例３＞
有効な治療法が開発されて，回復するまでにかかる期間が短くなった場合を想定して，γを大きくしてみてください．
＜検討方法例４＞
感染がある程度拡大した後に，緊急事態宣言が出るなどして屋外へ出る機会を減らした場合を想定して，感染率βを途中で低くしてみてください．たとえば，３日目からβを90％減らしてみるなど．これは，エクセルの計算欄を編集する必要があるので，少し大変ですがチャレンジしてみましょう．
注意：時間刻みは大きくすると計算の精度が落ちるのであまり大きくしないでください．
３．まとめ
本記事によって皆さんのマシンダイナミクスをはじめとした力学の理解が少しでも深まるのと同時に，計算をしてみることで外出自粛をする必要性を実感していただければと思います．何か面白い計算結果が得られたりモデルの改良を行ったりした場合は，ぜひ教えてください．本モデルは最も基本的なモデルであり，さらに高度なモデルも提案されているようです．また，初期に設定されているパラメータは今回の新型コロナウイルスを想定したものではありません．今回の感染拡大を表現できるパラメータについても検討してみると良いでしょう．
先が見えず不安な日々をお過ごしかと思いますが，出来ることを考えて，前向きに過ごしていきましょう！
参考文献
平山修：Excelで試す非線形力学，コロナ社，（2008），48-55.