Processing math: 0%
スキップしてメイン コンテンツに移動

感染症モデル入門: (2) SIRモデル(離散型)の場合。易しいです


感染症モデル入門 (2) SIRモデル(離散型)の場合。易しい!


1.SIR感染症モデル:一日ごとの人数変化

 時間間隔が1日という離散的な時間変化を考えると、前回示した微分方程式での記述は、次のような差分方程式に変換できる。
\begin{eqnarray}   \mbox{感染可能人数}  \qquad  {S_{t + 1}} &=& {S_t} - \beta {S_t}{I_t} \hfill \\   \mbox{感染者数} \qquad  {I_{t + 1}} &=&  {I_t} + \beta {S_t}{I_t} - \gamma {I_t} \hfill \\   \mbox{回復人数} \qquad {R_{t + 1}} &=&  {R_t} + \gamma {I_t} \end{eqnarray}

これは連続的な時間の変化dtを離散的な時間変化\Delta tでおきかえたものである。dtが「1秒間隔」であったのに対し、\Delta tは1日単位と考えてみよう。実際のCOVID-19のデータは1日単位で発表されるので、そのように考えたほうが実際に合っている。
 時間経過を1日と考えたのであるから、感染可能人数をdS_t \approx \Delta S_t、感染者数をdI_t \approx \Delta I_t, 回復者数をdR_t \approx \Delta R_tと表現する。これらも1日あたりの人数の変化である。
 そうすると、例えば、感染可能人数を示す微分方程式は\left( {\frac{{\Delta {S_t}}}{{\Delta t}} = \frac{{\Delta {S_t}}}{1} = \Delta {S_t} = {S_{t + 1}} - {S_t}} \right) =  - \beta {S_t}{I_t}と変換できる。つまり{S_{t + 1}} - {S_t} =  - \beta {S_{t}}{I_{t}}を得る。左辺のS_tを右辺に移項すれば上の式1を得る。これと同じ操作を残りの2つの式に適用すると、上の3つの連立非線形差分方程式を得ることができる。以下ではこれら3つの式が何を意味しているかを考えることにする。

2.感染可能人数を示す式1の意味

\mbox{式(1)  感染可能人数}  \qquad  {S_{t + 1}} = {S_t} - \beta {S_t}{I_t}

時間を表す添字tを今日と考え、上の式を言葉で表現するとつぎのようになる。
明日(t+1)の感染可能人数(S_{t+1}) = 今日(t)の感染可能人数(S_t) - 感染率\beta☓(今日(t)の感染可能人数(S_t)☓今日(t)の感染人数(I_t))
ここで注意すべき点は、右辺の第2項の今日の感染者が感染率☓(今日の感染可能人数☓今日の感染人数)と2つの変数を「掛け算」で計算されていることである。ここが感染症モデルの「キモ」であり、この点を理解することが大事である。微分方程式でなく、差分方程式を用いることにより、この部分の直感的な理解が容易になる。

2.1 \beta {S_t}{I_t}の意味を図で理解する。

変数の単位が人数であったので、2つの変数の積は次の様な図で表現することができる。今感染症が始まったばかりの時点t=0と考えよう。既に感染した人が2人(「い」さんと「ろ」さん)いたとしよう。変数で表すとI_0=2人である。また、感染していないが感染する可能性のある人が(A,B,C,D,E)の5人いたとしよう。つまり、S_0=5人である。掛け算をしているのであるから、それは次のように升目をもつ長方形で表現できる。
ここで升目の数は{I_0}\times{S_0}だけあることになる。この場合、升目の数はS_0I_0=5 \times 2=10個である。これは、感染可能者が感染者に接触する「可能性」つまり「場合の数」を表している。最大接触回数とでも表現しても良いかもしれない。  しかし、これは最大数であって、全ての感染可能者が感染者と出会うわけではない。また。接触したからといって適度な「ソーシャルディスタンス」と取っていれば、感染するわけではない。従って最大の接触回数をあらわす{S_0}{I_0}に感染率あるいは感染確率とも解釈できる\betaを掛け、\beta {S_t}{I_t}として今日の実際の感染者数を計算できるはずである。感染率あるいは感染確率\beta$は勿論推定する必要がある。推定問題は重要デリバティブズ、別途検討することにする。

数値例:今日の感染可能人数$S_0=100$人、今日の感染者数$I_0=1$人、感染()$\beta=1$% =0.01とすれば、明日t=1日目の感染可能人数は
{S_1} = {S_0} - \beta \left( {{S_0}{I_0}} \right) = 100 - 0.01\left( {100 \times 1} \right) = 99と99人になる。


3. 明日の感染者数の計算



\mbox{式(2)  感染者数}  \qquad  I_{t + 1} = I_t + \beta {S_t}{I_t} - \gamma {I_t}

これを言葉で表現すると

明日(t+1)の感染者数(I_{t+1}) = 今日(t)の感染者数(I_t) + 今日(t)感染者数今日(t)新規回復者数
= 今日(t)の感染者数(I_t) + 感染率(\beta)☓(今日(t)の感染可能人数(S_t)☓今日(t)の感染人数(I_t)-回復率(\gamma)今日(t)の感染人数(I_t)

数値例:今日の感染人数$I_0=1$人、感染()$\beta=1$% =0.01、回復確率\gamma=0.1=10パーセントとすれば、明日t=1日目の感染者人数は
{I_1} = {I_0} +\beta \left( {{S_0}{I_0}} \right) - \gamma I_0= 1- 0.01\left( {100 \times 1} \right) -0.1 \times 1= 1.9人になる。

3.1  水道方式で説明する。

ちょっと視点を変えて、初等数学の教育方法に革新をもたらした言われている遠山啓先生による「水道方式」を用いて説明をしてみよう。遠山先生は「関数」の概念を水道の蛇口を回すことで流れる水の量の変化で図式化して説明しようとした。SRIモデルは、こうすれば小学校高学年あるいは中学生にでも理解できるとおもう。これと次回以降で説明するEXCELプログラムを用いたシミュレーションによって感できるだろう。

3.1.1 5月7日(時点t=ゼロ)の感染者数と感染可能人数:初期値を設定する。

 時点ゼロでの感染可能人数と感染者との間の関係は、次のように図式化できる。前回のべた御蔵島の例をとって考えてみる。村議会の議長さん1人が帰島後感染をした。5月7日を時点ゼロとすると、感染者数は1人である。残りの島民320人は感染の可能性のある時点ゼロの「感染可能人数」である。上の水槽の水の量が5月7日の感染可能な島民の数、320人を表し、その下の水槽の水が5月7日の感染者(村議会議長)数1人を示している。一番右の水槽の水(回復者数)は最初はだれもいない、つまりゼロ人である。
この図で水道の蛇口はひねられていないので、水(人の移動)は存在し無いことに注意されたい。

3.1.2 5月8日(時点t=1)の感染可能人数と感染者数



 次の日(5月8日)になると上の水槽の水(感染可能人数)のうちで感染した人が発生する。それは\beta {S_0}{I_0}人と計算できる。感染率が1パーセンとすれば新規の感染者数は$\beta {S_0}{I_0}=0.01(300)(1)=3人となる。その分の水(感染可能人数)は上の水槽から下の水槽に「新規の」感染者数として流れ込む。水筒の蛇口のひねり具合が感染率\beta$を表すと言ってよいだろう。大きく(小さく)回せば感染者数が大きく(小さく)増えることになる。

蛇口をひねって流れだした水(新規感染者)は下の水槽に流れ込んで、既にある水(既に感染をした人の数、つまり議長の1人)に付け足される。新規感染者数は最初の式1の右辺で既に計算されているので、それを加えればよい。このことが式2の左辺と右辺の第1項で示されている。右辺の第二項は回復者数を示し、感染者数からの脱落人数に相当する。それは式3と関係するので、次の図で説明をする。

回復者が発生したことによる感染者数の減少とと新規回復者の増加(式2と式3)
今日の感染者I_0人の中から回復者が回復率\gammaパーセントを掛けた人数だけ発生をする。従ってその分感染可能者が減少し、新規の感染者数が増加する。言い換えれば感染者数を示す水槽の水が\gamma I_0だけ減る。減った水(回復人数)はその下の水槽に流れ込む。従って明日の回復者数は、今日の回復者数にこの新規の回復者数を足しただけ、つまり{R_1} = {R_0}+\boxed{\gamma {I_0}}人だけ増える。例えば、回復率が10パーセントとすれば、明日の回復者数は{R_1} = {R_0}+\gamma {I_0}=0+0.1(3)=0.3人となる。蛇口の回し具合が回復率を表している。

以上の計算が終われば、今日の数値(t=0の初期値)をもとにして、明日の感染可能数、感染者数、そして回復者数が計算できたことになる。これを新たな初期値(各式での右辺の値)にして、明後日の感染可能数、感染者数、そして回復者数が計算できる。これを繰り返すことによる、任意の将来時点tでの感染の流行や収束を予測できる。

SIRモデルを水道方式で考えると、すぐわかることは、3つの水槽の水(人数)は毎日変わるけれども、合計した水の水量は変わらないことである。時間がどのくらいかかるかは、感染率と回復率によるが、最終的には、全ての人が回復し、幸せになるのがこのモデルの特徴である。

次回では、こうした計算をExcelを用いて行う。初期値やパラメータ(感染率と回復率)を様々に変化させることによるシミュレーションを試みる。





コメント

このブログの人気の投稿

「8割おじさん」の数理モデル

米国の週刊誌Newwweek誌の日本語版『ニューズウィーク日本版』の最新の 2020年6月9日号 では「日本モデル」と題した特集号で、北海道大学医学部の西浦教授が、「8割おじさん」の数理モデル」という論文を寄稿している。冊子体の本誌でも読めるし、また楽天マガジン・dマガジンなどのオンラインの雑誌読み放題でも全てのページが配信されている。 ニューズウィークとしては異例のアカデミックな論文である。数式こそ無いものの、変数名や学術用語が散見される異例の論文に成っている。同誌で経済や政治関係の記事に親しんでいる人にとっては、ちょっと難しい論文かもしれない。 しかし、非常に示唆に富んだ論文である。特に、「8割削減」の意味を丁重に説明している。最近、8割削減に関して、そこまでする必要がなかったという批判があちこちから出ているが、これをよんでから批判をすべきであろう。 疫学における感染症モデルは勿論、マクロ経済や資産価格決定モデル(CAPMやブラック=ショールズモデルなど)も、所詮は複雑な現実の抽象化に過ぎない。しかし、よいモデルというのは、抽象の革新をつかみ、リスク管理を行う場合、大きな間違いをしない基準となりうる。 また、この号では、”日本のコロナ対策は過剰だったか”という、西浦教授と國井修(世界エイズ・結核・マラリア対策基金)との対談記事も掲載されている。ともにお医者さんでありかつモデリングを研究している二人の言葉は、経済やファイナンスのモデルを用いて仕事をしている、アナリスト、エコノミスト、研究者にとっても学ぶことが大きいとおもう。 國井氏はこう述べている「・・・私達モデラーはリスク評価についてはアンダーアクト(控えめに言う)よりは、オーバーリアクト(大げさに言う)して話をすべきと、肝に命じながらやってきた・・・・」 このことが経済や金融市場のリスク分析に当てはまるかどうか、人それぞれ異なる意見をもっているかもしれない、しかし、1つの教訓だろう。 ただし、西浦先生の論文は、数式を使わないようにしたために、高校数学を理解している人にとっては、帰って難しくなっているきらいがある。むしろ、数式を使って説明したほうが良いかもしれない。ところが、多くの感染症に関する論文や本では、特に日本語で書かれた本や論文では連立微分方程式をつかって説明している。そうす...

感染症モデル入門 (3) Excelで感染症 SIRモデルをシミュレーションする

感染症流行をExcelでシミュレーションする(ファイルがダウンロード出来ます)。 1.準備:ファイルのダウンロードをする。 SIRモデルのわかりやすい(つまり最も簡単な)Excelシートを作成した。Excelファイルを用いたシミュレーションは以下の二通りの方法で実行できる。 1.1 「Google スプレッドシート」を用いてオンラインでSIRモデルのシミュレーションを行う 。  オンラインでシミュレーションをする場合には、Google スプレッドシートをつかいます。「 ここ 」をクリックしてください。「感染症数理モデル_入門と応用_森平_2020-05-29」というファイルがGoogle スプレッドシート上で開くはずです。次のような画面を見ることができます。 携帯などから利用していて、通信環境などが良くないときには、オフラインでの実行をおすすめします。>ファイル?オフラインで実行する、とすれば、自分の携帯にダウンロードしてあるGoogle スプレッドシート上でこのプログラムが実行できるはずです。 1.2 Excelファイルとしてダウンロードして、自分のPC上でシミュレーションをおこなう。  Excelシートとして使い時は、上で示したGoogle スプレッドシートのメニュー上で >ファイル>ダウンロード>MicroSoft Excel (.xlsx)  としてファイルをExcelプログラムしてダウンロードしてください。ダウンロードしたファイルを実行すると、次のような画面が現れるはずです。 2.初期値とパラメータ値を設定する。 2.1 Excelのシートで説明しよう。1行目は変数名とパラメー名前である。2行目の薄緑色のセルに具体的な数値を与える必要がある。ここでは次のような値を初期値とパラメータ値として設定している。 1)  感染可能人数の初期値:S_0=100人 2)  感染者人数の初期値:I_0=1人 3)  回復人数の初期値:R_0=0人 4)  感染率:\beta=0.01,つまり1パーセント 5)  回復率:\beta=0.1,つまり10パーセント (削減率については今回は計算結果に影響の無いようにしている) これらの...

感染症モデル入門 (5) 再生産数Rtとそれに基づくリスク管理

再生産数R_tを理解する。 1.再生産数とはなにか? SIRモデルから導出する。 離散的なSIR感染症モデルの2番目の式(2)を次のように変形する。 \begin{eqnarray}   {I_{t + 1}} &=& {I_t} + \beta {S_t}{I_t} - \gamma {I_t} \hfill \label{eqtn:1} \\    \Rightarrow \,\,\,\,\,\,{I_{t + 1}} - {I_t} &=& \left( {\beta {S_t} - \gamma } \right){I_t} \hfill \nonumber \\    \Rightarrow \,\,\,\,\,\,\frac{{{I_{t + 1}} - {I_t}}}{{{I_t}}} &=& \left( {\beta {S_t} - \gamma } \right) \end{eqnarray}  最後の式1の左辺の「感染者数の伸び率」が0を超える(下回る)ことは、感染人数I_tが増加(減少)することを意味する。言い換えれば、右辺が、 \begin{gathered}   \beta {S_{t}} - \gamma  > 0 \,\,\,\,\,\,\Rightarrow \,\,\,\,\,\,{S_{t}} > \frac{\gamma }{\beta } \hfill \\ \end{gathered} であれば、感染者数は増加する。 両辺を \gamma/\beta で割ると, 明日の時点における再生産数を得ることができる。 \begin{equation} \boxed{\,\,{{\bar R}_t} \equiv {S_t}\frac{\beta }{\gamma } > 1\,\,} \end{equation}  これを時点tにおける実行再生産数と呼ぶ。回復人数を示すR_tと区別するために、Rの上にバーを付けて再生産数\bar{R}_tとしている。この意味について以下で議論する。その前に感染症...