さんげつ 雑纂

6.1　SEIARモデル

　発症しない人がいる場合を取り扱うため、長田直樹氏がSIRモデルを拡張したSEIARモデルを提唱した（区画に無症状感染者を持つ感染症の数理モデル）。これはグーグルの予測モデルでも使用されている（グーグルの予測モデル(2)）。

　SIRモデルは、可感染者数$S$ 、感染者数$I$、隔離数$R$を使って記述する。SEIARモデルは$I$を細分化する。感染後しばらく感染力を持たない潜伏状態にある（数を$E$、Exposed）と考え、その後それが感染力を持つ発症者（数を従来通り$I$）および感染力を持つが発症しない者（数を$A$、Asymptomatic)に移行すると考える。前者は検査を経て隔離され、後者は感染力喪失という形で隔離される。このモデルを変数の頭文字をとって、SEIARモデルと言う。
　図示すると図6.1のようになる。

図6.1　SEIARモデルの流れ図

可感染者Sは感染するが、しばらく感染力を持たない状態E（潜伏状態）になる。その後分離して、発症者Iと無症状者Aに分離する。前者は検査を経て、後者は感染力を失って隔離状態になる。

　方程式を書き下せば
\begin{align}
\frac{dS(t)}{dt} &= - \beta S(t) (I(t) + \alpha A(t)) \tag{6.1}\\
\frac{dE(t)}{dt} &= \beta S(t)(I(t) + \alpha A(t)) - \epsilon E(t) \tag{6.2}\\
\frac{dI(t)}{dt} &= \epsilon \zeta E(t) - \gamma I(t) \tag{6.3}\\
\frac{dA(t)}{dt} &= \epsilon(1-\zeta) E(t) - \lambda A(t) \tag{6.4}\\
\frac{dR(t)}{dt} &= \gamma I(t) + \lambda A(t) \tag{6.5}
\end{align}
　(6.1)式において、右辺第2項は発症しないが感染力をもつ感染者$A$の感染力が発症者$I$の$\beta $値の$\alpha$倍になっていることを表す。$\alpha$については、CDCが0.75という数字を出しているのでこれを採用する。
　(6.2)式の右辺第3項の$\epsilon$は感染力の無い状態から感染力のある状態への転移を表すが、その時間が4.5日だったので、$\epsilon = 1/4.5\sim 0.2　$　とする。　感染－＞発症－＞ 検査のルートでの隔離に要する平均日数は10日なので、発症－＞検査に要する平均日数は$10- 1/0.2 =5日$であるから、$\gamma = 1/5 = 0.2 （=\gamma_2とおく）$である。
　(6.3)と(6.4)の右辺第1項は、$E$の状態から離脱したとき割合$\zeta$が$I$に、残りの$1-\zeta$が$A$に転移することを示す。ダイヤモンドプリンセス号の例では、$\zeta = 0.44 $　、他の例では$0.3 \le \zeta \le 0.5$であるので、$\zeta＝0.4$とする。
　(6.4)の右辺第2項は、「PCR検査で陽性と判断されてから陰性の結果が出るまでの期間は3〜21日で、最頻値は9日前後」であるので、$\lambda = 1/10 = 0.1$とする。
　残るパラメタ$\beta$は基本再生産数から算出する（第6.2節参照）。

　これらのパラメタは、$\alpha と\gamma$を除いて、長田論文の数字を使っている。同論文では$\alpha$は値が不明として、0.5, 1, 2の場合を検討している。また$\gamma=0.1$を使用しており、感染－＞ 発症－＞検査のルートでの隔離に要する平均日数を15日とみていることが分かる。これは緊急事態宣言が発せられたころの数字である。本ブログではこの日数をより短い10日としている。

6.2　典型例

　まず、典型的な変化を示しておく。$\mathcal R_0= 2.5$の場合を図6.2に示す。

図6.2　SEIARモデルの典型的な変化

感染者数$E, I, A$が一山を作る単純な形をしており、潜伏状態$E$のピークが先行し、その後発者$I$と非発症者$A$が続く。非感染者$S$は単調減少、回復者$R$は単調増加。 $S (0) = 10^7$ （人口1千万人）、 $E(0) = 1, I(0)=A(0)=0 $（初期感染者 1 人）、基本再生産数 $\mathcal R_0$=2.5 の場合。

　SIRモデルと同じように感染者数$E, I, A$が一山を作る単純な形をしており、潜伏状態$E$のピークが先行し、その後発症者$I$と非発症者$A$が続く。これらは最終的には0に近づく。SEIARモデルの感染者の総数$E+I+A$は、SIRモデルの感染者数$I$より高くなっている。

図6.3　SEIARモデルとSIRモデルの比較

SEIARモデルの潜伏状態$E$、感染状態（発症）$I$、感染状態（無症状）$A$の和とSIRモデルの感染状態$I$を比較した。SEIARモデルの感染者の総数$E+I+A$はSIRモデルの感染者数$I$より高い。初期条件は図6.2と同じ。

　SEIARモデルの一般的振る舞いは詳しくは調べられていないと思うが、基本的に図6.2のような単純なものだろうと予想される。一般的には変数の数が3以上ではいわゆるカオスと呼ばれる複雑な挙動を示すことがあるが、SEIARモデルでは図6.1の流れ図から分かるように、下流への動きしかないからである。

6.3　基本再生産数と実効再生産数

　SIRモデルでは基本再生産数$\mathcal R_0$は感染係数$\beta$と

\begin{align}
\mathcal R_0 = \beta S(0)/\gamma \tag{再掲2.7}
\end{align}
で結びついていた。SEIARモデルでは$E,\ I,\ A$の3変数が絡むので関係式は複雑になる。(6.2),(6.3),(6.4)を線型化して$I(t) , E(t), A(t) \ll S(0)$を考慮すると
\begin{align*}
\frac{dE(t)}{dt} &= \beta S(0)(I(t) + \alpha A(t)) - \epsilon E(t) \\
\frac{dI(t)}{dt} &= \epsilon \zeta E(t) - \gamma I(t) \\
\frac{dA(t)}{dt} &= \epsilon(1-\zeta) E(t) - \lambda A(t)
\end{align*}
この解は$\exp(\rho t)$に比例して増大するが、$\rho$はある3次方程式を満たす。この方程式の根で実部が最大になるものに対し（注2）、その実部を$\rho_1$とおく。
　基本再生産数は、指数関数的増大部分を
\begin{align}
\exp(\Gamma (\mathcal R_0 -1)t) \tag{6.6}
\end{align}
と表現する。ここで$\Gamma $は感染してから隔離されるまでの平均時間の逆数として定義され、SEIARモデルでは
\begin{align}
\Gamma = 1/(1/\epsilon + 1/\gamma) = \epsilon\gamma/(\epsilon + \gamma) = 0.1 \tag{6.7}
\end{align}
である。(6.6)から
\begin{align}
\mathcal R_0 = \rho_1 /\Gamma +1 \tag{6.8}
\end{align}
の関係がなりたつ。

（再掲2.7）の関係では基本再生産数$\mathcal R_0$と感染係数$\beta$は比例しているから、接触率の削減はそのまま基本再生産数$\mathcal R_0$の削減に効いた。たとえば接触率を半減すれば基本再生産数も半減する。
　ところが、無症状感染者が存在するSEIARモデルでは、両者の関係は単純な比例関係ではない。この結果、接触の削減の効きが悪いことがおきる。

　図6.4は両者の関係を示した（$\mathcal R_0=2.5$の場合）。この図によれは、例えば接触率をもとの値の5割にしたとき、基本再生産数はSIRモデルでは元の値の5割であるが（灰色直線）、SIEARモデルでは6割程度になる（黒色曲線）。つまり基本再生産数の削減率は4割程度しかない。

　潜伏状態$E$を作り出す段階では、$\beta$の削減は直接（比例関係）で効く。しかし、その後$I$と$A$を作り出す段階には直接的影響を及ぼしていない。このことが基本再生産数削減への効きの悪さの原因になっていると考えられる。簡単に言うと、無発症感染者$A$の存在が効きの悪さの原因である。

図6.4　基本再生産数$\mathcal R_0$と感染係数$\beta $の関係

例えば接触率をもとの値の5割にしたとき、基本再生産数はSIRモデルでは元の値の5割であるが（灰色直線）SIEARモデルでは6割程度になる（黒色曲線）。つまり基本再生産数の削減率は4割程度しかない。

　実効再生産数はベースに$I$、$I+A$、$I+A+E$のいずれを取るかで、3通り定義できる。ここでは$I+A$を取ることにする。
基本再生産数の議論にならって
\begin{align}
\frac{d \log(A+I)}{dt}= \Gamma (\mathcal R_t -1) \tag{6.9}
\end{align}
で$\mathcal R_t$を定義すると、(6.3,6.4)を使って、
\begin{align}
\mathcal R_t = \frac{\epsilon E(t) - \lambda A(t) - \gamma I(t)}{ \Gamma (I(t) + A(t))} + 1 \tag{6.10}
\end{align}
で与えられる。変化の様子を図6.5に示す。

図6.5　実効再生産数の時間変化

基本再生産数$\mathcal R_0=2.5$とした場合。初期条件の影響が消えて、基本再生産数の値に近づく。その後感染拡大で値が低下する。初期条件として$E(0)=I(0)=A(0)=1$とした。

　これをSIRモデルの場合の実効再生産数のグラフ（再掲図2.5、「8割削減とは何だったのか」(2)）と、初期フェイズを除けば非常によく似た変化をしていることが分かる。

再掲図2.5　実効再生産数 $\mathcal R_t $の時間的変化(SIRモデル）

初期フェイズを除けば、SEIARモデルとよく似た変化をしている。基本再生産数$\mathcal R_0=2.5$、初期条件$I(0)=1$。

　なお、$I$または$I+A+E$をベースに取った場合も、初期フェイズを除けば、ここに述べた$I+A$をベースにしたものと大差ない。

6.4　削減問題

　$\beta$の削減率の効きが悪い（基本再生産数$\mathcal R_0$の減少の度合いが小さい）ことは、大きな意味を持つことがある。図6.5はSIRモデルにおいて$\mathcal R_0=2.5$で出発し、途中50日目に$\beta$に6割5分削減と8割削減かけた結果である。6割5分削減では、減衰の大きさは8割削減に比べて小さいものの感染の抑え込みには成功している。

図6.6　感染係数$\beta$の削減を行った結果（SIRモデル）

SIRモデルにおいて$\mathcal R_0=2.5$で出発し、途中50日目に$\beta$に6割5分削減と8割削減かけた結果である。50日以降、6割5分削減をかけた場合（灰色と橙色）は、8割削減（赤色と青色）ほどではないが、減少している。

　同じことをSEIARモデルで行ったのが、図6.7である。この場合、6割5分削減では、感染の抑え込みに失敗している。

図6.7　感染係数$\beta$の削減を行った結果（SEIARモデル）

SIRモデルにおいて$\mathcal R_0=2.5$で出発し、途中50日目に$\beta$に6割5分削減と8割削減かけた結果である。50日以降、6割5分削減をかけた場合（灰色と橙色）は、8割削減（赤色と青色）と異なり、抑え込みに失敗している。

6.5　市中感染者の推定

　実際に観測されるのは日々報告される感染者数である。これから、市中にどの程度の感染者がいるかという推測は重要な意味を持つと考えられる。SIRモデルでは日々報告される感染者数は$\gamma I$である。$\gamma =0.1$であるから、市中の感染者数$I$は日々報告される感染者の$1/\gamma$倍、すなわち10倍である。
　SEIARモデルでは、日々報告される感染者数はSIRモデルと同じく$\gamma I$であるが、市中の（他人を感染させ得る）感染者数は$I + A$である。したがって、倍率（拡大率）は$(I + A)/(\gamma I)$で与えられる。これをプロットしたのが図6.8である。この図には参考のため$I$のグラフも書き入れてある。
　倍率は時期によって異なるのだが、まだ本格的な感染爆発は経験しておらず、このグラフの左側の部分に相当すると考えらえる。倍率は10数倍というところである。

図6.8　拡大率の推移

拡大率＝市中感染者数を日々報告される感染者数で割ったもの。SIRモデルでは常に一定の10である。SEIARモデルでは、時期によって異なるが、感染爆発に至る前では10数倍程度である。

6.6　パラメタの妥当性

　パラメタは、一応根拠はあるものの、かなりばらつきがある中でエイヤと決めたであろうものを使用した。その妥当性は、さまざまな観点からチェックを入れる必要がある。そのひとつが、市中の無症状感染者の占める割合である。

　図6.9は無症状感染者数$A$を感染力を有する感染者数$I+A$で割った　$r_1 = A/(A+I)$　である。時期を判断するために発症感染者数のグラフも書き入れてある。　$r_1$ は0.7～0.8と見積もられる。$E$から$A$と$I$に分離する時、$\zeta = 0.4$なので、$r_1=0.6$である。この値から0.7以上に増えたのは、$A$から$R$への移行に要する日数が10日で、$I$から$R$への5日の倍あり、その分$A$に滞留する数が大きくなるからである。
　$r_1$に関して、CDCのファウチは5割と言い、中国のデータは6割としており（「見えない感染者」と抗体検査）、これらの数字と比較すると、モデルの値はやや過大である。このことから$\zeta$の値を小さくする、$A$から$R$への移行に要する日数を短くする、といった調整が必要に見える。
　なお、非発症者に$E$を入れると、この比は$ r_2= (E+A)/(E+A+I)$となり、これも図6.9に書き入れてある。値は$r_1$より大きく、0.8台になる。$E$の人は、初期にはPCR検査で陽性判定されず、後期になってから陽転すると思う。したがって状感染者の占める割合は$r_1とr_2$の中間に来ることになる。モデルの見積もり値はやや過大に見える。

図6.9　無症状者の占める割合

黒色：無症状感染者数$A$を感染力を有する感染者数$I+A$で割った$r_1 = A/(A+I)$の時間変化、
灰色：潜伏状態の感染者数 $E$ も計算に入れた場合で $r_2 = (E+A)/(E+A+I)$ 。

注　釈

（注1）

　$E(t), I(t), A(t)$が$\exp(\rho t)$に比例した変化をするとして次の特性方程式を得る。
\begin{align*}
\rho E &= \beta S(0)(I + \alpha A) - \epsilon E \\
\rho I &= \epsilon \zeta E - \gamma I \\
\rho A &= \epsilon(1-\zeta) E - \lambda
\end{align*}
　これが自明でない解を持つために
\begin{align*}
\left |
\begin{array}{ccc}
\rho + \epsilon & - \beta S(0) & - \alpha \beta S(0) \\
-\epsilon \zeta & \rho + \gamma & 0 \\
-\epsilon(1-\zeta)& 0 & \rho + \lambda \\
\end{array}
\right |
= 0
\end{align*}
　この行列式を展開して3次方程式
\begin{align}
\rho^3 + a \rho^2 + b\rho + c =0 \tag{6.11}
\end{align}
を得る。ここで
\begin{align*}
a &= \gamma + \epsilon + \lambda \\
b &= \gamma \epsilon + \gamma \lambda + \epsilon\lambda - \alpha \beta \epsilon(1-\zeta) S(0)
-\epsilon \zeta \beta S(0) \\
c &= \gamma \epsilon \lambda - \{\lambda\zeta - \alpha (1-\zeta)\} \beta\epsilon S(0)
\end{align*}

（注2）

　第6.1節のパラメタでは実数になる。虚数部を持つ場合振動しながら増大することになるが、モデルの単純さから見てその可能性は低いように思う。SIRモデルの基本再生産数$\beta S(0)/\gamma$が1より大きい場合、3次方程式(6.11)の根で実部が最大となるものは実根だと予想するが、未確認。