スキップしてメイン コンテンツに移動

感染症モデル入門 (1) 難しい！けど安心して


　感染症モデル入門： (1) SIRモデル（連続型）の場合

１．連立非線形常微分方程式で表現したSIRモデル

　8割おじさんで有名になった接触率「8割削減」は、その計算の基礎を100年以上も前の「SIRモデル」によっている(Kermack and  McKendrick (1927))。第１回目は連立微分方程式で記述したSIRモデルについて、それが何を意味しているか考えてみる。SIRモデルは次のような3本の微分方程式から成り立っている。

\begin{eqnarray}
\mbox{感染可能人数} \qquad  \frac{d{S_t}}{dt} &=&  - \beta {S_t}{I_t} \hfill \\
\mbox{感染した人数} \qquad \frac{dI_t}{dt} &=& \beta {S_t}{I_t} - \gamma {I_t} \hfill \\
\mbox{回復した人数} \qquad \frac{dR_t}{dt} &=& \gamma {I_t} \hfill \\
\end{eqnarray}

　ここで、$S_t$, $I_t$, $R_t$はそれぞれ時間$t$とともに「連続的」に変化するみっつの変数であり、$\beta$と$\gamma$は時間が経過しても変わらない一定の値をとる定数（パラメータ）である。$S_t$は感染症の専門書や論文では「感受性人口」と呼ばれる。これでは馴染みが無い学術用語なので、感染する可能性のある人の数、つまり「$t$期の感染可能人数」と呼ぶことにする。 $I_t$は、感染可能人数($S_t$)のなかから、発生した$t$期の感染者数である。パラメータ$\beta$は感染率をしめす。$R_t$は感染した人$I_t$のなかから回復した$t$期の人数である。右辺のパラメータ$\gamma$は回復率示している。(SIRモデルでは、全ての感染者は、回復すると考える。ハッピーなケースから始めることにしよう。死者が出るというようなことを議論すると気が滅入りますから！)

　感染可能人数を英語ではSusceptibleと呼び、感染者をInfected、回復者をRecoveryと呼ぶので、それぞれの頭文字をとって、この３本の方程式で表される感染症モデルをSIRモデルと呼ぶ。

　このモデルは「連立非線形常微分方程式」である。多くの人にとっては一生でおめにかかることとのないものであろう。しかし、今度のパンデミックリスクではこのモデルの結果が日常の生活に大きな影響をあたえたのだ。
　第1に「連立」の意味は、３本の方程式からなりたっていることと、それぞれが互いに関係していることを意味する。第2に「常微分方程式」というのは、この連立方程式で説明される左辺にくる３つの変数$S_t$, $I_t$, $R_t$が時間に関する微分（偏微分でなくて）の形をとっているからである。3つの式は、時間が微小経過した時に、これら３つの変数$S_t$, $I_t$, $R_t$がどのくらい変化するかをしめしている。第3に「非線形」であるとは、式(1)と(2)の右辺に、２つの変数$S_t$と$I_t$を掛け算$S_t \cdots I_t$した項があるからである。
　モデルが非線形であるがゆえに、このSIRモデルは、閉じた解を持たない。閉じた解というのは、右辺がわかれば左辺の任意の時点の値が計算できる形の解が得られるということである。$S_t$, $I_t$, $R_t$の時間経路をもとめるためには、数値解法に頼らざるを得ない。次回以降でExcelでこの「連立非線形常微分方程式」を数値シミュレーションによって解いた結果を示す。

２．コンパートメントで表現したSIRモデル

　上の式１から式３までは、次のようなコンパートメント(Compartment)あるいは変換図（Transfer diagram）と呼ばれる図によって、わかりやすく表現できる。コンパートメントというのは「区切り」という意味で、オリエンタル急行列車での個室のようなものだと思ったらよい。
 それぞれの個室は上の３つの変数$S_t$, $I_t$, $R_t$を表しており、矢印は３つの変数間の「関係」つまり関数式を示している。




次のような実際にあった事例を考えてみよう。東京都に「御蔵島（みくらじま）」という小さな島がある。人口は2020年6月現在321人である。４月に都心に出張をしていた村議会議長が帰島後の20年5月7日に新型コロナウィルスに感染していることがわかった。
　3月7日を時点ゼロ$t=0$としよう。そうするとこの時点の感染可能人口は$S_0=320$人（議長を除く村民の数、321人-1人, 感染人口$I_0=1$人（議長その人）である。もちろん、この時点の回復人数はゼロであるので$I_0=0$である。添え字がゼロであるので、これらを３つの変数の初期値($S_0$, $I_0$, $R_0$)と呼ぶ。
　もし何も対策を講じないと、これ以降$t=1,2,3,\cdots$に、感染人口、回復人数、そして感染可能人口はどのように推移していくのかを示すのがSIRモデルである。ただし、微分方程式で、時間の推移は「連続」なので（厳密ではないか直感的には1秒！ごとと考えてよいだろう）その様な人数の推移を考える。

式２から説明を始めよう。式２は感染者数が単位時間（この場合は１秒単位で！）で何人増えるかを示している。感染人数の「純増加分」は、１）右辺第１項の「その時の感染可能人口$S_t$」に感染者数$I_t$と掛けたものに感染率というパラメータを乗じたものから、２）第２項の回復した人数（回復率☓感染者数）、を差し引いたものとして、計算できる。

感染者数の増加分は、式２の右辺の第１項部分でしめされている。従って、その人数分、感染可能な人数は「減少」する。式１はその点を表している($dS_t/dt=－\beta S_t I_t$)。この式１の右辺はマイナスがついているので、感染可能人口が「減少」することを示している。

式３($dR_t/dt=\gamma I_t$)は回復者の増加人数を示している。式２の右辺第２項で示されているように、感染者の内で回復する人数は、回復率☓感染者数($\gamma I_t$)として計算できので式３の回復者数をしめすしきは$dR_t/dt=+\gamma I_t$となる。

ここで重要なことは、なぜ感染者数の「増加」が、感染可能人数と感染者数の「積」の一定割合として表現できるかである。そうせずに、感染者数に一定割合を掛けたものとすれば、つまり、

\[ \frac{dI_t}{d_t} = \beta I_t - \gamma I_t \]

とすれば、このときのSIRモデルは、線形微分方程式になり、しかも、実際には3つの常備分方程式を連立させずに、一つ一つ解いていけばよい。簡単に閉じた解を得ることができる。そうしないのは重要な理由があるからである。そのことは次回このモデルを差分方程式として再構成したモデルで議論することにする。

次回では、一秒ごとの変化(連続)でなく、１日ごとの変化(離散)を示す差分方程式になおして説明する。えー「差分方程式」なって知らない、という人も多いかと思うが、心配しないでほしい。高校数学IIIの数列を知っていれば大丈夫。「数列」なって忘れた！という人も心配ありません。数列の難しいことをしらずとも、SIRの本質がわかるように説明します。また、Excel（あるいはネット上でつかえる無料のGoogleスプレッドシート）で、いろいろな状況を考えてシミュレーション分析をしてみることにします。安心してください。Excelのシートも公開することにします。

ここで示した常微分方程式の理解は、経済学部で学んだことのあるひとならそれほど難しいことはない。経済学部では、経済数学とか経済動学(economic dynamics)を学ぶ。経済数学は必修科目であろう。もし忘れた！というようなひとがあれば、チャンの有名なテキスト(Chiang and Wainwright (2005))、邦訳では下巻の第5部を読み返してください。

なにか質問があれば、下記のコメントにかきこんでください。

参考文献
  1. Kermack, W. O., and  McKendrick, A. G. (1927). A contribution to the mathematical theory of epidemics. Proceedings of the royal society of London. Series A, Containing papers of a mathematical and physical character, 115(772), 700-721
  2. Chiang, Alpha C., and Kevin Wainwright (2005). Fundamental Methods of Mathematical Economics, Fourth Edition, Boston, McGraw-Hill/Irwin (『現代経済学の数学基礎(下）』, A.C.チャン、K.ウエインライト著、 小田 正雄、高森 寛、森崎 初男、森平 爽一郎　訳)の第5部を参照のこと
  3. デヴィッド・バージェス,モラグ・ボリー. (1990).『微分方程式で数学モデルを作ろう.』日本評論社　（Kermack and  McKendrick(1927)の説明、微分方程式によるSRIモデルの定式化）
ラベル:

コメント

コメントを投稿

ここに感想をお願いいたします。

このブログの人気の投稿

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

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

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

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

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

再生産数$R_t$を理解する。 １．再生産数とはなにか？　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} 　最後の式１の左辺の「感染者数の伸び率」が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$としている。この意味について以下で議論する。その前に感染症...
コメントを投稿
続きを読む