感染症数理モデルで面白いものを見つけたので紹介します。

こんにちは。寝起き(@teimon0)です。

これはョョョAdvent Calender 2020の6日目の記事です。
adventar.org



もともとョョョAdCには卍関数しか書くつもりはなかったのですが、なんか空いているのが目に入ったので、とりあえず場繋ぎ的にこの記事を書いています。

22時45分に執筆を始め、慣れないはてなブログに悪戦苦闘しながら、もうすぐ3時にもなろうとしていますがやっと書き上げました。(なお明日は午前6時起きである。絶起の未来ィ、、、)(追記12/6 06:30|絶起回避です。)

ョョョねこ要素はありませんが、この新型コロナウイルス騒ぎと若干関係がある(ほぼない)ことを書きますので、まあ適当に読んでください。

 

1.数理モデリングってなーに

数理モデリングとは、現象を数理モデルを用いて記述することを指しますが、とくに時間変化を見る時に用いられます。言ってることは面倒くさそうですが、要するに物理法則みたいに時間の連立微分方程式を立てたら、現象を説明できるんじゃない???ということです。
たとえば、物理法則の中で、高校生でも習う式としていわゆる運動方程式ニュートンの第二法則)があります。
これは、さまざまな書き方がありますが、次のように書くことができます。
\displaystyle{\frac{d\overrightarrow{v(t)}}{dt}=\frac{\overrightarrow{F(t)}}{m}}

これは、当然ですが、ある瞬間の与えている力に比例して速度が変化する、ことを意味します。

たとえばある物体について見てみると、\displaystyle{ \overrightarrow{F(t)}=c \overrightarrow{v(t)}}となっていたとしましょう。ただし\displaystyle{c}は定数とします。すると、先述の運動方程式から、この物体の速度は、その瞬間の速度に比例して変化するんだなあ、ということがわかります。

これを、今回議論したい「数理モデル」に落とし込んでみましょう。たとえば\displaystyle{v(t)}のかわりに、単位時間あたりの分裂率\displaystyle{\alpha}である(つまり、ある瞬間に100個のバクテリアのうち、一定時間経過後に100\displaystyle{\alpha}個のバクテリアが分裂しているということ。)バクテリアの数\displaystyle{n(t)}とおけば、ある瞬間のバクテリアの個体数\displaystyle{n(t)}は次の式で書き表すことができます。

\displaystyle{\frac{dn(t)}{dt}=\alpha n(t)}

この式はいま運動の式から適当な代入によって導くことができましたが、これをもっと原理的に導くことも可能です。
微小時間\displaystyle{\Delta t}を考えます。単位時間あたりの増加率は\displaystyle{\alpha}でしたから、ある瞬間に\displaystyle{n(t)}個あったバクテリアから生み出されるバクテリアの数(新しい個体の数)は\displaystyle{\alpha n(t) \Delta t}です。ここで、\displaystyle{\Delta t}はごく短い時間ですから、この間に「分裂によってできた個体が再度分裂して個体が増える」数は無視できるほど少ないことに注意してください。

すると、\displaystyle{t+\Delta t}における個体数\displaystyle{n(t+\Delta t)}は以下のように書くことができます。
\displaystyle{ n(t+\Delta t)= n(t)+ \alpha n(t) \Delta t}
右辺第1項を移項し、全体を\displaystyle{\Delta t}で割れば、
\displaystyle{\frac{n(t+\Delta t)-n(t)}{\Delta t}=\alpha n(t)}
となり、高校で数学IIをやった方ならなんとなく見知った式が出てきます。

もう一歩進んでみましょう。\displaystyle{\Delta t}はごくごく短い微小時間を指していますから、\displaystyle{\Delta t\rightarrow 0}の極限をとれば、右辺は\displaystyle{\Delta t}を含まず、その影響を受けませんから
\displaystyle{\lim_{\Delta t \to 0} \frac{n(t+\Delta t)-n(t)}{\Delta t}=\alpha n(t)}
と書くことができます。

そう、これは微分の定義式にほかなりませんから、結局右辺を微分形式に書き換えることで
\displaystyle{\frac{dn(t)}{dt}=\alpha n(t)}
が得られることになります。

この式は解析的に解くことができます(\displaystyle{n(t)=n(0){\rm exp}(\alpha t)}になります。いわゆる指数関数的増加ってやつですね)が、必ずしも微分方程式は解析的に解くことはできません。しかし、それは厳密解を得られないというだけで、今の時代スパコンで殴れば、いい感じのシミュレーションができますから、ここに現実世界の現象を数理モデルで表現する旨味があるわけですね。

2.SIRモデルってなーに

生物学者のKermackとMcKendrickは1927年に感染症の流行を説明する数理モデルを発表しました*1。これは感染症の急速かつ短期的な流行をよく表すことが知られています。議論の対象となる人口をSusceptible(健常)、Infected(感染)、Recovered(回復)の3クラスに分けたことから、SIRモデルと呼ばれます。

この数理モデルは、以下の3式より成ります。
\displaystyle{\begin{align}
\frac{dS}{dt}=&-\beta SI\\
\frac{dI}{dt}=&\beta SI - \gamma I\\
\frac{dR}{dt}=&\gamma I
\end{align}}

それぞれの式を見ていきましょう。

まず、第1式は\displaystyle{S}の時間変化を表示します。どう時間変化するかというと、感染者\displaystyle{I}と遭遇したら、確率\displaystyle{\beta}で感染し、\displaystyle{I}へ移行します。感染者との遭遇の可能性を\displaystyle{SI}であらわしています。これはよくわからないかと思いますが、私見では伝染リスクのある遭遇の回数を示しています。

 

第2式は\displaystyle{I}の時間変化を表示します。第1項は\displaystyle{S}からの流入を表し、第2項では感染者\displaystyle{I}が確率\displaystyle{\gamma}で回復、または死亡などにより、感染症から隔離され、\displaystyle{R}に移行することを示します。単位時間当たり確率\displaystyle{\gamma}感染症から離脱するということは、感染継続期間が平均すると単位時間の\displaystyle{\frac{1}{\gamma}}倍であることがわかります。つまり、\displaystyle{\gamma}感染症の症例持続時間という意味も持っています。

第3式は\displaystyle{R}の時間変化を表示します。これは\displaystyle{I}からの流入を示しているにほかありません。

ここで、適当な代入によって、この数理モデルをもっと簡単にすることを試みます。

まず、人数を表す\displaystyle{S,I,R}をいずれも総人口\displaystyle{N}あたりの数\displaystyle{s=\frac{S}{N},i=\frac{I}{N},r=\frac{R}{N}}とします。ここで人口が系外に出たり、あるいは系外から入ってくることはなく、\displaystyle{N=S+I+R\;(={\rm const.})}です。時間もせっかく\displaystyle{\gamma}が感染継続期間の意味を持っているということがわかったので、これを基準に、\displaystyle{\tau=\gamma t}とおきましょう。

これらを第1式、第2式に代入してみると、次のような方程式が得られます。

\displaystyle{\begin{align}\frac{ds}{d\tau}=&-\frac{\beta N}{\gamma}si\\\frac{ds}{d\tau}=&\left(\frac{\beta N}{\gamma}s-1\right)i\end{align}}

第2式に注目してください。これは感染者の増減ペースをあらわしていますが、\displaystyle{i}はゼロを下回ることはありませんから、\displaystyle{\frac{\beta N}{\gamma}s-1}が正か負かで感染者が増加するのか減少するのかが決まります。

はじめは当然感染者が増加しますが、感染が進み、ある時点で

\displaystyle{\frac{\beta N}{\gamma}s-1=0}

となり、それ以降はこの項は負の値をとることになります。この境界点の条件を見てみれば、\displaystyle{s+i+r=1}を用いて

\displaystyle{i+r=1-\frac{\gamma}{\beta N}}

が得られます。

この式は、導出の経緯を鑑みれば、「総感染者数\displaystyle{i+r}\displaystyle{1-\frac{\gamma}{\beta N}}を超えると、感染者数は減少に転じる」ことを示します。これがいわゆる集団免疫です。つまり、感染拡大という観点では、\displaystyle{\frac{\beta N}{\gamma}}が1よりも大きいか、小さいかが非常に大きな影響をもっていることがわかります。

この\displaystyle{\frac{\beta N}{\gamma}}は表式からわかる通り、1人の感染者が感染期間中に\displaystyle{N}人の系内で伝染させる人数を表しています。このことから、\displaystyle{\frac{\beta N}{\gamma}}を「基本再生産数\displaystyle{(R_0)}」と呼びます。

3.ほかの感染症モデル

先述のSIRモデルは、系の内外での人口の出入りを無視できるような非常に短期間で、かつ急速に進行するような短期間のモデルとしては単純かつ優秀ですが、場合によってはこれを修正したモデルが必要になります。

そこで、たとえばいわゆる潜伏期間にある感染者を\displaystyle{E}(Exposed)とし、Infectedとは別に扱うSEIRモデルや、回復者が免疫を獲得しないとし、回復者はRecoveredではなくSusceptibleに戻るとするSISモデル、性感染症を扱うため男女それぞれのSuceptibleとInfectedが出現し、異性のSuceptibleにしか伝染させることのないとするSI並列モデル、そのほか死亡率が無視できない感染症についてDiedを導入することものや、季節による感染率の変化、当該地域の人口の年齢層などを加味するものなど、さまざまなものが考案されています。

4.2014年の西アフリカ地域におけるエボラ出血熱の流行と数理モデル

2014年に西アフリカ地域でエボラ出血熱が流行したというニュースは、6年も前とはいえある程度覚えている方もいらっしゃると思います。WHOが2014年8月11日に発表したデータ*2によれば、感染の疑いも含めて把握された1975件のうち、死亡例は実に1069件にのぼるという、大変なものでした。Wikipediaにも独立した記事がありますから、詳細はそちらをご確認ください。

さて、感染症の流行予測にあたって、途中のデータ(要するに、何月何日のデータ、とか)をうまく取り出すことができれば、それを根拠に数理モデルの計算の一助にしたり、あるいは現在検討している数理モデルの正当性をシミュレーション結果と実測値を重ねることで確かめることができます。

日本は先進国ですから、こういう感染症の問題が発生した場合、衛生当局もガンガン動いていますし、健康保険制度のおかげで病院に気軽に(気軽に?)かかることができます。そのため、現在の感染者数あるいは死亡者数として出てくる値の数え落としのリスクはある程度低く抑えられています。

しかし、本件のような西アフリカ地域ではどうでしょうか。公衆衛生は整っておらず、医療も普及していない地域が多く、教育も十分に提供されていない可能性があり、衛生当局どころか政情不安定のため内戦や無政府に陥っている地域もあるような状況です。

このような地域で発生した感染症について、たとえばこの感染症を発症した場合も、当局等に発見されない可能性があります。発見されないまま死に至った場合、それが当局等に発見されないまま、公表値に出現しない可能性があります。そのため、当局等が公開した感染者数や死亡者数を先述のSIRモデルやSEIRモデルにそのまま代入することはのぞましくなく、発見された発症者と発見されていない発症者、また確認された死亡者と認知外の死亡者は別の変数として扱うべきです。

これに対処するため、山本ら*3は次のようなモデルを提案しました。

\displaystyle{\begin{align}\frac{dS}{dt}&=-\beta S \frac{I_f}{S+E+I_f+R}\\\frac{dE}{dt}&=\beta S \frac{I_f}{S+E+I_f+R} -\kappa E\\\frac{dI_f}{dt}&=\kappa E- \gamma I_f-\tau I_f\\\frac{dI_h}{dt}&=\tau I_f-\gamma I_h\\\frac{dR}{dt}&=\gamma(1-\delta)I_f+\gamma(1-\mu\delta)I_f\\\frac{dD_d}{dt}&=\gamma\delta\nu I_f+\gamma\mu\delta I_h\\\frac{dD_n}{dt}&=\gamma\delta(1-\nu)I_f\end{align}}

 ここで、\displaystyle{S}はSusceptible(未感染)、\displaystyle{E}はExposed(潜伏)、\displaystyle{I_f}はInfected(発症)のうちfree(未隔離)、\displaystyle{I_h}はInfectedのうちhospitalized(隔離済)、\displaystyle{R}はRecovered(回復)、\displaystyle{D_d}はDied(死亡)のうちdetected(認知)、\displaystyle{D_n}はDiedのうちnon-detected(未認知)を指します。実際に当局等によって発表される数値は\displaystyle{I_h}\displaystyle{D_d}のみです。

また、係数については、\displaystyle{\beta}が伝達係数(\displaystyle{I_f}接触したときの感染確率)、\displaystyle{\kappa}\displaystyle{E}が発症し\displaystyle{I_f}へ変化する時間当たりの遷移率、\displaystyle{\gamma}\displaystyle{I_f}および\displaystyle{I_h}が最終的なステージ(\displaystyle{R,D_d,D_n}のいずれか)に変化する時間当たりの遷移率、\displaystyle{\tau}\displaystyle{I_f}が隔離され\displaystyle{I_h}に変化する時間当たりの遷移率(といっても当局等に発見されながら隔離されない発症者がいるとも思えないので、実質的には当局等によって発見される割合と解釈してよい)、\displaystyle{\delta}は発症者の時間当たり死亡率、、また隔離環境で医療機関や当局による治療を受けられる場合は死亡率がそうでない場合より下がると考えられるので\displaystyle{\mu}はその低減率、\displaystyle{\nu}は生前に当局等に発見されず死に至ったケースのうち死亡の届出などによって当局等に認知されるに至った割合をあらわします。

このモデルを用いて、適当な係数を設定しシミュレーションをした結果と、実際の発表データが論文内で比較されているので引用します。

f:id:neoki_teimon0:20201206023004j:plain
原論文:図4

非常によく実際値とシミュレーションが一致していることが確かめられます。

5.さいごに

なぜ4.でこれを取り上げたかというと、「見えない数字がある?見えないことにしちゃえ!!」という仮定が非常に面白いと感じたためです。数理モデリング、式立てるだけならたのしそうですが、、、

あと、数学は専門じゃないので、間違ったことを言っていたらごめんなさい。

ョョョAdCも1週間が経過しようとしていますが、毎日読んでは「うにゃーーーーわかんねーーー面白いこと書いてるんだろうなーーーすげーーー」と言っています。僕自身ももう一通書く予定ですので、その際はまた読んでいただけると非常にうれしいです。

あと、はてなブログで太字のベクトルをどうやったら書けるか教えてください、、、

ではでは。

*1:W. O. Kermack and A. G. McKendrick (1927). “A Contribution to the Mathematical Theory of Epidemics”. Proc. Roy. Soc. of London. Series A 115 (772): 700-721

*2:https://www.who.int/csr/don/2014_08_13_ebola/en/

*3:山本健久, 早山陽子, 筒井俊之「数理モデルを用いたエボラウイルス感染症の流行の解析」, 獣医疫学雑誌,19(1),p48-54(2015)