モデルの修正
家庭内など比較的固定した集団内での感染経路に比べて、会食など動的に変化する集団内の感染経路では追跡できる割合が低い可能性を考慮してモデルを修正します。発症した陽性者は感染経路不明者として観測される感染者数に数えられますが、無症状者は未確認キャリアの集団に入ることになります。
今まで同様に時間(t)の関数として、未確認キャリアの数 f(t) と、観測される感染者数 g(t) の二つの関数を考えます。
モデルに登場するパラメータとしては
潜伏期間: Te
感染力のある発症者が自己隔離するまでの時間: Ti
無症状者が感染力を持ったまま活動している時間: kTi
の三つを考えます。全体の感染者数が少ないため、集団免疫により感染が広がらなくなる効果は無視します。
p を感染者の中の無症状者の割合、Co を感染力のある人が固定した集団内で他の人にどれだけ伝染すかを表す定数とします。固定した集団内の感染者は追跡調査で100%つかまるものとします。
さらに C を動的に変化する集団内で他の人にどれだけ伝染すかを表す数とし、d を伝染した人の中で追跡調査で判明する割合とします。
すると微分方程式は以下のようになります。
f'(t) = p(Co+C){f(t-Te)-f(t-Te-kTi)} + dpC{g(t-Te)-g(t-Te-Ti)}
g'(t) = (1-p)(Co+C){f(t-Te)-f(t-Te-kTi)} + {Co+(1-dp)C}{g(t-Te)-g(t-Te-Ti)}
前回同様に Te=3, Ti=2, k=5 とします。
行列を分解すれば解析的に解けそうですが、試行錯誤で近い結果が得られるパラメータを求めるためと、時間経過とともにパラメータを変化させるために、数値計算で計算することにしました。
以下が比較的良くあった結果をグラフにしたものになります。
p=0.31, Co=0.197572, d=0.8 として C は以下のように時間的に変化させています。
～2/13 C=0
2/13～2/20 C=0 から C=0.04 に線形に変化
2/20～2/28 C=0.04
2/28～3/10 C=0.04 から C=0.0936 に線形に変化
3/10～3/25 C=0.0936
3/25～3/30 C=0.0936 から C=0.136 に線形に変化
3/30～ C=0.136
パラメータ数が増えた分、計算結果は当然合いやすくなってます。
本日以降の C は単なる予想ですが、実際にどうなったか 4月に振り返りたいと思います。
予測結果の振り返り(3)
時間の経過とともに実効再生産数も当初に比べてかなり変わって来たので、陽性者数と感染経路不明者の割合の7日間の移動平均のグラフを使って、パラメータの時間的変化を見て行きたいと思います。
まず最初の頃の減少傾向のグラフですが、無症状者の割合0.32、実効再生産数0.7で比較的近い結果が得られます。
この後、減少傾向がかなり鈍化したグラフは、無症状者の割合0.34、実効再生産数0.8で比較的近い結果が得られます。
その後、ほぼフラットになったグラフは、無症状者の割合0.38、実効再生産数0.99(式の都合上1,0は特異点となるため)が比較的近い結果を示します。
時間的変化を見ると、感染経路不明者の割合から推定される無症状者の割合が、時間的に増加する結果となっています。しかし感染者のデータを見ると、この期間の感染者の年齢分布に大きな変化は見られないため、無症状者の割合が変化する理由は見当たりません。最初の単純なモデルでは現実をうまく説明できないという結論にならざるを得ません。最初のモデルでは、実行再生産数が増えると、感染経路不明者の割合は小さくなる傾向があります。今回の結果を説明するには感染経路不明者の割合が大きくなる要因を導入する必要があります。
毎日発表されるデータを見ると感染経路が判明している陽性者の中で、家庭内感染した人数に比べて、会食などで感染した人数が小さい結果となっています。家庭内など比較的固定した集団内での感染経路に比べて、会食など動的に変化する集団内の感染経路では追跡できる割合が低い可能性が考えられます。この場合発症した陽性者は感染経路不明者に数えられますが、無症状者は未確認キャリアの集団に入ることになります。次回はこの影響を考慮してモデルを少し修正した計算結果について説明したいと思います。
予測結果の振り返り(2)
前回に引続き予測結果について振り返ってみたいと思います。
まずは以下のグラフが前回の p=0.33 で計算した予測値と実績値を突き合わせたグラフになります。
週末と祝日の検査数が少ない影響により時々値が外れていますが、それなりに合っているように見えます。ただこの一週間はやや減少が鈍ってきているように見えるのが少し心配です。
予測結果の振り返り
前回の予測結果について振り返ってみたいと思います。
まずは以下のグラフは実績値と突き合わせたグラフになります。
感染経路不明者の割合は予測のように上がらず50%からほぼ変化していません。また陽性者数は予測値に比べやや低くなってます。予測で使用した式は第二波の時の推定値 p=0.42 と感染経路不明者の初期値50%、6日間で75%程度に減少するという仮定で、以下のパラメータを用いました。
p=0.42, C=0.191124, x=-0.02709, y=-0.23806, a=11280.77, b=27141.74, c=456.7896
感染経路不明者の割合の収束値 u(∞) に効くのは p の値になりますが、この値が実際より大きかったためずれた可能性があります。
第二波と比較すると現在の状況には以下の特徴があると思います。
・30代以下の比率が低い(高齢者の比率が高い)ため無症状者の割合が低い可能性がある。
・会食による感染者の数が少なく、家庭内感染や施設内感染による感染者おｎ数が多いため、感染経路が追いやすい。逆に第二波の時は感染経路不明の割合が高めにずれている可能性がある。
6日間で75%程度に減少し、感染経路不明者の割合が50%に収束するように値を選ぶと、以下のパラメータとなります。
p=0.33, C=0.204526, x=-0.04795, y=-0.22144, a=4348.831, b=17509.07, c=32.43973
これをグラフにすると以下の結果となります。週末の検査数が少ない影響により日、月、火の値はずれてますが、前回よりは実測値に近い計算結果が得られるようです。
新型コロナ第三波の東京の感染者数
緊急事態宣言が出て新型コロナウィルス第三波の東京の感染者数もようやく少し減ってきました。前回計算した式を使って今後の陽性者数を予測してみたいと思います。
計算するためには最初にパラメータを決める必要があります。第三波は12月末に感染者数が急激に増え、それと同時に感染経路不明者の割合も上昇しました。その後1月に緊急事態宣言が出た後、陽性者数が徐々に減少するとともに、感染経路不明者の割合も現在は50%程度に下がっています。
一説では12月末は帰省する人等がPCR検査を受信したことにより陽性者数が一気に増えたと言われています。この状態は前回のf(t)に含まれる人が検査により陽性が判明しg(t)に移ったと考えることができます。この段階ではf(t)の人の感染経路は不明なため、感染経路不明者の割合は上昇します。しかしその後はf(t)が減りg(t)が増えることで感染経路不明者の割合は低い状態が続きます。
パラメータは現在の状態を参考に、陽性者が6日間で75%程度に減少するとして、1月29日の陽性者数868人、感染経路不明者の割合50%をもとに、前回のx, y, a, b, c等を求めて計算します。
結果のグラフは以下のようになりました。陽性者数が徐々に減少するのと合わせて、感染経路不明者の割合は少しずつ上昇する結果となっています。
微分方程式にして計算してみた
以前の計算では世代間で感染するモデルを元に、数列として計算していましたが、今回はこれを微分方程式に直してみました。
まず時間(t)の関数として、未確認キャリアの数 f(t) と、観測される感染者数 g(t) の二つの関数を考えます。
モデルに登場するパラメータとしては
潜伏期間: Te
感染力のある発症者が自己隔離するまでの時間: Ti
無症状者が感染力を持ったまま活動している時間: kTi
の三つを考えます。全体の感染者数が少ないため、集団免疫により感染が広がらなくなる効果は無視します。
p を感染者の中の無症状者の割合、C を感染力のある人が他の人にどれだけ移しやすいかを表す定数とすると、微分方程式は以下のようになります。
f'(t) = pC{f(t-Te)-f(t-Te-kTi)}
g'(t) = (1-p)C{f(t-Te)-f(t-Te-kTi)} + C{g(t-Te)-g(t-Te-Ti)}
これを解くと、一般解は a, b, c を定数、Ta, Ts を時定数として、以下の形の式になります。
f(t) = a・exp(t/Ta)
g(t) = b・exp(t/Ta) + c・exp(t/Ts)
これを f'(t) に関する式に代入すると
a/Ta・exp(t/Ta) = apC[exp(-Te/Ta)-exp{-(Te+kTi)/Ta}]・exp(t/Ta)
より、
1/Ta = pC[exp(-Te/Ta)-exp{-(Te+kTi)/Ta}]
exp(-Te/Ta)-exp{-(Te+kTi)/Ta} = 1/pCTa --- (1)
g'(t) に関する式の内 Ts に関する項だけを取り出して考えると
c/Ts・exp(t/Ts) = cC[exp(-Te/Ts)-exp{-(Te+Ti)/Ts}]・exp(t/Ts)
より、
1/Ts = C[exp(-Te/Ts)-exp{-(Te+Ti)/Ts}]
exp(-Te/Ts)-exp{-(Te+Ti)/Ts} = 1/CTs --- (2)
g'(t) に関する式の内 Ta に関する項だけを取り出して考えると
b/Ta・exp(t/Ta) = a(1-p)C[exp(-Te/Ta)-exp{-(Te+kTi)/Ta}・exp(t/Ta) + bC[exp(- Te/Ta)-exp{-(Te+Ti)/Ta}]・exp(t/Ta)
より (1) を使って、
b/Ta = a(1-p)/p/Ta + bC[exp(-Te/Ta)-exp{-(Te+Ti)/Ta}]
bC[exp(-Te/Ta)-exp{-(Te+Ti)/Ta}] = b/Ta - a(1-p)/p/Ta
k > 1 の条件では、Ta < Ts となることを利用すると、t→∞ では Ta に関する項が支配的となるため、t→∞ における感染経路不明者の割合を u(∞) とすると
g'(t→∞) ≒ a(1-p)/p/Ta・exp(t/Ta) + {b/Ta - a(1-p)/p/Ta}・exp(t/Ta)
より
u(∞) = a(1-p)/b/p
となります。
また
g'(0) = (1-p)aC[exp(-Te/Ta)-exp{-(Te+kTi)/Ta}] + bC[exp(-Te/Ta)-exp{-(Te+Ti)/Ta}] + cC{exp(-Te/Ts)-exp(-(Te+Ti)/Ts}]
= a(1-p)/p/Ta + {b-a(1-p)/p}/Ta + c/Ts
= b/Ta + c/Ts
t=0 における感染経路不明者の割合を u(0) とすると
a(1-p)/p/Ta/u(0) = b/Ta + c/Ts
c = {a(1-p)/p/u(0) - b}・Ts/Ta
p, k, Te, Ti が決まっている場合、実行再生産数から C が求まり、これから (1) と (2) を使って、Ta と Ts が計算でき、時間 0 と 無限大における感染経路不明者数の割合から、a, b, c の比率が求まり、方程式が解けます。
数列として計算していた場合と異なるのは、u(∞) が C を含む関数となり、実効再生産数が 1 より大きい場合は感染経路不明者の割合は低くなり、実効再生産数が 1 より小さい場合は感染経路不明者の割合は高くなりまる。
新型コロナ(COVID-19)の感染経路不明の割合と無症状者の実効再生産数
細かい計算の話題が続いていて元の考え方の説明が埋まってきてしまったので、ここで一回まとめておきます。
最初は東京都が発表しているデータで、感染経路不明の割合がだいたい一定の値になっていたのが不思議だったため考え始めたアイディアになります。
感染経路不明の感染者が出続けていることから、検査で捕まっている感染者の他に「未確認キャリア」が全体の中の q という割合で存在すると考えました。「未確認キャリア」は発症していたら捕まるはずなので、おそらく無症状と推測されます。
この「未確認キャリア」が伝染す実効再生産数を Ra、検査で捕まる感染者が伝染す実効再生産数を Rs とします。Ra は Rs に比べて小さな値かもしれないし、大きな値かもしれませんが、一般的に異なる値とします。Ra と Rs の効果が合わさって、全体の実効再生産数 (R) となります。
感染者の中の無症状者の割合を p とした時、一つの世代が次の世代に伝染すのを図で表すと以下のようになります。
# 図中のpと1-pが入れ替わっていたので訂正しました。
充分に時間が経過してそれぞれの比率が一定の値に収束した定常状態を考えると、感染経路不明者の割合と全体の実効再生産数(R)から、Ra と Rs がそれぞれ R に比例する形で求まります。
実際に計算してみると Ra は Rs の 5倍程度の値となります。発症する場合、発症前二日間ほど無症状だけど人に伝染す期間があると言われています。発症した後は、通常は自己隔離するため、それ以降人に伝染す可能性は下がります。
新型コロナウィルスの場合、唾液による飛沫感染と言われています。唾液によるPCR 検査では、発症後九日間以上経つと、唾液中のウィルスの量が減少して感度が悪くなると言われています。もし無症状の場合も移す仕組みが発症者と同じで
症状が出ていないだけとすると、2+9日の間、人に伝染する期間があることになります。
この比率はだいたい 5 程度であり、無症状者の実効再生産数が高い理由は、無症状者が社会に滞留しやすいからと考えられます。
