感染症の流行モデル SIRモデル
つい最近知ったのだけれど、感染症の流行過程をあらわす古典的なモデルにSIRモデルというものがあるそうだ。
Sを感染感受者（非感染者）、Iを感染者、Rを回復者（隔離者）とするとき、SIRモデルは次の連立常微分方程式で表される。
ここで、βは感染率、γは回復率（隔離率）であり、tはある時刻を起点にした時間。
（１）式の第１式、第２式、第３式の辺々を足し合わせると、
になるので、
そこで、
とおこう。Kは、人の出入りがない、孤立した地域の総人口だと考えればいい。
特に、感染症の流行初期においては、感染者が少なく、
と近似することができ、感染者数は次のマルサスモデルに従う。
よって、流行初期において、感染者が増加するためには
でなければならない。
このパラメータR₀を基本再生産数とよび、R₀>1のとき、流行初期の段階で感染は拡大し、感染者数は指数関数的に増加し、R₀<1ならば、感染の流行は発生せず、感染者数は下の図のように自然に減衰してゆく。
連立常微分方程式（１）は、一般に、解けない。より正確に表現すると、（１）の解であるS(t)、I(t)、R(t)を初等的なtの関数の組み合わせで表すことができない。
しかし、γ=0で、R(0)=0のとき、
となるので、これから
となり、これを（１）の第２式に代入すると、
ここで
という解が得られる。
したがって、
回復率（隔離率）γ=0のときは、住民の全てが感染しまう。
一度感染すると誰一人として回復することがなく、しかも、感染者は非感染者に永遠にうつし続けるのだから、この結果は当然。
なお、（６）式の分子分母をN₀で割り、βN=r、さらに、I(t)=Nとおくと、（３）式は、次のロジスティック関数になる。
このとき、Kは環境収容力とか何とか呼ばれるにゃ。
ネムネコは、少し前に、新型コロナウイルスの感染者数は、概ね、ロジスティック関数（７）に従って増加し続けると記事に書いたと思うけれど、これは回復率γ=0としたSIRモデルの解ということになる。
この記事はうっかり消してしまったようだけれど…。
上の図ではK=10⁵と書いてあるけれど、K=5×10⁵が正しいにゃ。
中国当局の発表によると、中国での感染者数は約８万人なので、ネムネコの予想より随分と少ないけれど、実際の感染者は発表の１０倍ともいわれているので、と苦しい言い訳をする。
S(t)とI(t)とR(t)の関係ならば、
（１）の第２式を第１式で割ると、
この両辺をSで積分すると、
t=~のとき
となるので、
よって、
また、（２）式より
となるので、
これを解いたといっていいのか、微妙だけれど(^^ゞ。
また、感染者R(t)が最大になるのは、
の時だにゃ。
この時の感受性保持者の数は、γとβの値だけで決まるんだケロよ。
さて、γ=0といった特殊なケースでないと、非線形の連立常微分方程式（１）の解をtの関数として求めることはできないので、これを数値的に解こうじゃないか。
そこで、
という初期条件で、
感染率β=0.05に固定し、回復率γ=0、1、２の場合について――この数値が適当なものであるかどうかについては知らない！！ーー、
オイラー法を用いて、Δt=0.1として、数値的に解いてみたにゃ。
ここでは、ウィキペディアの記事にならって、γを回復率としたけれど、隔離率、R(t)を隔離された人口（病気からの回復による免疫保持者、隔離者、死亡者）とするのが正式みたいだケロね。
なお、ここでいう隔離者人口は感染症に発症した人、病院や隔離施設に隔離された人のことではないみたいだから、注意するにゃ。
また、感染者数R(t)は、時刻tの時点において感染症に罹っている人数のことで、時刻tまでに感染症の人数、つまり、累積感染者数ではないので、この点も注意するにゃ。
ところで、
先述したように、γ=0のときI(t)はロジスティック関数になるので、この数値計算の結果の妥当性を調べるために、β=0.025、Δt=0.1ときのオイラー法に数値解とロジスティック関数との比較を以下に示す。
ロジスティック曲線は、tがあまり大きくないとき、指数関数的な挙動をするので、Δt=0.1という粗い計算だと、どうしても誤差が大きくなってしまうけど、定性的にはロジスティック曲線を再現しているにゃ。
Δt=0.05くらいにすると
と、結構、合うようになるにゃ。
ちなみに、オイラー法を用いた計算法は、
とし、連立微分方程式をつぎの漸化式に近似し、k=1,2,3,・・・に対して、逐次計算すればいいにゃ。
計算精度はともかく、オイラー法ならば、プログラムを作ることなく、表計算ソフトで、簡単に、非線形の連立常微分方程式を解くことができるにゃ。
そんな物好きがいるかどうかわからないけれど、計算に使用したスプレッドシートをネットで公開したので、興味のある奴はダウンロードし、β、ガンマの値などを変化させたとき、計算結果がどうなるか確かめてみるといいにゃ。
アドレスは
オイラー法を用いた数値解は、何故か知らないけれど、
を満たしいて、定性的に、そしてある程度定量的に（１）の解の挙動を知るにはこれで十分でしょう。
より高精度な計算をしたい奴は、Δtを小さくするのではなく、２次や４次のルンゲ・クッタ法を用いて解くといい。
精度には４次のルンゲ・クッタ法を用いたほうがいいけれど、この程度の微分方程式ならば、計算精度的に２次のルンゲ・クッタ法で十分だと思う。