DEV Community: Tomoaki Masuda

隠れSIRモデルを試作して、東京都の感染者を予測する

Tomoaki Masuda — Tue, 07 Apr 2020 01:48:32 +0000

実地な値を使って検証するのは避けていたのだけど、現在の状況を見て解説することにします。

隠れ感染者が一定量いる

現在、東京都の新規患者が100人/日を超えています。指数関数的増加（差分が1.0以上となるので複利的に増加する）が理解できていれば、次の増加が100人ではなく、それ以上にとなることは容易に推測ができます。

ただし、単純に「指数関数的増加」とはいえ、いったいどのくらいの増加数なのか？を予測しておくことは大切なことです。

新型コロナウィルスの厄介な点は、

感染率がインフルエンザと同程度に低い（麻疹などよりもかなり低い）
潜伏期間が2週間と長い
感染しても無症状な人が一定量（8割程度）いるらしい

ところです。このため、麻疹のように一気に拡散＆重症化する場合は、発病後すぐに家や病院にこもるために、あまり外側に拡散されることがないとありません。潜伏期間が長いと、いつ罹患したのかがやっかいです。クラスター対策班がクラスターの追跡ができなくなってきているのはこのためです。同時に、感染しても無症状であれば、感染したまま普通の生活を送ることになります。無症状であるために、病院に行くこともなく家で休むこともなく PCR 検査や検体検査を受けることしません。このため、本人の意図とは関係なく、ウィルスが蔓延してしまいます。

進化論的に言えば、これが新型コロナウィルスの「生存戦略」ということになります。

隠れ感染者を推測する

潜伏期間中であったり無症状な感染者を「隠れ感染者」と呼びましょう。この「隠れ感染者」は数値の上ではでてきません。毎日、定期的に発表される新規感染者とは別に扱う必要があります。

では、どうやって「隠れ感染者」を推測するのでしょうか？

感染の因果関係として、

新規感染者が発生する
新規感染者は潜伏期間(2週間)前に何処かで感染している
新規感染者を感染させた元の感染者がいる

ことが条件になります。実際は海外から帰国した感染者もいるので誤差でるのですが、おおまかな東京都の感染者の状況はこの感染の因果関係にあてはめられます。

これを SIR モデルに当てはめてみます。

新規感染者 dIt/dt を発症した日（14日前）に「隠れ新規感染者」として設定します。
隠れ新規を累積して、隠れ累積患者数を計算します。このとき、その日の新規患者を取り除きます。これは新規患者としてカウントして隔離されるためです。
隠れ累計と隠れ新規の関係から、感染率が計算できます。感染率 β は、もとの SIR モデルの感染率や基本再生産数とは異なる定義となるため、このモデル独自のものです。もう少し式変形をして、基本再生産数と比較できるようにしていきます。

既に新規感染者が計測されている日から14日さかのぼったところの感染率βが算出できます。
これをもとにして、予測感染率β*を推測します。
一般に、何らかの環境の操作（外出の自粛など）が行われない限り、予測感染率は前日の感染率を踏襲するはずです。

この予測感染率 β* から、予測隠れ新規患者数を計算します。
予測隠れ新規の数は、14日後に予測新規患者数となり、新規の感染者数が予測できるというモデルです。

Excel で予測する

これを Excel 使って計算してみましょう。

東京都新型コロナウイルス陽性患者発表詳細 - 東京都_新型コロナウイルス陽性患者発表詳細 - 東京都オープンデータカタログサイト

4/6 付けの東京都の患者データを使って予測してみます。

青いセルが予測の部分です。「患者累積」と「隠れ累積」が大きくずれているのは、新規患者をうまく拾えていないのと、発症までの潜伏期間が2週間と長いために潜在的な患者が滞留してしまっているのが理由です。
単純計算いけば、新規患者数100人×14日間 = 1400人以上は滞留していることになります。

2週間ほどさかのぼって、感染率の予測値を入れます

0.160 は、3/22 までの大まかな平均値です。4/6 付けで新規患者数が83人と若干減っているので実効感染率が一時的に下がっていますが、いまのところ外れ値として処理しています。継続して下がっていれば、この感染率をさげます。
4/4 からの 0.080 は土日で自粛が効果的に行われたとみなして下げてあります。

感染率の絶対値は、特に意味はありません。いずれ実効再生産数のような割り出しが必要になります。
ここでは、相対的な値となっています。

それぞれの値をプロットしたの下の図です。

予測患者数累積が非常にあがっていることが分かります。

対数軸にプロットしたとき、患者累積と隠れ累積の差が大きくならないように注意が必要です。

現時点でプロットしてみたところ

4/11 頃に患者数の累積が2,000人を超える
新規患者は500人のピークがある
4/4 の自粛は、4/18 頃に効果が現れるので、それまでは同じペースで自粛が必要となる
増え方は急ではあるが、きちんと自粛をすれば、患者数が1万人を超えるのは4月末以降の予想となる

問題は、新規患者数が500人/日程度発生するが、医療崩壊を起さずにこれを受け入れらるか？ということです。
このモデルで計算したときには、「隠れ累積患者数」が4,000人程度いることと2週間遅れで発病することをあわせると、かなり厳しい状態と言えます。ただし、軽症/重症者が混在するので、この比率が問題になります。

データ

Excel データは以下でダウンロードが可能です。

参考文献

感染症数理モデルでシミュレーションする（その２）

Tomoaki Masuda — Thu, 26 Mar 2020 01:39:10 +0000

3週間ほど経ち状況が変わってきていますが、感染症数理モデルの話を引き続き

SEIR モデルを JS で書く

新型コロナウィルスの専門家委員会から「基本再生産数 R0」や「実効再生産数 R」という言葉でてきます。他にも専門用語がでてきますが、これらを一般的に解説するかどうかの前に「その専門分野ではどのような使い方をされているのか？」を正確に把握するする必要があります。まあ、それが「専門分野」なのですから、同じ言葉を使っていても専門分野ごとに多少の違いはあるし、逆に同じ用語を他の分野でも同じ様に使っている場合もあります。どちらにせよ、その分野での専門用語の把握が必要です。

そんな訳で、２週間ほど前になりますが、SEIR モデルを Javascript で書きました。元の Python からそのまま JS に書き替えて、そのあと R0 を取り込めるようにしています。ちょっと、R0 の意味が分からなくてコードが間違えていたのですが、先日直したので、ここにコードの解説をしておきます。

SEIR モデルの数式

実は、Wikipedia を見ると、日本語のものと英語のものと SEIR モデルの式が若干違います。感染率 β の部分で総数 N が出てくるのですが、日本語の場合は、β(N) のように総数に依存するのですが、英語版のほうは β は N から独立しています。
感染率が総数 N に依存するのは変な感じがするので、英語版のほうを使うことにしました。

The SEIR model

出生率 Λ と死亡率 μ は短期間なので、0 とみなせます。これを書き方えたのが下の式になります。

N: 全体の総数
S: 未感染者
E: 潜伏期間中
I: 発症者
R: 回復者/免疫を得た人

になります。イギリスで発表された全国民が感染して免疫を得るという方式は、この SEIR モデルの、R になるということです。SEIR モデルの場合、感染して潜伏→発症→回復という順番になるのですが、かならず回復するというモデルなので、今回のような新型コロナウィルスのように発症後に死亡するケースがある場合には、致死率が2％程度だとしても相当の死亡者が予想されます。

即効撤回されたのは、良い判断でした。

この式の中でのパラメータとして、以下の3つがあります。

α: 潜伏率
β: 感染率
γ: 回復率

これらは、潜伏期間 lp と発症期間 ip で書き替えることができます。

α = 1/lp : 潜伏帰還の逆数
γ = 1/ip : 発症期間の逆数

肝心の感染率 β の値ですが、基本再生産数 R0 との関係が以下のようになります。

死亡率 μ が 0 となるので、R0 は感染率 β と回復率 γ との比になります。
式の中では感染率 β を使うことになるので、次のように書き替えます。

この式から分かることは、

感染率 β と回復率 γ が等しいときに R0 が 1.0 になる
感染率 β が回復率 γ が大きいときに R0 は 1 を超えて、感染が拡大する
感染率 β が回復率 γ が小さいときに R0 は 1 より小さくなり、感染は縮小していく

という現象になります。これが専門家会議の言う「基本再生産数 R0」の意味です。
なので、できるだけ R0 < 1.0 の状態を保ちながら、ピークを低く抑えましょうという対策の根拠になります。

SEIR モデルを JS に直す

この数式を Javascript に直すと次のようになります。
v のところが配列になっているのは Python コードの名残りです。

function seir_eq(v,t,alpha,beta,gamma,N) {

    S = v[0]
    E = v[1]
    I = v[2]
    R = v[3]
    ds = - beta * I / N * S             // dS/dt = -βI/N*S
    de = beta * I / N * S - alpha * E   // dE/dt = βI/N*S-αE
    di = alpha * E - gamma * I          // dI/dt = αE - γI
    dr = gamma * I                      // dR/dt = γI 

    return [ds,de,di,dr];
}

これを100日間繰り返すのが次の calc 関数です。
SEIR モデルに、試しに隔離率(T) を埋め込んでみて、発症者の一定割合を病院に隔離するというシミュレーションも入れておきます。
あと、S,E,I,R の数はマイナス値にはならないので繰り返し計算をするときに補正します。

function calc(state,alpha,beta,gamma) {
    var t_max = 100 ;
    var dt = 1 ;
    lst = []
    var N = Sinit + Einit + Iinit + Rinit
    console.log( state );

    for ( var i=0; i<t_max; i++ ) {

        var d =  seir_eq( state, i, alpha,beta,gamma, N )
        var Si = state[0]+d[0]
        var Ei = state[1]+d[1]
        var Ii = state[2]+d[2]
        var Ri = state[3]+d[3]
        // 感染者を発見して隔離する
        dx = Ii * T
        Ii = Ii - dx
        Ri = Ri + dx // 免疫者に加算

        // マイナス値を調節する
        if ( Si < 0 ) {
            Ei = Ei + Si; Si = 0;
        }
        if ( Ei < 0 ) {
            Ii = Ii + Ei; Ei = 0;
        }
        if ( Ii < 0 ) {
            Ri = Ri + Ii; Ii = 0;
        }

        state = [ Si, Ei, Ii, Ri ]
        // console.log( state );
        lst.push( state );
    }
}

グラフ化する

これに vue.js とグラフツールの c3.js を追加したが次の形です。

SEIR モデルシミュレータ

このグラフは基本再生産数 R0 が 10 のときのものです。10 というのは非常に多い値なのですが、いわゆるダイヤモンドプリンセス号のような閉鎖空間のクラスター（患者集団）の場合がこれに相当します。時間を追うごとに加速度的に感染者が広がります。

注目して欲しいのは、以下の３点です。

発病者のピークは感染スタートよりも左にある
潜伏期間中（黄色）のピークは、発病者（緑）のピークよりも左にある
最終的にほとんどの人が回復者（赤）となる

このグラフでは潜伏期間２週間で考えているので、発病者（緑）のピークは感染スタートから２週間以降になります。感染スタート時にすべての人が感染する訳ではないので、幅を取って２週間から３週間というところです。これが、現在の花見や春休みのシーズがスタートとなったときに、2,3週間後に流行するのでは？という理由です。

当然のことながら、発病期間のピークの前に、潜伏期間のピークがあります。SEIR モデルは潜伏期間中には罹患させないモデルになっていますが、新型コロナウィルスは潜伏期間中も罹患させる可能性がある、または無症状の場合も感染さえる場合があるという「可能性」が指摘されています。
このため、潜伏期間中の人が活発に活動すると、その後に発症のピークが出るかもしれません。

SEIR モデルの場合は、最終的にほとんどの人が感染して回復します。というか、無限に計算を続ければ100%の人が感染し回復して終わるようになっています。
先に書きましたが、SEIR モデルの場合は発症→回復となるので、死亡者がいません。実際の新型コロナウィルスでは死亡するので、最終的な回復者数×死亡率が全体の死亡数になってしまいます。
このため、無防備な策は誤りです。

基本再生産数 R0 を 1 に近づける

基本再生産数 R0（あるいは実効再生産数 R）は、感染率 β と回復率 γ との比率になります。

R0 を 1.0 に近づける（1.0以下にする）というのはどういうことかというと、数式からは以下の2つの手段があります。

感染率 β を下げる
回復率 γ を上げる

感染率を低く押さえるというとは、いわゆる「誰かに罹患させる確率を減らす」ということです。移動を減らすとか、3つの条件が重なる場所を作らないあるいは行かない、というのがその手段です。
もうひとつは、回復率を上げるというのは、病床を増やす、体力を蓄えておいて軽症で済ませるなどがあります。

基本再生産数 R0 を 10,5,2 と変化させてみましょう。

発症者（緑）のピークが限りなく右にずれることがわかります。
つまり「ピークを平らにする」ということは、基本再生産数 R0 を限りなく1.0に近づける（あるいは1.0以下にする）ことであり、同時に感染者のピークを限りなく未来に追いやることができます。

そして、新型コロナウィルスに対しての有効なワクチンや予防接種ができるまでの時間稼ぎができます。

基本再生産数 R0 が 1.0 のときは？

実験として、基本再生産数 R0 が 1.0 のときはどうなるのかを見ていきましょう。

感染者数（緑）のグラフが非常に平らになりほとんど直線になります。
罹患した後に回復した人も全然増えていないように見えますが、実際は回復者が100日目で25人になります。

これから分かることは、

感染がまったく広がらない訳ではない。わずかずつ感染は広がる
しかし、対処可能な感染者数となる

ことを示しています。確率的に感染者は発生するのですから、何処から感染するかは分かりません。しかし、感染して発病したとしても、十分に治療可能な病床が用意できる、ことを意味しています。

github

コードは github moonmile/seir-model: SEIR model simulator に公開しています。

おまけ実効再生産数 R とは？

ここからは私的なメモです。

基本再生産数 R0 と実効再生産数 R の違いを記述しておきます。
専門家会議では、最初の NHK での解説で「基本再生産数 R0」という用語を使っていました。再生産数という言葉から一般的に「基本再生産数を 1.0 以下に抑えれば、なにかグラフが平たくなる」ことが理解できたと思います。

このグラフが出たのが 2/24 の NHK です。

この感染症の流行モデルの元ネタは当時の NHK の解説では何処にも出てきませんが、専門家会議に北大の西浦教授がいることと、西浦教授の2017の論文から、SEIR モデルがベースであることがわかります。

そのあと、色々な批判もあったのち、海外でもこの「平たくする」図がでてくるようになりました。

Flatten the curve | These guidelines are intended to help Flatten the Curve with the COVID19 outbreak, to help limit spread and reduce the load on hospitals and other healthcare.

ただし、これらの「ピークを平たくする」というイメージには、なんら数式的な根拠が示されていません。SEIR モデルを理解しているのか、それとも読者が SEIR モデルを理解できないと思っているのかは定かではないのですが、肝心の数理モデルはでてきません。twitter で #FlattenTheCurve というハッシュタグで拡散されつつありありますが、根拠は示されてません。

結果的に、発症時のピークを下げる＝医療崩壊を防ぐことになるので、効果は同じだと思われるのですが、私が懸念するのは、

ピークを下げた時期の概算が示されていない
ピークを過ぎた後の発症数の裾野が、いつまで続くのか示されていない

のが問題だと思っています。

また、新型コロナウィルスの特徴として無症状な人が多いことから、新しい感染者を0人にすることはほぼ不可能です。天然痘のように全ての人にワクチンをうって撲滅しない限り 0 にはなりません。

同時に「クラスター対策」として言われるように、基本再生産数 R0 は状況によって大きく変わります。全体の平均として 1.0 以下におさえると、マクロの全体として発症者が押さえられるのは確かなのですが、個々の感染しやすい場所、あるいは感染しやすい季節や場所などのミクロの視点で見れば、基本再生産数は大きく変わり分布を持ちます。

このため「実効再生産数 R」として「実効」を付けて、基本再生産数 R0 になんらかのパラメータを掛けます。いわゆる事前条件を加味します。
どうやら疫学の場合の実効再生産数は、二次感染の平均を考慮するようです。このあたり、最近の専門家会議での解説をみると「実効再生産数 R」と区別されていることがわかります。

詳細は2020年3月19日付けの提言を読むとわかります。

「感染拡大地域では自粛検討を」専門家会議が提言【全文】

基本再生産数(R0:すべての者が感受性を有する集団において1人の感染者が生み出した二次感染者数の平均値)
実効再生産数(感染症の流行が進行中の集団のある時刻における、1人の感染者が生み出した二次感染者数の平均値)

基本再生産数は、一般的な都市での感染率 → 都市機能を保ったまま低く保つ必要がある（一般的なイベントの開催、通勤や通学の正常化など）
実効再生産数は、感染してしまった都市での感染率 → 緊急で下げる必要がある（北海道の緊急事態宣言、ロックダウン、都市封鎖）

のように区別します。

カレーの作り方を学ぼう

Tomoaki Masuda — Thu, 26 Mar 2020 01:16:29 +0000

カレーの作り方を学ぼう

var seed = [肉,人参,ジャガイモ,カレールー] ;
var c = new Curry( seed ) ;
c.cut();
c.hot();

カレーができた。

DEV Community: Tomoaki Masuda

隠れSIRモデルを試作して、東京都の感染者を予測する

隠れ感染者が一定量いる

隠れ感染者を推測する

Excel で予測する

データ

参考文献

感染症数理モデルでシミュレーションする（その２）

SEIR モデルを JS で書く

SEIR モデルの数式

SEIR モデルを JS に直す

グラフ化する

基本再生産数 R0 を 1 に近づける

基本再生産数 R0 が 1.0 のときは？

github

おまけ 実効再生産数 R とは？

カレーの作り方を学ぼう

おまけ実効再生産数 R とは？