【感染症モデル入門】フィッティングして遊んでみた♬

今回は、SIRモデルを拡張してSIHRモデルでフィッティングを試みた。
本記事も素人が書いたものであり内容は自己責任で取り扱ってください

前回のSIRモデルでは、S-I-Rがある一定割合で連鎖的に進むため、現在のIからRに遷移する部分の遅れを再現できない。
そこで、実際にはI-RのプロセスでHospitalやHomeなど、一定の隔離状態の後、治癒すると仮定する。
時間発展の微分方程式は以下のとおりとなる。
一応、HとRの状態では感染者を増加させないという仮定である。

{\begin{align}
\frac{dS}{dt} &= -\beta \frac{SI}{N} \\
\frac{dI}{dt} &=  \beta \frac{SI}{N} -\gamma I \\
\frac{dH}{dt} &=  \gamma I-\delta I \\
\frac{dR}{dt} &=  \delta I \\
\end{align} 
}

【参考】感染症の土地勘つけるために、参考を３篇上げておく.最初の記事はすごく面白い（pdf）
・感染症の数理 - 東京大学稲葉寿＠平成 20 年度年第 7 回例会「医療とアクチュアリー」講演
・感染症流行の予測：感染症数理モデルにおける定量的課題@西浦博、稲葉寿
・微分方程式と感染症数理疫学@稲葉寿

やったこと

・日本の状況
・終息する各国の状況
・増加する各国の状況

・日本の状況

まず、ターゲットとなる日本の4／30現在のグラフを示す。
横軸は、days_from_22_Jan_20から開始しており、丁度100日目である。

このグラフの現存感染数（赤プロット）、治癒数(緑プロット)、そして死亡数（青プロット）を片対数で示すと以下のようになる。なんとも美しい平行曲線となっている。つまり、一定割合でI-RやI-Dの遷移が起こっていることが分かる。少し波線となっている部分が感染の急激な拡大などを示していると分かる。さらに断定はできないが、同時にどのあたりの時期にそれぞれの遷移が発生しているかが推測できる。

グラフを見ると、新規感染数の棒グラフを見ると、完全にピークアウトしてきており、この後、現存感染数が減少に転じそうな状況である。
確かに、昨日5／1に全国では入院数がマイナスに転じた。しかし、本日はまたプラスになってしまっている。
この先、きちんと終息に向かうかどうかが最大の関心事である。
結果（条件は後で述べる）は、以下のとおりとなった。
日々の新規感染数のフィッティング結果と予測がday_est赤プロットのようになり、現存感染数、治癒数がそれぞれフィッティング出来た。
なお、今回は死亡数のフィッティングはモデル上放棄している。ただし、死亡数と治癒数はしばらくは似たような推移になると予測できるので治雄数の推移からまだまだ増加することが予測できる。

新規感染数も片対数でプロットすると、昨夜の有識者会議の会見で尾身さんがお話しされていたが、傾きが見える。そしてそれが増加している部分より減少に転じた後の方が傾きが小さいことが分かる。つまり、長期化の可能性がある。大体、一番下のグラフから、あとひと月半すると現存感染数（入院数）が半減するという結果となっている。
これを早めるには、西浦さんのお話にもあったように感染者との接触率を下げることである。上式でいうと、$\beta$を下げることである。もう一つの方法は、$\gamma$を上げて治癒率を上げることである。
実質的に実行再生産数（$\frac{\beta}{\gamma}\frac{S}{N}$）を下げることになる。
この数値が１以下になると、感染伝播は発生できず、感染者は自然消滅する。
まあ、言葉で言えば、まず感染者が減少し、かつ感染してもどんどん治せる状態だと、感染は自然消滅するということである。

実は、Nとしての値が実際の観測値から変更しているが、議論が面倒なので考え方をまとめに記すこととする。
ちなみに、日本の昨夜までの実効再生産数R＝1.03であり、感染率$\gamma (R-1)=7.7e^{-4}$であり、いよいよ新規感染者０が実現まじかという状況である。
なお、基本再生産数5.09となっており、かなり高い値なので無防備な状態の集団に感染者が来ると、感染拡大するスピードは速いと言える。

・終息する各国の状況

同じように各国を見ていく

スイス

一番きれいにフィッティング出来たのはスイス。以下のように終息が見えてきた。
実効再生産数0.21であり、感染率－0.049となっており、新規感染数もピークでは1000人を超えていたが、100人程度に減少してきた。スイスでも増加時と減少時の傾きは増加時の方が4倍程度になっており、減少時の傾きは日本と同程度（1桁／40日）である。

ドイツ

綺麗に終息しつつあるのは、ドイツである。
実効再生産数0.22であり、感染率－0.054となっており、新規感染数もピークでは7000人だったが、現在は1000人程度まで減少した。日本のピーク時よりまだまだ多いのだが、減少率は大きい。
ドイツの問題は、その新規感染数の減少の速度である。日本やスイスと比較すると若干傾きが小さく、かつ裾野が残ってしまっており、新たな感染の火種になりかねない。まだまだ予断を許さない状況と言える。
このグラフを見ると、死亡数曲線の様子は治癒数曲線を平行移動したものと同じだと分かる。

韓国

韓国が優秀な国だと報道されているので見てみよう。
以下のとおりです。
実効再生産数0.04であり、感染率－0.051と終息しています。しかし、新規感染数の棒グラフはほぼひと月間100人から減少しなかったことを物語っています。そして現存感染数と治癒数が交差する頃、約ひと月前に減少に転じてやっと10人程度まで減少しました。つまり学ぶべきことは、この経過は終息が始まってからもまだまだ油断は禁物だという事実です。

タイ

あまり知られていませんが比較的優秀なのがタイです。
覚えている人も多いと思いますが、日本同様発生当初から日本、シンガポールと並んで感染者を競っていました。
この3か国ではタイが最初に終息しました。
実効再生産数0.49であり、感染率－0.055と終息しています。そして、新規感染数の棒グラフの減少率も（1桁／20日）と韓国の最近と同程度の大きな傾きで減少しており、新規感染数10人まで減少しました。
特徴はピーク近傍で平たい台地になってしまったことです。これはピークを抑えたという評価もできますが、ピークを長引かせたと言えなくもありません。内容は何だったのか振り返る必要があるでしょう。

イタリア

なんともゆっくり進行していて、医療崩壊だと騒がれたイタリアもどうにか終息を向かえそうです。
イタリアもピークがだらだら長引いたのでフィッティングは不十分です。
しかし、一応実効再生産数0.13、感染率－0.024と終息しています。ただし、あまりに減少もゆっくりなのでまだまだ十分に終息できたとは言えません。ピーク時新規感染数6500人以上でしたが、現在は2000人程度に減少しました。何より現存感染者数が減少に転じ、10万人位まで減少してきました。
そして、何より治癒数が順調に増加してきており、たぶん医療崩壊は避けられたものと思われます。
※死亡率も世界4位で（参考13.61％28,236人）と非常に高いですが。。。

米国

うーん、こっちに入れていいかどうか悩む国です。
結果は以下のとおりです
つまり、新規感染数がピークで止まってしまっていて、減少していません。
なので、フィッティングは不十分であり、結果は参考で載せる程度です。
一応、実効再生産数0.53、感染率－0.016なので減少するはずですが、現在５回目（ほぼ毎週）のピークを迎えていて、感染拡大を繰り返しているようです。しかも、リニアースケールで見ると、累計治癒数20万人ですが、増加がゆっくりに見えます。はっきり言って、この国の感染が止まらないと世界のコロナも消えないと言ってもいいと思います。あとひと頑張り我慢が必要な時期だと思います。

・増加する各国の状況

やはり、インドとロシアを見てみましょう。

ロシア

なんかすごい拡大期です。
実効再生産数2.85、感染率0.063と計算できました。ただし、ちょっとS0が大きすぎるかもしれません。まだまだピーク前であり、拡大期なのでどこまで感染数が大きくなるか不明な状況なのでこの程度の値としています。それでも、全体としては飽和しつつあるようです。あと数日で新規感染数ピークとなり、約２週間で感染数ピークになりそうです（ただし、ちょっとパラメータ変更すると変化するので、保証はできません）。

インド

次に心配なのは、インドです。
インドは前回少しピークが見え始めそうという時期でした。
実効再生産数1.57、感染率0.030と計算できました。ほぼ新規感染数がピークを迎え、現存感染数もあと1週間でピークになりそうです。治癒数が現存感染数曲線に近づいてきており、ロシアよりは進んでいると見えます。

まとめ

・モデルをSIHRに拡張してフィッティング精度を上げた
・日本他各国のフィッティングを実施して考察した
・日本もあと一息で終息に向かえるところまで来たと言えるが、各国の状況を見るにつけ、最近の東京や北海道の振る舞い、そして何より院内感染の状況は改善しないと問題が大きくなると思われ、予断を許さない状況であると言える。
・終息しつつある各国もそれなりに大変な状況であることがうかがえる

・いよいよ、分布関数の導入や得られた諸量の精度について、統計的な取り扱いを実施したいと思う

おまけ

今回の微分方程式の問題はS0が不明な点である。終息してしまえば、一定の評価ができるが、終息前の状況だと適当なS0を見極めることが難しい。
以下、簡単にS0を変化したときの日本のフィッティング状況を掲載しておく。
さらに、行動変容が実際に効果あるとすると、減少の傾きすなわち微分方程式の各係数は一定とは言えないので、この解析方法はやはりだいたいの傾向を示す程度のものと考えた方がよい。