感染症の流行モデル SIRモデル

感染症の流行モデル　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,・・・に対して、逐次計算すればいいにゃ。

　　

計算精度はともかく、オイラー法ならば、プログラムを作ることなく、表計算ソフトで、簡単に、非線形の連立常微分方程式を解くことができるにゃ。

 

そんな物好きがいるかどうかわからないけれど、計算に使用したスプレッドシートをネットで公開したので、興味のある奴はダウンロードし、β、ガンマの値などを変化させたとき、計算結果がどうなるか確かめてみるといいにゃ。

アドレスは

 

https://docs.google.com/spreadsheets/d/e/2PACX-1vQbAtk0oqiyJWLnIxqSrbIS8kMPoLX4D78IWk6hXsnhjSMryHs-9liX5Y4R68Sq-D6YyYfMGzkxJj2V/pub?output=xlsx

 

 

オイラー法を用いた数値解は、何故か知らないけれど、

　　

を満たしいて、定性的に、そしてある程度定量的に（１）の解の挙動を知るにはこれで十分でしょう。

 

より高精度な計算をしたい奴は、Δtを小さくするのではなく、２次や４次のルンゲ・クッタ法を用いて解くといい。

精度には４次のルンゲ・クッタ法を用いたほうがいいけれど、この程度の微分方程式ならば、計算精度的に２次のルンゲ・クッタ法で十分だと思う。

 

 

コメント 3

　お久しぶりです。いやっ本当に、コロナでダウンかと心配しましたよ(^^)。

　「そう言うなら、オイラに数学記事をよこせ！」と言われそうですが、現在またラプラス型の境界要素法に戻ってまして（流体解析というか水理学）、結果が出るのにしばらく時間がかかりそうなのです。

　結果が出たら、投稿しようかなと思ってます。もちろん、よろしければですが・・・(^^;)。
by ddtddtddt (2020-03-13 20:15) 

お久しぶりでございます。
色々とありまして、新潟にほとんどいないので、ブログを更新できないでおります。

ネタに困っているので、どのような記事でも大歓迎です。
by nemurineko (2020-03-17 18:47) 

上品なブランドの新作が期間限定セール発売上品なルイヴィトン、グッチやエルメスなどのブランドコピーの新作が大量入荷しました。種類が豊富で、勝手に選べます。激安の上に、品質には保証があって、末長くご愛用いただけます。新年とクリスマスを迎えで、期間限定セールが開催中で、早く手に入れましょう。モンクレールコピーの新作入荷人気沸騰なモンクレールコピーの新作が大量入荷しました。ユニークなデザインと優れた機能性を兼ね備えた逸品はいま、割引の活動を進行しています。寒い冬に、一着あれば、十分暖かいです。ぜひこちらでチェックしましょう}}}}}}
https://www.gmt78.com/product/detail/1414.htm

by BrAfspor (2023-05-25 15:37)