制御系CAD

制御を意識したモノつくりを目指して

感染制御

投稿日時: 2020年9月5日投稿者: cacsd

1. SIRモデルに基づく感染制御
1.1 SIRモデルの基礎式

感染ダイナミックスを表すSIRモデルでは、当該地域の人口 $N$ を

● $S(t)$ ：感受性(susceptible)状態の人口（未感染者数）
● $I(t)$ ：感染性(infectious)状態の人口（感染者数）
● $R(t)$ ：回復・隔離または免疫状態(recovered/removed/immune)状態の人口（回復者数）

に分け（ $N=S+I+R$ ）、これらが次の連立微分方程式に従うとします。

$\displaystyle{(1)\quad \left\{\begin{array}{l} \frac{dS(t)}{dt}=-\beta S(t)I(t) \\ \frac{dI(t)}{dt}=-\gamma I(t)+\beta S(t)I(t)\\ \frac{dR(t)}{dt}=\gamma I(t) \end{array}\right. }$

ここで、 $\beta$ は感染率、 $\gamma$ は回復率と呼ばれます。

いま未感染者数 $S$ 、感染者数 $I$ 、回復者数 $R$ の当該地域の人口 $N$ に対する比率を考えると次式が成り立ちます。

$\displaystyle{(2)\quad \underbrace{\frac{I(t)}{N}}_{i(t)}+\underbrace{\frac{R(t)}{N}}_{r(t)}+\underbrace{\frac{S(t)}{N}}_{s(t)}=1 }$

(1)と(2)から

$\displaystyle{(3)\quad \left\{\begin{array}{l} \underbrace{\frac{d\frac{I(t)}{N}}{dt}}_{\frac{di(t)}{dt}}=-\gamma \underbrace{\frac{I(t)}{N}}_{i(t)}+\beta N(\underbrace{\frac{N}{N}}_{1}-\underbrace{\frac{I(t)}{N}}_{i(t)}-\underbrace{\frac{R(t)}{N}}_{r(t)})\underbrace{\frac{I}{N}}_{i(t)}\\ \underbrace{\frac{d\frac{R(t)}{N}}{dt}}_{\frac{dr(t)}{dt}}=\gamma \underbrace{\frac{I(t)}{N}}_{i(t)} \end{array}\right. }$

すなわち、次の正規化（無次元化）されたSIRモデルを得ます。

$\displaystyle{(4)\quad \left\{\begin{array}{l} \frac{di(t)}{dt}=\gamma(-1+\underbrace{\frac{\beta N}{\gamma}}_{\bar{R}_0} \underbrace{(1-i(t)-r(t))}_{s(t)}) i(t)\\ \frac{dr(t)}{dt}=\gamma i(t) \end{array}\right. }$

ここで、自粛率 $q(t)$ を導入し、 $\beta$ を $(1-q(t))\beta$ に置き換えて、次式を得ます。

$\displaystyle{(5)\quad \left\{\begin{array}{l} \frac{di(t)}{dt}=(-1+(1-q(t))\underbrace{\frac{\beta N}{\gamma}}_{\bar{R}_0} \underbrace{(1-i(t)-r(t))}_{s(t)}) \gamma i(t)\\ \frac{dr(t)}{dt}=\gamma i(t) \end{array}\right. }$

図1にSIRモデル(5)に基づく感染シミュレーションの例を示します。ただし、自粛率 $q$ の時系列は８割削減と５割削減を適当に設定しています（ $1-q$ は削減率と呼ばれます）。

図1 正規化SIRモデルに基づく感染シミュレータの例(MATLAB/Simulink)
（ $\gamma=1$ ; $\bar{R}_0=2.5$ ; $q(t)=0,0.8,0.5$ ; $i(0)=0.003$ ）

1.2 SIRモデル（全国版）に基づく感染予測

図2に全国の新規感染者数の移動平均値を示します。これに基づく基本再生産数R0の値を図3に示します。これから7/25時点でのR0の値は $1.5$ 程度（ $\gamma=0.16$ ）と推測されます。

図2 全国の新規感染者数の実績値

図3 全国の基本再生産数の推定値

さらに感染予測結果を図4に示します。左図は7/1を初期時刻にとって、自粛は行われていないものとして、7/25現在の感染者数を予測したものです。右図は、7/25を初期時刻にとって、自粛率 $q$ を５割削減と２割削減を適当に設定して、向こう300日間の感染者数を予測したものです。新規感染者数の実績値は2000人を超えることはなかったので、もう少し厳しい自粛が行われていたと言えます。

図4 SIRモデルに基づく全国の感染予測の例

1.3 正規化SIRモデルに基づく感染制御

適当な累積感染者数と収束期間を実現するように自粛率 $q(t)$ の時系列を制御理論を用いて導出することを考えます。(5)の行列表示は次式となります。

$\displaystyle{(6)\quad \left[\begin{array}{c} \dot{r}(t)\\ \dot{i}(t) \end{array}\right] =\gamma\left[\begin{array}{cc} 0 & 1\\ 0 & -1+(1-q(t))\bar{R}_0s(t) \end{array}\right] \left[\begin{array}{c} r(t)\\ i(t) \end{array}\right] }$

さらに自粛率 $q(t)$ を１次遅れのダイナミックスを通して間接的に操作できるように政策変数（仮想的な操作変数） $u_q(t)$ を導入します。

$\displaystyle{(7)\quad \left\{\begin{array}{l} \left[\begin{array}{c} \dot{r}(t)\\ \dot{i}(t) \end{array}\right] =\left[\begin{array}{cc} 0 & \gamma\\ 0 & -\gamma+\gamma\bar{R}_0s(t) \end{array}\right] \left[\begin{array}{c} r(t)\\ i(t) \end{array}\right] +\left[\begin{array}{c} 0\\ -\gamma\bar{R}_0s(t)i(t) \end{array}\right]q(t)\\ \dot{q}(t)=-\frac{1}{T_q}q(t)+\frac{1}{T_q}u_q(t) \end{array}\right. }$

すなわち、次式を得ます。

$\displaystyle{(8)\quad \left[\begin{array}{c} \dot{r}(t)\\ \dot{i}(t)\\ \dot{q}(t) \end{array}\right] =\left[\begin{array}{ccc} 0 & \gamma &0\\ 0 & -\gamma+\gamma\bar{R}_0s(t)&-\gamma\bar{R}_0s(t)i(t)\\ 0 & 0 &-\frac{1}{T_q} \end{array}\right] \left[\begin{array}{c} r(t)\\ i(t)\\ q(t) \end{array}\right] +\left[\begin{array}{c} 0\\ 0\\ \frac{1}{T_q} \end{array}\right]u_q(t) }$

これは２つの変動パラメータ $\alpha(t)=\bar{R}_0s(t)$ と $\beta(t)=\bar{R}_0s(t)i(t)$ をもつ、次のようなLPVモデル（線形パラメータ変動モデル）とみなすことができます。

$\displaystyle{(9)\quad \begin{array}{l} \underbrace{ \left[\begin{array}{c} \dot{r}(t)\\ \dot{i}(t)\\ \dot{q}(t) \end{array}\right] }_{\dot{x}(t)} =\left(\underbrace{\left[\begin{array}{ccc} 0 & \gamma &0\\ 0 & -\gamma &0\\ 0 & 0 &-\frac{1}{T_q} \end{array}\right] }_{A_0} + \underbrace{\bar{R}_0s(t)}_{\alpha(t)} \underbrace{\left[\begin{array}{ccc} 0 & 0 & 0\\ 0 & \gamma & 0\\ 0 & 0 & 0 \end{array}\right] }_{A_1} \right.\\ \left.+ \underbrace{\bar{R}_0s(t)i(t)}_{\beta(t)} \underbrace{\left[\begin{array}{ccc} 0 & 0 & 0\\ 0 & 0 & -\gamma\\ 0 & 0 & 0 \end{array}\right] }_{A_2} \right) \underbrace{\left[\begin{array}{c} r(t)\\ i(t)\\ q(t) \end{array}\right] }_{x(t)} +\underbrace{\left[\begin{array}{c} 0\\ 0\\ \frac{1}{T_q} \end{array}\right] }_{B} \underbrace{u_q(t)}_{u(t)}\\ (0<\alpha(t),\beta(t)<\bar{R}_0) \end{array} }$

いま、２つの変動パラメータを適当な値に $\alpha^*=\bar{R}_0s^*$ と $\beta^*=\bar{R}_0s^*i^*$ に固定すると、次のLTIモデル（線形時不変モデル）を得ます。

$\displaystyle{(10)\quad \left[\begin{array}{c} \dot{r}(t)\\ \dot{i}(t)\\ \dot{q}(t) \end{array}\right] =\left[\begin{array}{ccc} 0 & \gamma &0\\ 0 & -\gamma+\gamma\alpha^* &　-\gamma\beta^*\\ 0 & 0 &-\frac{1}{T_q} \end{array}\right] \left[\begin{array}{c} r(t)\\ i(t)\\ q(t) \end{array}\right] +\left[\begin{array}{c} 0\\ 0\\ \frac{1}{T_q} \end{array}\right]u_q(t) }$

そこで、２次形式評価関数

$\displaystyle{(11)\quad J=\int_0^\infty(r^2(t)+i^2(t)+q^2(t)+\rho^2u_q^2(t))dt }$

を最小化する状態フィードバック（LQ制御則）

$\displaystyle{(12)\quad u_q(f)=-f_1r(t)-f_2i(t)-f_3q(t) }$

の効果を確かめてみます。その結果を図5に示します。左図は自粛なしの場合、右図は時刻 $t=3$ 以降、LQ制御則を適用しています。ここで、 $s^*=0.85$ 、 $i^*=0.15$ 、 $T_q=0.2$ としています。時刻 $t=3$ ではすでに人口の15%（0.15 $N$ ）の感染者がいますので、強力な自粛が必要となり、 $q$ の値は１を超える大きな値となっています。そのためシミュレータ上では飽和要素を入れてこれを制限しています。結果的にいきなり10割削減相当、いわゆるロックダウンの自粛となっています。しかし最終的な累積感染者数は9割から３割程度に抑えられています。

図5 正規化SIRモデルに基づく感染制御（LQ制御）の例
（ $\gamma=1$ ; $\bar{R}_0=2.5$ ; $i(0)=0.003$ ）

さて、 $\alpha_1\le\alpha(t)\le\alpha_2$ 、 $\beta_1\le\beta(t)\le\beta_2$ に対する次の内分式に注目します。

$\displaystyle{(13)\quad \begin{array}{l} \left[\begin{array}{c} \alpha(t)\\ \beta(t) \end{array}\right]= \underbrace{\frac{\alpha_2-\alpha(t)}{\alpha_2-\alpha_1}\frac{\beta_2-\beta(t)}{\beta_2-\beta_1}}_{p_{11}(\alpha(t),\beta(t))}\left[\begin{array}{c} \alpha_1\\ \beta_1 \end{array}\right]+ \underbrace{\frac{\alpha_2-\alpha(t)}{\alpha_2-\alpha_1}\frac{\beta(t)-\beta_1}{\beta_2-\beta_1}}_{p_{12}(\alpha(t),\beta(t))}\left[\begin{array}{c} \alpha_1\\ \beta_2 \end{array}\right]\\+ \underbrace{\frac{\alpha(t)-\alpha_1}{\alpha_2-\alpha_1}\frac{\beta_2-\beta(t)}{\beta_2-\beta_1}}_{p_{21}(\alpha(t),\beta(t))}\left[\begin{array}{c} \alpha_2\\ \beta_1 \end{array}\right]+ \underbrace{\frac{\alpha(t)-\alpha_1}{\alpha_2-\alpha_1}\frac{\beta(t)-\beta_1}{\beta_2-\beta_1}}_{p_{22}(\alpha(t),\beta(t))}\left[\begin{array}{c} \alpha_2\\ \beta_2 \end{array}\right] \\ \end{array} }$

$\displaystyle{(14)\quad p_{11}(\alpha(t),\beta(t))+p_{12}(\alpha(t),\beta(t))+p_{21}(\alpha(t),\beta(t))+p_{22}(\alpha(t),\beta(t))=1 }$

状態方程式(9)は、端点モデル

$\displaystyle{(15)\quad \begin{array}{l} \dot{x}(t)=\underbrace{A(\alpha_1,\beta_1)}_{A_{11}}x(t)+Bu(t)\\ \dot{x}(t)=\underbrace{A(\alpha_1,\beta_2)}_{A_{12}}x(t)+Bu(t)\\ \dot{x}(t)=\underbrace{A(\alpha_2,\beta_1)}_{A_{21}}x(t)+Bu(t)\\ \dot{x}(t)=\underbrace{A(\alpha_2,\beta_2)}_{A_{22}}x(t)+Bu(t) \end{array} }$

を、 $p_{11},p_{12},p_{12},p_{22}$ によって重み付けして、次式のように表されます。

$\displaystyle{(16)\quad \dot{x}(t)=\underbrace{(p_{11}(t)A_{11}+p_{12}(t)A_{12}+p_{21}(t)A_{21}+p_{22}(t)A_{22})}_{A(\alpha(t),\beta(t))}x(t)+Bu(t) }$

このとき４つの端点モデルを安定化する状態フィードバック

$\displaystyle{(17)\quad \begin{array}{l} u(t)=-F_{11}x(t)\\ u(t)=-F_{12}x(t)\\ u(t)=-F_{21}x(t)\\ u(t)=-F_{22}x(t) \end{array} }$

を同時安定化して求め（共通のリャプノフ関数を求め）、LPV制御則を次式で決めます。

$\displaystyle{(18)\quad u(t)=-\underbrace{(p_{11}(t)F_{11}+p_{12}(t)F_{12}+p_{21}(t)F_{21}+p_{22}(t)F_{22})}_{F(\alpha(t),\beta(t))}x(t) }$

この制御則の一例の効果を図6に示します。自粛率は１を超えることはありませんが、累積感染者数は９割から６割程度にしか抑制されていません。

図6 正規化SIRモデルに基づく感染制御（LPV制御）の例
（ $\gamma=1$ ; $\bar{R}_0=2.5$ ; $i(0)=0.003$ ）

2. SIQRモデルに基づく感染制御
2.1 SIQRモデルの基礎式

SIRモデルは、次の３つの変数で表されていました。

● $S(t)$ ：感受性状態の人口（未感染者数）
● $I(t)$ ：感染性状態の人口（感染者数）
● $R(t)$ ：回復・隔離または免疫状態状態の人口（回復者数）

いま、感染者数を検査の有無によって、次のように区別します（ $I=H+Q$ ）。

● $H(t)$ ：未検査感染性人口（市中感染者数）
● $Q(t)$ ：検査後隔離状態にある感染性人口（隔離感染者数）

このとき、次式で表されるSIQRモデルを考えます。

$(101)\quad\left\{\begin{array}{l} \frac{dS(t)}{dt}=-\beta S(t)H(t) \\ \frac{dH(t)}{dt}=-\gamma H(t)-pH(t)+(1-q)\beta S(t)H(t)\\ \frac{dQ(t)}{dt}=-\gamma Q(t)+pH(t)+q\beta S(t)H(t)\\ \frac{dR(t)}{dt}=\gamma (H(t)+Q(t)) \end{array}\right.$

ここで、SIQRモデルの第２式と第３式は、SIRモデルの第２式を、 $I=H+Q$ を用いて分けたものとなっています。その際、自粛率 $q$ の他に隔離率 $p$ も導入されていることに注意してください。 $pH(t)$ が日々発表される感染者数となります。

これら未感染者数、市中感染者数、隔離感染者数、回復者数の当該地域の人口 $N$ に対する比率を考えます。

$(102)\quad \underbrace{\frac{H(t)}{N}}_{h(t)}+\underbrace{\frac{Q(t)}{N}}_{q_u(t)}+\underbrace{\frac{R(t)}{N}}_{r(t)}+\underbrace{\frac{S(t)}{N}}_{s(t)}=1$

このとき

$(103)\quad\left\{\begin{array}{l} \underbrace{\frac{d\frac{H(t)}{N}}{dt}}_{\frac{dh(t)}{dt}}=-\gamma \underbrace{\frac{H(t)}{N}}_{h(t)}-p \underbrace{\frac{H(t)}{N}}_{h(t)}+(1-q)\beta N(\underbrace{\frac{N}{N}}_{1}-\underbrace{\frac{H(t)}{N}}_{h(t)}-\underbrace{\frac{Q(t)}{N}}_{q_u(t)}-\underbrace{\frac{R(t)}{N}}_{r(t)})\underbrace{\frac{H(t)}{N}}_{h(t)}\\ \underbrace{\frac{d\frac{Q(t)}{N}}{dt}}_{\frac{dq_u(t)}{dt}}=-\gamma \underbrace{\frac{Q(t)}{N}}_{q_u(t)}+p \underbrace{\frac{H(t)}{N}}_{h(t)}+q\beta N(\underbrace{\frac{N}{N}}_{1}-\underbrace{\frac{H(t)}{N}}_{h(t)}-\underbrace{\frac{Q(t)}{N}}_{q_u(t)}-\underbrace{\frac{R(t)}{N}}_{r(t)})\underbrace{\frac{H(t)}{N}}_{h(t)}\\ \underbrace{\frac{d\frac{R(t)}{N}}{dt}}_{\frac{dr(t)}{dt}}=\gamma \underbrace{\frac{H(t)}{N}}_{h(t)}+\gamma \underbrace{\frac{Q(t)}{N}}_{q_u(t)} \end{array}\right.$

すなわち、次の正規化（無次元化）されたSIQRモデルを得ます。

$(104)\quad\left\{\begin{array}{l} \frac{dh(t)}{dt}=-(\gamma+p)h(t)+(1-q)\underbrace{\frac{\beta N}{\gamma}}_{\bar{R}_0} \underbrace{(1-h(t)-q_u(t)-r(t))}_{s(t)} \gamma h(t)\\ \frac{dq_u(t)}{dt}=-{\gamma}q_u(t)+{p}h(t)+q\underbrace{\frac{\beta N}{\gamma}}_{\bar{R}_0} \underbrace{(1-h(t)-q_u(t)-r(t))}_{s(t)} \gamma h(t)\\ \frac{dr(t)}{dt}=\gamma(h(t)+q_u(t)) \end{array}\right.$

図7にSIQRモデル(104)に基づく感染シミュレーションの例を示します。ここで自粛率 $q$ の時系列は図1の半分にしています。これを隔離率 $p$ を導入することにより、図1と同程度の感染者数に抑え込むことができています。

図7 正規化SIQRモデルに基づく感染シミュレータの例(MATLAB/Simulink)
（ $\gamma=1; \bar{R}_0=2.5; q=0,0.8,0.5; p=0.1,0.5; i(0)=0.003$ ）

2.2 SIQRモデル（全国版）に基づく感染予測

図8 SIQRモデルに基づく全国の感染予測シミュレーション

2.3 正規化SIQRモデルに基づく感染制御

適当な累積感染者数と収束期間を実現するように自粛率 $q(t)$ と隔離率 $q(t)$ の時系列を制御理論を用いて導出してみます。(104)の行列表示は次式となります。

$(105)\quad\left[\begin{array}{c} \dot{r}(t)\\ \dot{h}(t)\\ \dot{q}_u(t) \end{array}\right] =\gamma\left[\begin{array}{ccc} 0 & 1 & 1\\ 0 & -1-\frac{p}{\gamma}+(1-q)\bar{R}_0s(t) & 0\\ 0 & \frac{p}{\gamma}+q\bar{R}_0s(t) & -1\\ \end{array}\right] \left[\begin{array}{c} r(t)\\ h(t)\\ q_u(t) \end{array}\right]$

自粛率 $q(t)$ と隔離率 $q(t)$ を１次遅れのダイナミックスを通して間接的に操作できるようにそれぞれの政策変数（仮想的な操作変数） $u_q(t)$ と $u_p(t)$ を導入します。

$(106)\quad\left\{\begin{array}{l} \left[\begin{array}{c} \dot{r}(t)\\ \dot{h}(t)\\ \dot{q}_u(t) \end{array}\right] =\left[\begin{array}{ccc} 0 & \gamma & \gamma\\ 0 & -\gamma+\gamma\bar{R}_0s(t) & 0\\ 0 & 0 & -\gamma\\ \end{array}\right] \left[\begin{array}{c} r(t)\\ h(t)\\ q_u(t) \end{array}\right]\\+ \left[\begin{array}{cc} 0 & 0\\ -\gamma\bar{R}_0s(t)h(t) & -1\\ \gamma\bar{R}_0s(t)h(t) & 1 \end{array}\right] \left[\begin{array}{c} q(t) \\ p(t) \end{array}\right]\\ \dot{q}(t)=-\frac{1}{T_q}q(t)+\frac{1}{T_q}u_q(t)\\ \dot{p}(t)=-\frac{1}{T_p}p(t)+\frac{1}{T_p}u_p(t) \end{array}\right.$

すなわち、次式を得ます。

$(107)\quad\begin{array}{l} \left[\begin{array}{c} \dot{r}(t)\\ \dot{h}(t)\\ \dot{q}_u(t)\\ \dot{q}(t)\\ \dot{p}(t) \end{array}\right] =\left[\begin{array}{ccccc} 0 & \gamma & \gamma&0 & 0\\ 0 & -\gamma+\gamma\bar{R}_0s(t) & 0&-\gamma\bar{R}_0s(t)h(t) & -1\\ 0 & 0 & -\gamma&\gamma\bar{R}_0s(t)h(t) & 1\\ 0&0&0&-\frac{1}{T_q}&0\\ 0&0&0&0&-\frac{1}{T_p} \end{array}\right] \left[\begin{array}{c} r(t)\\ h(t)\\ q_u(t)\\ q(t) \\ p(t) \end{array}\right]\\ +\left[\begin{array}{cc} 0&0\\ 0&0\\ 0&0\\ \frac{1}{T_q}&0\\ 0&\frac{1}{T_p} \end{array}\right] \left[\begin{array}{c} u_q(t) \\ u_p(t) \end{array}\right] \end{array}$

これは２つの変動パラメータ $\alpha(t)=\bar{R}_0s(t)$ と $\beta(t)=\bar{R}_0s(t)h(t)$ をもつ次のようなLPVモデル（線形パラメータ変動モデル）とみなすことができます。

PCR検査

投稿日時: 2020年9月5日投稿者: cacsd

[1] 現在、PCR検査拡大について賛否が分かれています。たとえば
●左翼メディアと煽り系専門家たちの罪
●無症状にPCR検査して隔離する世田谷モデルはどこがアホなのか説明する
●日本がコロナ第3波を避けられない決定的弱点
などを読むと、どれも説得力があって迷ってしまいます。ただ、一般的に言われている賛否の理由を整理してみると次のようになるかと思います。

●検査隔離拡大の反対理由：
(1) 偽陰性者が市中感染を広める（症状のある者・濃厚接触者に限定すべき）
(2) 偽陽性者数は偽陰性者数を上まわり、検査の無駄が多い（厚労省の特異度は99%、世界標準はほぼ100%）

●検査隔離拡大の賛成理由：
(3) 無症状者に繰り返し検査を行ない、検査陽性者を隔離することが、市中感染抑制につながる
(4) 無症状で検査陰性であることが経済を回すための国際ルールとなりつつある

個人的には(3)が重要だと思うのですが、今一つ理屈がよく分からないというのが、筆者ばかりでなく大方の率直な印象ではないでしょうか？その理解の一助となることを目指して、ここでは感染ダイナミックスに基づいて、一体どれくらいの人数を検査・隔離すればよいのか、具体的には市中感染者の何割ぐらいを検査・隔離すればよいのかを調べてみたいと思います。

[2] いま、当該地域の人口を $N$ 、感染率を $\beta$ 、回復率を $\gamma$ 、削減率を $q$ 、隔離率を $p$ とすると、感染が収束に向かうには、次式が成り立つ必要があります。

$\displaystyle{(5)\quad p>\gamma((1-q)R_0-1) }$

ここで、 $R_0=\frac{\beta N}{\gamma}$ は基本再生産数と呼ばれています。自粛をしていない場合は削減率 $q=0$ なので次式となります。

$\displaystyle{(6)\quad p>\gamma(R_0-1) }$

また市中感染者数を $H$ で表すと、これに隔離率 $p$ をかけた $p\times H$ が日々報告される感染者数となります。したがって、日々報告される感染者数を $p$ で割れば、検査すべき市中感染者数 $H$ を見積もることができます。

いま(5)の右辺を

$\displaystyle{(7)\quad p^*(R_0,q,\gamma)=\gamma((1-q)R_0-1) }$

とおいて、 $R_0$ を1.5, 2.0, 2.5, 3.0の４通り、 $q$ を0, 0.2, 0.5, 0.8の４通り、 $\gamma$ を0.04, 0.08, 0.12, 0.16の４通り、すべての組合せ $4^3=64$ 通りについて、 $p^*(R_0,q,\gamma)$ の値を計算したものを、表1の左側に示します（負値を与える組合せはすでに収束に向かっていることを意味します）。また対応する $1/p^*(R_0,q,\gamma)$ を表1の左側に示します。これは検査拡大率を意味します。その棒グラフを描いた図1から次が分かります。

●基本再生数 $R_0$ が小さいほど検査拡大率は大きい
●削減率 $q$ が大きいほど検査拡大率は大きい

したがって、検査拡大は、基本再生数 $R_0$ が小さい（感染力が小さい）ときほど、削減率 $q$ が大きい（自粛レベルが高い）ときほど必要と言えます。すなわちよく言われるように感染が一旦収まったときに徹底した検査を行うべきという主張を裏付けるものと思われます。

表1 $p^*(R_0,q,\gamma)$ と $1/p^*(R_0,q,\gamma)$

図1 表1の $1/p^*(R_0,q,\gamma)$ の棒グラフ

注感染ダイナミックスの著名なものはSIRモデルと言われます。これに自粛率 $1-q$ （ $q$ は削減率）ばかりでなく、隔離率 $p$ を導入し、SIQRモデルとして提案されたのは小田垣先生です。詳しくは感染ダイナミックスをご覧ください。そこで示している初期成長率(31)式から上の不等式(5)が得られます。あくまで感染初期に適用できる式です。

[3] さて、図1に全国の新規感染者数の実績値の時系列を示します。これを見ると、５月下旬から６月までは第１波が収束していた時期です。ところが７月から急速に感染者数が増加していきます。７月は正に自粛なしの期間と言ってよいかと思います。図2に対応する基本生産数を図3に示します。７月の１か月間では、 $R_0=1.5$ 、 $\gamma=0.12$ と読み取ることができ、(7)は次式となります。

$\displaystyle{(7)\quad p^*(1.5,0,0.12)=\gamma(R_0-1)=0.12\times 0.5=0.06 }$

したがって、全国レベルで日々1000人前後の感染者数が報告されていたので、この約17倍(1/0.06=17)の検査が毎日必要だったと言えそうです。さらに３か月間週１回のペースで繰り返し12回検査するとすると、約20万人規模の検査が必要だったと言えそうです。

以上は、少し大雑把で強引な議論ですが、経済を回すために自粛せずに、感染収束に向かわせるためにはかなりの検査規模が必要になることが分かります。このことが、反対理由(1),(2)から生じるデメリット（検査の無駄）を打ち消すことになるのかは明確ではありませんが、ご参考になればと思います。

図3 全国の新規感染者数の実績値

図4 全国の基本再生産数の推定値

感染ダイナミックス（2020/7/28)

投稿日時: 2020年7月28日投稿者: cacsd

感染ダイナミックスはSIRモデルによって表すことができます。図1は、これをMATLAB/Simulinkで表したものです。

図1 正規化SIRモデルに基づく感染シミュレータの例(MATLAB/Simulink)
（ $\gamma=1$ ; $\bar{R}_0=2.5$ ; $q(t)=0,0.8,0.5$ ; $i(0)=0.003$ ）

図2に全国の新規感染者数の移動平均値を示します。これに基づく基本再生産数 $\bar{R}_0$ の値を図3に示します。これから7/25時点での $\bar{R}_0$ の値は $1.5$ 程度（ $\gamma=0.16$ ）と推測されます。さらに予測結果を図4に示します。左図は7/1を初期時刻にとって、自粛は行われていないものとして、7/25現在の感染者数を予測したものです。右図は、7/25を初期時刻にとって、自粛率 $q$ を５割削減と２割削減を適当に設定して、向こう300日間の感染者数を予測したものです。

図2 全国の新規感染者数の実績値

図3 全国の基本再生産数の推定値

図4 SIRモデルに基づく全国の感染予測の例

これから早急に自粛要請をかけないと感染者数は爆発的に増加していくことが分かります。ただ自粛率という数値を行動変容にどう結び付けるかは不明です。我々の感染拡大への想像力と感染対策の知恵が問われています。

感染ダイナミックス（2020/6/26）

投稿日時: 2020年6月26日投稿者: cacsd

ここでは次の文献にもとづいて「感染ダイナミックス」について述べます。
[1] 稲葉寿編著：「感染症の数理モデル」,培風館（2008）
[2] 牧野淳一郎：「3.11以後の科学リテラシーno.89」,科学, vol.90, no.5, 岩波書店(2020)
[3] 小田垣孝：「新型コロナウイルスの蔓延に関する一考察」
なお、感染ダイナミックスの数理モデルは、次の文献が先駆けとなっています。
Kermack, W. O. and McKendrick, A. G.（1927）. Contributions to the mathematical theory of epidemics – I, Proceedings of the Royal Society Series A, 115, 700–721

1 SIRモデル
1.1 基本再生産数の概念
感染ダイナミックスの数理モデルのうち最も基本的なものはSIRモデルとして知られています。ここでは文献[1]に従って、その概要を述べます。

いま当該地域の人口 $N$ は次のような構成であるとします。

$(1)\quad S(t)+I(t)+R(t)=N$

ただし、時刻 $t$ において

● $S(t)$ ：感受性(susceptible)状態の人口
● $I(t)$ ：感染性(infectious)状態の人口
● $R(t)$ ：回復・隔離または免疫状態(recovered/removed/immune)状態の人口

以下では $S$ 、 $I$ 、 $R$ をそれぞれ未感染者数、感染者数、回復者数と呼びます。いくつかの仮定のもと、感染ダイナミックスの支配方程式として、次のSIRモデルを考えます。

$(2)\quad\left\{\begin{array}{l} \frac{dS(t)}{dt}=-\beta S(t)I(t) \\ \frac{dI(t)}{dt}=-\gamma I(t)+\beta S(t)I(t)\\ \frac{dR(t)}{dt}=\gamma I(t) \end{array}\right.$

ここで、 $\beta$ は感染率、 $\gamma$ は回復率と呼ばれています。

いま時刻 $t$ における感染状況を３項組 $(S(t),I(t),R(t))$ で表します。感染前の自明な平衡状態は

$(3)\quad (S,I,R)=(N,0,0)$

ですが、これが初期時刻 $t=0$ において、感染者がごく少数 $I_0$ だけ現れ

$(4)\quad (S_0,I_0,R_0)=(N-I_0,I_0,0)$

となったとします。このときの関心事は、感染流行が始まるか、それとも感染収束するかどうかです。これは(2)の第２式において、感染初期では $S(t)\simeq S_0$ として

$(5)\quad \frac{dI(t)}{dt}=\underbrace{(\beta S_0-\gamma)}_{\lambda_0}I(t), \quad I(0)=I_0$

の解

$(6)\quad I(t)=\exp(\lambda_0 t)I_0$

の振舞いを調べることになります。明らかに、初期成長率 $\lambda_0$ の正負によって、感染流行か感染収束かが決まります。これを

$(7)\quad \lambda_0=\beta S_0-\gamma=\gamma(\frac{\beta S_0}{\gamma}-1)$

と変形して、ここで現れた定数

$(8)\quad R_0=\frac{\beta S_0}{\gamma}$

を基本再生産数と名付け

● $R_0>1$ のとき感染流行
● $R_0<1$ のとき感染収束

となると判断します。この方法は閾値原理と呼ばれています。

さて、基本再生産数は次のように表されます。

$(9)\quad R_0=\beta S_0\underbrace{\frac{1}{\gamma}}_{T}=1+\underbrace{(\beta S_0-\gamma)}_{\lambda_0}\underbrace{\frac{1}{\gamma}}_{T}$

ここで、 $T$ は感染性時間（感染力をもつ時間）を表すと解釈できます。実際、感染力がなくなったと仮定すると( $\beta=0$ )、(2)の第２式は

$(10)\quad\frac{dI(t)}{dt}=-\gamma I(t), \quad I(0)=I_0$

となり、この解は

$(11)\quad I(t)=\exp(-\gamma t)I_0\rightarrow0\quad(t\rightarrow\infty)$

と表され、 $T=1/\gamma$ の5～10倍の時間が経てば、感染者数 $I_0$ は十分零になるからです。

そこで(9)において、 $\beta S_0$ は１人が単位時間に２次感染させる人数、これに感染性時間 $T$ をかけるので、基本再生産数は感染者１人が２次感染させる人数を表すと言えます。

その算出には初期成長率 $\lambda_0$ が必要ですが、次式を用います。

$(12)\quad \frac{I(t+\Delta t)}{I(t)}=\frac{\exp(\lambda_0 (t+\Delta t))I_0}{\exp(\lambda_0 t)I_0}=\exp(\lambda_0 \Delta t) \Rightarrow \lambda_0=\ln(\frac{I(t+\Delta t)}{I(t)})/\Delta t$

1.2 基本再生産数の推定例
専門家による基本再生産数の推定法は、感染学上の知見に基づく統計推測手法によるものと思われますが、詳細は明らかではありません。そこで、次のような基本再生産数の推定手法を考えてみます。

基本再生産数 $R_0$ の推定手法１：
$1^\circ$ 当該地域の感染者数データの移動平均をとる
$2^\circ$ 次式を用いて時刻 $t$ における（初期）成長率 $\lambda_0$ を計算する

$(13)\quad\lambda_0(t)=\ln(\frac{I(t)}{I(t-\Delta t)})/\Delta t$

$3^\circ$ 次式を用いて時刻 $t$ における基本再生産数 $R_0$ を計算する

$(14)\quad R_0(t)=1+\lambda_0(t) T$

図1に全国、東京都、福岡県、NYの感染数を7日間移動平均したデータを示します。図2に日ごとに基本再生産数 $R_0$ を推定した例を示します。ただし、 $2^\circ$ で $\Delta t=7$ としています。 $3^\circ$ で回復率は $\gamma=0.04,0.08,0.12,0.16$ の４通りを試し、 $R_0<0$ となる場合は0としています。

図1 感染者数の実データ（移動平均後）
国内データソース：NHK 特設サイト新型コロナウイルス
NYデータソ－ス：New York State Department of Health COVID-19 Tracker

図2 基本再生産数 $R_0$ の推定例

図2の妥当性を検討するために、「新型コロナウイルス感染症対策の状況分析・提言」（2020年5月1日）に記載されている次の公表値に注目します。

●全国の実効再生産数：　 $R_0(3/25)=2$ , 　 $R_0(4/10)=0.7$
●東京都の実効再生産数： $R_0(3/14)=2.6$ , $R_0(4/10)=0.5$

ここで基本再生産数の代わりに実効再⽣産数という術語が用いられています。ここでは両者の理論上の違いまでには踏み込まず、実効再⽣産数は基本再生産数 $R_0$ に実質的に等しいとして話を進めます。上記の値は上記資料に掲載されている次の図版から読み取ることができます。

図3 専門家会議資料における実効再生産数の推定

図3の左側は全国について、右側は東京都について、感染者数と実効再生産数が示されています。下側に示されている棒グラフは上側の感染者数の移動平均をとったものと思われ、これに基づいて実効再生産数が計算されています。ただし、移動平均後のデータは２週間ほど前倒しされています（ピーク値と $R_0$ が１を切る時点が少しずれています）。実際、全国（東京都）の感染者数のピークは4/11（4/9）ですが、移動平均後は3/28（3/27）にピークが現れています。図2において全国の基本再生産数が0.7まで下がるのは早くて4/20当たり、東京都の基本再生産数が0.5まで下がるのは早くて4/23当たりですから、ほぼ２週間遅れとなっています。また、回復率 $\gamma$ の値に依らず、感染者数がピークを取る時点で $R_0$ が１を切り、収束に向かっている傾向が捉えられています。

問題はどの回復率 $\gamma$ が妥当かですが、文献でよく見るのは $\gamma=0.04$ です。しかしこれでは直近で $R_0$ が負となりました（ $\gamma=0.12$ , $0.16$ では常に正でした）。回復率の逆数は感染性時間 $T=1/\gamma$ ですから、 $\gamma=0.04$ , $0.08$ , $0.12$ , $0.16$ に対応して、 $T=25$ , $12.5$ , $8.3$ , $6.25$ です。指数関数的に減衰する場合は、初期値の10%まで減衰する期間は $2.3T$ 、2%まで減衰する期間は $3.9T$ ですから、たとえば、 $\gamma=0.16$ , $T=6.25$ でもそうおかしくはないかと思います。

2 感染シミュレーション
2.1 無次元化SIRモデル
前章では、SIRモデルのうち第２式のみを用いて議論しました。ここでは３本の微分法方程式を連立させて、 $S$ 、 $I$ 、 $R$ の時系列のシミュレーションを行います。

$S$ 、 $I$ 、 $R$ の間には(1)の拘束があるので、実質的には２本の微分方程式を連立させて解きます。(2)の第２式と(1)からSを消去して、第３式を合わせると次式を得ます（第１式はこれらから導出できます）。

$(15)\quad\left\{\begin{array}{l} \frac{dI(t)}{dt}=-\gamma I(t)+\beta \underbrace{(N-I(t)-R(t))}_{S(t)}I(t)\\ \frac{dR(t)}{dt}=\gamma I(t) \end{array}\right.$

ここで３つのパラメータ $N$ 、 $\beta$ 、 $\gamma$ がありますが、無次元化という操作を行って、１つのパラメータにまとめます。これは文献[2]を参考にしています。

まず、人口 $N$ への依存を除くために、各変数の割合（人口比）、すなわち、感染者数の割合、回復者数の割合、未感染者数の割合を、それぞれ

$(16)\quad i(t)=I(t)/N, r(t)=R(t)/N, s(t)=S(t)/N$

のように小文字で表します。(1)より次式が成り立ちます。

$(17)\quad \underbrace{\frac{I(t)}{N}}_{i(t)}+\underbrace{\frac{R(t)}{N}}_{r(t)}+\underbrace{\frac{S(t)}{N}}_{s(t)}=1$

これを用いて、(15)から次式を得ます。

$(18)\quad\left\{\begin{array}{l} \underbrace{\frac{d\frac{I(t)}{N}}{dt}}_{\frac{di(t)}{dt}}=-\gamma \underbrace{\frac{I(t)}{N}}_{i(t)}+\beta N(\underbrace{\frac{N}{N}}_{1}-\underbrace{\frac{I(t)}{N}}_{i(t)}-\underbrace{\frac{R(t)}{N}}_{r(t)})\underbrace{\frac{I}{N}}_{i(t)}\\ \underbrace{\frac{d\frac{R(t)}{N}}{dt}}_{\frac{dr(t)}{dt}}=\gamma \underbrace{\frac{I(t)}{N}}_{i(t)} \end{array}\right.$

次に、感染性時間 $T=1/\gamma$ が単位時間となるように、次のように変形します。

$(19)\quad\left\{\begin{array}{l} \underbrace{\frac{di(t)}{d\frac{t}{T}}}_{\frac{di(t)}{d\tau}}=-i(t)+\underbrace{\frac{\beta N}{\gamma}}_{\bar{R}_0} \underbrace{(1-i(t)-r(t))}_{s(t)} i(t)\\ \underbrace{\frac{dr(t)}{d\frac{t}{T}}}_{\frac{dr(t)}{d\tau}}=i(t) \end{array}\right.$

すなわち、 $\tau=t/T$ とおいて、次式を得ます。

$(20)\quad\left\{\begin{array}{l} {\frac{di(\tau)}{d\tau}}=-i(\tau)+\bar{R}_0 s(\tau)i(\tau)=\underbrace{(-1+\bar{R}_0 s(\tau))}_{\lambda(\tau)}i(\tau)\\ {\frac{dr(\tau)}{d\tau}}=i(\tau) \end{array}\right.$

ここで、次の定数が定義されていることに注意してください。

$(21)\quad \bar{R}_0=\frac{\beta N}{\gamma}>\frac{\beta S_0}{\gamma}$

すなわち、 $\bar{R}_0$ は基本再生産数 $R_0$ の上界を与えています。ただし、感染初期の実データでは $S_0\simeq N$ が成り立ち、 $\bar{R}_0\simeq R_0$ です。以下では、 $R_0$ は $\bar{R}_0$ の意味で使用しています。

注1：(20)における $i(\tau)$ , $s(\tau)$ , $r(\tau)$ , $\lambda(t)$ はそれぞれ $i(\tau T)$ , $s(\tau T)$ , $r(\tau T)$ , $\lambda(\tau T)$ として、離散時間の差分方程式として表すべきかもしれませんが、以下では上式を用います。

2-2 感染シミュレーション
SIRモデル(20)の解の振舞いは、数値計算を行って容易にチェックできます。図3は $i(0)=0.001$ のとき、基本再生産数 $R_0$ の値によって感染ダイナミックスがどのように変わるかを示しています。青線が未感染者数 $S$ 、赤線が感染者数 $I$ 、緑線が回復者数（＝累積感染者数） $R$ の割合です。これから、基本再生産数 $R_0$ の値によって、最終的にどれくらいの割合で感染者が出るか、どれくらいの収束期間がかかるかを読み取ることができます。ただし、どの図も $R_0$ が一定の場合であることに注意してください。これらの図から、もし $R_0$ の値を順次下げていくことができれば、最終的な感染者数の割合は小さくなることが分かります。

図4 基本再生産数 $R_0$ に依存する感染ダイナミックス

注2：図4において、時間軸の単位は感染性時間 $T=1/\gamma$ です。したがって、収束期間を見積もるためには、回復率 $\gamma$ の推定が必要になります。

図5は、いま基本再生産数は $R_0=2.5$ であるとし、ある時点からこれを $(1-q)R_0$ に低減したとき、感染者数がどのように変化するかを示しています( $i(0)=0.001$ の場合です)。 $q=0.8$ のときが8割低減、 $q=0.2$ のときが2割低減です。

図5 基本再生産数 $R_0=2.5$ を $(1-q)R_0$ に低減したときの効果

注3：「８割削減の妥当性」についてですが、これはツイート_R020405が基礎となっているようです。ここでは基本再生産数 $R_0=2.5$ の場合であることが表明されていますので図5が参考になります。これは「基本再生産数を８割減することが必要で、２割減では効果がない」という主張をほぼ裏付けていると言えます。

2-2 感染初期における予測
実データによる感染予測ではSIRモデル(20)を解くわけですが、その際 $R_0$ を推定する必要があります。そのために次の問題設定を考えます。

基本再生産数 $R_0$ の推定手法２：
いま人口を $N$ 、感染者数の時系列を $I(t_1),\cdots,I(t_n)$ とします。またパラメータ $R_0$ を含むSIRモデル(20)の解の $t_1,\cdots,t_n$ における値を $i(t_1),\cdots,i(t_n)$ とします。このとき評価関数

$(22)\quad J(R_0)=(I(t_1)/N-i(t_1))^2+\cdots+(I(t_n)/N-i(t_n))^2$

を最小化するようなパラメータ $R_0$ の値 $R_0^*$ を、数値最適化の手法を用いて求めます。

この手法を用いて、全国の1/24～5/5の新規感染者数の一部のデータ（1/24～4/7、赤色〇）をもとに、基本再生産数を推定したところ $R_0^*$ =1.08となりました。図6は、感染者数の実績値（〇）と、SIRモデルによる計算値（青色実線）を重ねてプロットしたものです。感染初期の感染者数の指数関数的増加傾向がよく捉えられています。
図6 日本における感染者数の推移

首相は4/7日夜の記者会見で『このペースで感染拡大が続けば（感染者が）２週間後には１万人、１カ月後には８万人を越える』と指摘しました。図7は図6の感染者数の累積値を示しています。起点となる日を4/7（75日目）としますと、２週間後（88日目）より早く１万人を超えています（実際に１万人を超えたのは4/18（86日目）でした）。また８万人に到達するのは112日目あたりで、１か月後(105日目)より少し遅くなっています。
図7 日本における累積感染者数の推移

注4：上の問題を解くとき、少数の感染者が出た日を第１日として、 $t_1=1$ 、以後 $t_2=2,t_3=3\cdots$ とします。一方、SIRモデル(20)の解も $t_1=1,t_2=2,\cdots$ において求めます。この方法で感染初期の指数関数的増大の様子をよく捉えることができます。一方、4/11放映のNHKスペシャルより、基本再生産数 $R_0$ は大体1.5程度であることが表明されています。上で求めた $R_0^*$ の値1.08は公表値1.5より低めに出ています。そこで、日ごとの実データに回帰をとる場合 $T=1$ と考えられるので、(9)より $初期成長率\lambda_0$ を求めると

$(23)\quad \underbrace{1.08}_{R_0^*}=1+\lambda_0\times\underbrace{1}_{T}\quad\Rightarrow\quad \lambda_0=0.08}$

となります。これをたとえば感染性時間 $T=1/\gamma=1/0.16=6.25$ を用いて

$(24)\quad R_0=1+\underbrace{0.08}_{\lambda_0}\times\underbrace{6.25}_{T}=1.5}$

のように修正すると、公表値と一致します。

3. SIQRモデル
5/5に小田垣先生が新型コロナウイルスの蔓延に関する一考察を発表されました。現在指摘されている「PCR検査不足」の主張の妥当性を議論する際のベースを与えることのできるきわめて興味深い論文です。以下では、本稿の枠組みの中でフォローアップさせていただきたいと思います。

当該論文では、感染性人口 $I(t)$ を

● $H(t)$ ：未検査感染性人口（市中感染者数）　⇒　筆者の方で記号をIからHに変更
● $Q(t)$ ：検査後隔離状態にある感染性人口（隔離感染者数）

に分けて、次のSIQRモデルが提案されています。

$(25)\quad\left\{\begin{array}{l} \frac{dS(t)}{dt}=-\beta S(t)H(t) \\ \frac{dH(t)}{dt}=-\gamma H(t)-pH(t)+(1-q)\beta S(t)H(t)\\ \frac{dQ(t)}{dt}=-\gamma Q(t)+pH(t)+q\beta S(t)H(t)\\ \frac{dR(t)}{dt}=\gamma (H(t)+Q(t)) \end{array}\right.$

これをSIRモデル

$(2)\quad\left\{\begin{array}{l} \frac{dS(t)}{dt}=-\beta S(t)I(t) \\ \frac{dI(t)}{dt}=-\gamma I(t)+\beta S(t)I(t)\\ \frac{dR(t)}{dt}=\gamma I(t) \end{array}\right.$

と比較します。それぞれの第 $i$ 式を、SIRi、SIQRiと参照します。

●SIQR4はSIR3に $I=H+Q$ を代入しただけです。
●SIQR1はSIR1の $I=H+Q$ のうち隔離感染者数 $Q$ の影響はないとしたものです。
●SIQR2とSIQR3は、SIR2の $I=H+Q$ を分けたものに相当します。SIR2の $\beta S(t)I(t)$ については

$(26)\quad \beta S(t)I(t)=(1-q)\beta S(t)I(t)+q\beta S(t)I(t)$

と分解し、 $I=H+Q$ のうち隔離感染者数 $Q$ の影響はないとしています。ここでパラメータ $q$ は感染率 $\beta$ に対する削減率とみなせます。またSIQR2とSIQR3には、 $pH(t)$ が差し引きしてあります。この $pH(t)$ は検査後陽性判定を受け隔離される人口で、日々発表される感染者数とされています。ここでパラメータ $p$ は市中感染者に対する隔離率とみなせます。

このとき、SIRモデルの場合の(5)と同様に

$(27)\quad \frac{dH(t)}{dt}=\underbrace{((1-q)\beta S_0-\gamma-p)}_{\lambda_0(q,p)}H(t), \quad H(0)=H_0$

が成り立ち、初期成長率 $\lambda_0(q,p)$ は

$(28)\quad \lambda_0(q,p)=(\gamma+p)((1-q)\underbrace{\frac{\beta S_0}{\gamma+p}}_{\tilde{R}_0}-1)$

と表すことができます。これをSIRモデルの場合の

$(7)\quad \lambda_0=\gamma(\underbrace{\frac{\beta S_0}{\gamma}}_{R_0}-1)$

と比較してみますと、 $\gamma$ が $\gamma+p$ に代わっており、基本再生産数が小さくなっています。したがって、SIQRモデルによれば、隔離率 $p$ を上げれば基本再生産数が小さくなる可能性を理解できます。

(27)の解は

$(29)\quad H(t)=\exp(\lambda_0(q,p)t)H_0\rightarrow0\quad(t\rightarrow\infty)$

と表され、初期値 $H_0$ の10%まで下がる時間は( $\lambda_0(q,p)<0$ に注意して)

$(30)\quad \exp(\lambda_0(q,p)t)=0.1\quad\Rightarrow\quad t=\frac{\ln(0.1)}{|\lambda_0(q,p)|}=\frac{2.3}{|\lambda_0(q,p)|}$

で求められます。

当該論文の重要な貢献は、初期成長率を

$(31)\quad \lambda_0=\beta S_0-\gamma\quad\Rightarrow\quad \lambda_0(q,p)=(1-q)\beta S_0-(\gamma+p)$

のように、隔離率 $p$ と削減率 $q$ で変化させることができることを明らかにしていることです。そして、当該論文では、全国の感染者数データについて、次の結論を示しています。

Case 1) 現状： $\lambda_0(0.6,p)=-0.075$ 　⇒　収束期間31日
Case 2) ８割自粛： $\lambda_0(0.8,p)=-0.101$ 　⇒　収束期間23日
Case 3) ロックダウン： $\lambda_0(1,p)=-0.126$ 　⇒　収束期間18日
Case 4) 検査４倍増： $\lambda_0(0,4p)=-0.288$ 　⇒　収束期間8日
Case 5) ５割自粛＋検査２倍増： $\lambda_0(0.5,2p)=-0.159$ 　⇒　収束期間14日

これらは次を(31)と(30)に代入して確認できます。

$(32)\quad \beta S_0=0.126$ 、 $\gamma=0.03$ 、 $p=0.096$

たとえば、Case 1)については

$(33)\quad \underbrace{0.4}_{1-q}\times \underbrace{0.126}_{\beta S_0}-\underbrace{0.096}_{p}-\underbrace{0.03}_{\gamma}=\underbrace{-0.075}_{\lambda_0}\quad\Rightarrow\quad \underbrace{2.3/0.075}_{2.3/|\lambda_0|}=31$

当該論文では、基本再生産数を

感染頂上期間： $\tilde{R}_0=\frac{\beta S_0}{\gamma+p}=\frac{0.126}{0.03+0.096}=1$
感染減少期間： $\tilde{R}_0=\frac{(1-q)\beta S_0}{\gamma+p}=\frac{0.051}{0.03+0.096}=0.405$

のように計算しています。これから現状の削減率は $q=1-0.051/0.126=0.595$ であることが確かめられます。また、回復率としては $\gamma=0.03$ 、感染性期間が $T=1/\gamma=33$ が採用されています。さらに、感染減少期間に入る前の感染頂上期間では $\tilde{R}_0=1$ となるように、隔離率 $p=0.096$ 、 $1/p=10.4$ が定められています。このように、感染者数データから隔離率 $p$ を見積もることができ、 $1/p$ 倍して市中感染者数を推定できるという意味で大変興味深いと思います。

注5：本稿の計算では、基本再生産数は

感染増大期間： $R_0=\frac{\beta S_0}{\gamma}=\frac{0.128}{0.04}=3.191$
感染減少期間： $R_0=\frac{\beta S_0}{\gamma}=\frac{0.022}{0.04}=0.545$

となります。そこで当該論文と同様の計算をしてみますと

感染頂上期間： $\tilde{R}_0=\frac{\beta S_0}{\gamma+p}=\frac{0.128}{0.04+0.088}=1$
感染減少期間： $\tilde{R}_0=\frac{(1-q)\beta S_0}{\gamma+p}=\frac{0.022}{0.04+0.088}=0.545$

となります。これから削減率 $q=1-0.022/0.128=0.829$ 、隔離率 $p=0.088$ 、 $1/p=11.4$ が分かります。これらを用いると次の計算結果を得ます。

Case 2) ６割自粛： $\lambda_0(0.6,p)=-0.077$ 　⇒　収束期間30日
Case 1) 現状： $\lambda_0(0.829,p)=-0.106$ 　⇒　収束期間22日
Case 3) ロックダウン： $\lambda_0(1,p)=-0.128$ 　⇒　収束期間18日
Case 4) 検査４倍増： $\lambda_0(0,4p)=-0.263$ 　⇒　収束期間9日
Case 5) ５割自粛＋検査２倍増： $\lambda_0(0.5,2p)=-0.151$ 　⇒　収束期間15日

本稿のデータ分析では、すでに８割強の自粛が行われている結果になっていますが、PCR検査を増やすことの効果を確認できたと言えます。ただし、偽陽性と偽陰性の可能性を勘案して、陽性の出る可能性の高い場合に検査を行う必要があることは述べるまでもありません。

感染ダイナミックス（2020/5/17)

投稿日時: 2020年5月17日投稿者: cacsd

このサイトはコロナ禍が始まって最初にまとめました。ただし、図1と図2は随時更新しています。また感染制御についてもまとめようとしています。

1 はじめに
新型コロナウィルスに関する検討結果を公開しています。「なぜ外出自粛が必要かを説得力をもって説明することは難しい」という問題意識から検討を始めました（ツイート_R020330、ツイート_R020402）。また「８割削減」に示される数値の根拠に感心が集まっています。筆者の専門は感染学ではなく、システム制御工学ですが、一定の説明はできると思っています。

ポイントは「感染ダイナミックス」への理解です。ここでは次の文献を参照します。
[1] 稲葉寿編著：「感染症の数理モデル」,培風館（2008）
[2] 牧野淳一郎：「3.11以後の科学リテラシーno.89」,科学, vol.90, no.5, 岩波書店(2020)
なお、感染ダイナミックスの数理モデルは、次の文献が先駆けとなっています。
Kermack, W. O. and McKendrick, A. G.（1927）. Contributions to the mathematical theory of epidemics – I, Proceedings of the Royal Society Series A, 115, 700–721

2 SIRモデル
2.1 基本再生産数の概念
感染ダイナミックスの数理モデルのうち最も基本的なものはSIRモデルとして知られています。ここでは文献[1]に従って、その概要を述べます。

いま当該地域の人口 $N$ は次のような構成であるとします。

$(1)\quad S(t)+I(t)+R(t)=N$

ただし、時刻 $t$ において

● $S(t)$ ：感受性(susceptible)状態の人口
● $I(t)$ ：感染性(infectious)状態の人口
● $R(t)$ ：回復・隔離または免疫状態(recovered/removed/immune)状態の人口

$(2)\quad\left\{\begin{array}{l} \frac{dS(t)}{dt}=-\beta S(t)I(t) \\ \frac{dI(t)}{dt}=-\gamma I(t)+\beta S(t)I(t)\\ \frac{dR(t)}{dt}=\gamma I(t) \end{array}\right.$

ここで、 $\beta$ は感染率、 $\gamma$ は回復率と呼ばれています。

いま時刻 $t$ における感染状況を３項組 $(S(t),I(t),R(t))$ で表します。感染前の自明な平衡状態は

$(3)\quad (S,I,R)=(N,0,0)$

ですが、これが初期時刻 $t=0$ において、感染者がごく少数 $I_0$ だけ現れ

$(4)\quad (S_0,I_0,R_0)=(N-I_0,I_0,0)$

$(5)\quad \frac{dI(t)}{dt}=\underbrace{(\beta S_0-\gamma)}_{\lambda_0}I(t), \quad I(0)=I_0$

の解

$(6)\quad I(t)=\exp(\lambda_0 t)I_0$

の振舞いを調べることになります。明らかに、初期成長率 $\lambda_0$ の正負によって、感染流行か感染収束かが決まります。これを

$(7)\quad \lambda_0=\beta S_0-\gamma=\gamma(\frac{\beta S_0}{\gamma}-1)$

と変形して、ここで現れた定数

$(8)\quad R_0=\frac{\beta S_0}{\gamma}$

を基本再生産数と名付け

● $R_0>1$ のとき感染流行
● $R_0<1$ のとき感染収束

となると判断します。この方法は閾値原理と呼ばれています。

さて、基本再生産数は次のように表されます。

$(9)\quad R_0=\beta S_0\underbrace{\frac{1}{\gamma}}_{T}=1+\underbrace{(\beta S_0-\gamma)}_{\lambda_0}\underbrace{\frac{1}{\gamma}}_{T}$

ここで、 $T$ は感染性時間（感染力をもつ時間）を表すと解釈できます。実際、感染力がなくなったと仮定すると( $\beta=0$ )、(2)の第２式は

$(10)\quad\frac{dI(t)}{dt}=-\gamma I(t), \quad I(0)=I_0$

となり、この解は

$(11)\quad I(t)=\exp(-\gamma t)I_0\rightarrow0\quad(t\rightarrow\infty)$

と表され、 $T=1/\gamma$ の5～10倍の時間が経てば、感染者数 $I_0$ は十分零になるからです。

その算出には初期成長率 $\lambda_0$ が必要ですが、次式を用います。

$(12)\quad \frac{I(t+\Delta t)}{I(t)}=\frac{\exp(\lambda_0 (t+\Delta t))I_0}{\exp(\lambda_0 t)I_0}=\exp(\lambda_0 \Delta t) \Rightarrow \lambda_0=\ln(\frac{I(t+\Delta t)}{I(t)})/\Delta t$

2.2 基本再生産数の推定例
専門家による基本再生産数の推定法は、感染学上の知見に基づく統計推測手法によるものと思われますが、詳細は明らかではありません。そこで、次のような基本再生産数の推定手法を考えてみます。

$(13)\quad\lambda_0(t)=\ln(\frac{I(t)}{I(t-\Delta t)})/\Delta t$

$3^\circ$ 次式を用いて時刻 $t$ における基本再生産数 $R_0$ を計算する

$(14)\quad R_0(t)=1+\lambda_0(t) T$

図2 基本再生産数 $R_0$ の推定例

図2の妥当性を検討するために、「新型コロナウイルス感染症対策の状況分析・提言」（2020年5月1日）に記載されている次の公表値に注目します。

●全国の実効再生産数：　 $R_0(3/25)=2$ , 　 $R_0(4/10)=0.7$
●東京都の実効再生産数： $R_0(3/14)=2.6$ , $R_0(4/10)=0.5$

図3 専門家会議資料における実効再生産数の推定

3 感染シミュレーション
3.1 無次元化SIRモデル
前章では、SIRモデルのうち第２式のみを用いて議論しました。ここでは３本の微分法方程式を連立させて、 $S$ 、 $I$ 、 $R$ の時系列のシミュレーションを行います。

$(15)\quad\left\{\begin{array}{l} \frac{dI(t)}{dt}=-\gamma I(t)+\beta \underbrace{(N-I(t)-R(t))}_{S(t)}I(t)\\ \frac{dR(t)}{dt}=\gamma I(t) \end{array}\right.$

まず、人口 $N$ への依存を除くために、各変数の割合（人口比）、すなわち、感染者数の割合、回復者数の割合、未感染者数の割合を、それぞれ

$(16)\quad i(t)=I(t)/N, r(t)=R(t)/N, s(t)=S(t)/N$

のように小文字で表します。(1)より次式が成り立ちます。

$(17)\quad \underbrace{\frac{I(t)}{N}}_{i(t)}+\underbrace{\frac{R(t)}{N}}_{r(t)}+\underbrace{\frac{S(t)}{N}}_{s(t)}=1$

これを用いて、(15)から次式を得ます。

次に、感染性時間 $T=1/\gamma$ が単位時間となるように、次のように変形します。

すなわち、 $\tau=t/T$ とおいて、次式を得ます。

ここで、次の定数が定義されていることに注意してください。

$(21)\quad \bar{R}_0=\frac{\beta N}{\gamma}>\frac{\beta S_0}{\gamma}$

3-2 感染シミュレーション
SIRモデル(20)の解の振舞いは、数値計算を行って容易にチェックできます。図3は $i(0)=0.001$ のとき、基本再生産数 $R_0$ の値によって感染ダイナミックスがどのように変わるかを示しています。青線が未感染者数 $S$ 、赤線が感染者数 $I$ 、緑線が回復者数（＝累積感染者数） $R$ の割合です。これから、基本再生産数 $R_0$ の値によって、最終的にどれくらいの割合で感染者が出るか、どれくらいの収束期間がかかるかを読み取ることができます。ただし、どの図も $R_0$ が一定の場合であることに注意してください。これらの図から、もし $R_0$ の値を順次下げていくことができれば、最終的な感染者数の割合は小さくなることが分かります。

図4 基本再生産数 $R_0$ に依存する感染ダイナミックス

図5 基本再生産数 $R_0=2.5$ を $(1-q)R_0$ に低減したときの効果

$(22)\quad J(R_0)=(I(t_1)/N-i(t_1))^2+\cdots+(I(t_n)/N-i(t_n))^2$

を最小化するようなパラメータ $R_0$ の値 $R_0^*$ を、数値最適化の手法を用いて求めます。

$(23)\quad \underbrace{1.08}_{R_0^*}=1+\lambda_0\times\underbrace{1}_{T}\quad\Rightarrow\quad \lambda_0=0.08}$

となります。これをたとえば感染性時間 $T=1/\gamma=1/0.16=6.25$ を用いて

$(24)\quad R_0=1+\underbrace{0.08}_{\lambda_0}\times\underbrace{6.25}_{T}=1.5}$

のように修正すると、公表値と一致します。

5 おわりに
「なぜ外出自粛は必要か」の問いには、感染爆発を防ぐためには基本再生産数を小さく抑え込むことが必要で、そのための間接的手段として外出自粛やロックダウンが考えられるからと答えられます。感染ダイナミックスの簡易モデル（SIRモデル）は発散から収束に切替る振舞いをうまくモデル化しています。ただ、基本再生産数というたった１つのパラメータで、すべての状況を反映させようとしているところには自ずと限界があります。それでも感染ダイナミクスの本質を捉えていることは間違いないと思います。

「８割削減の妥当性」についてですが、注3で述べたように、基本再生産数 $R_0=2.5$ の場合、図5から「基本再生産数を８割削減することが必要で、２割削減では効果がない」という主張は理解できます。ところが、4/11放映のNHKスペシャルより、基本再生産数 $R_0$ は大体1.5程度であることが表明されています。この場合、図5は違ったものになります。また、5/1の専門家会議では図3のように実効再生産数は4/10の時点で1を切っていることが表明されています。いつの時点で、どのくらいの基本再生産数であるかは、感染対策を議論する上で極めて重要です。もし4/7の時点で基本再生産数が１を切っているのなら、非常事態宣言は医療崩壊を防ぐ意味合いが強かったのかもしれません。

感染流行を制御する決め手は基本再生産数 $R_0$ をどう調整していくかにかかっています。これをSIRモデル(18)の非線形最適制御問題または線形パラメータ変動制御問題として定式化し、その解として望ましい $R_0$ の時系列を求めることが考えられます。もちろん、これが得られたとしても、どのような行動変容に結び付ければよいかは不明です。ただ、厳しい外出自粛をいつまでも続けるわけにはいきません。現在収束過程にあると思われる感染者数データから $R_0$ の時系列を求め、比較検討しておくことは、第２波への対応や出口戦略の検討に役立つと思っております。年内には検討結果をご報告したいと思います。

補遺
5/5に小田垣先生が新型コロナウイルスの蔓延に関する一考察を発表されました。現在指摘されている「PCR検査不足」の主張の妥当性を議論する際のベースを与えることのできるきわめて興味深い論文です。以下では、本稿の枠組みの中でフォローアップさせていただきたいと思います。

当該論文では、感染性人口 $I(t)$ を

● $H(t)$ ：未検査感染性人口（市中感染者数）　⇒　筆者の方で記号をIからHに変更
● $Q(t)$ ：検査後隔離状態にある感染性人口（隔離感染者数）

に分けて、次のSIQRモデルが提案されています。

これをSIRモデル

$(2)\quad\left\{\begin{array}{l} \frac{dS(t)}{dt}=-\beta S(t)I(t) \\ \frac{dI(t)}{dt}=-\gamma I(t)+\beta S(t)I(t)\\ \frac{dR(t)}{dt}=\gamma I(t) \end{array}\right.$

と比較します。それぞれの第 $i$ 式を、SIRi、SIQRiと参照します。

$(26)\quad \beta S(t)I(t)=(1-q)\beta S(t)I(t)+q\beta S(t)I(t)$

このとき、SIRモデルの場合の(5)と同様に

$(27)\quad \frac{dH(t)}{dt}=\underbrace{((1-q)\beta S_0-\gamma-p)}_{\lambda_0(q,p)}H(t), \quad H(0)=H_0$

が成り立ち、初期成長率 $\lambda_0(q,p)$ は

$(28)\quad \lambda_0(q,p)=(\gamma+p)((1-q)\underbrace{\frac{\beta S_0}{\gamma+p}}_{\tilde{R}_0}-1)$

と表すことができます。これをSIRモデルの場合の

$(7)\quad \lambda_0=\gamma(\underbrace{\frac{\beta S_0}{\gamma}}_{R_0}-1)$

(27)の解は

$(29)\quad H(t)=\exp(\lambda_0(q,p)t)H_0\rightarrow0\quad(t\rightarrow\infty)$

と表され、初期値 $H_0$ の10%まで下がる時間は( $\lambda_0(q,p)<0$ に注意して)

$(30)\quad \exp(\lambda_0(q,p)t)=0.1\quad\Rightarrow\quad t=\frac{\ln(0.1)}{|\lambda_0(q,p)|}=\frac{2.3}{|\lambda_0(q,p)|}$

で求められます。

当該論文の重要な貢献は、初期成長率を

$(31)\quad \lambda_0=\beta S_0-\gamma\quad\Rightarrow\quad \lambda_0(q,p)=(1-q)\beta S_0-(\gamma+p)$

これらは次を(31)と(30)に代入して確認できます。

$(32)\quad \beta S_0=0.126$ 、 $\gamma=0.03$ 、 $p=0.096$

たとえば、Case 1)については

当該論文では、基本再生産数を

感染頂上期間： $\tilde{R}_0=\frac{\beta S_0}{\gamma+p}=\frac{0.126}{0.03+0.096}=1$
感染減少期間： $\tilde{R}_0=\frac{(1-q)\beta S_0}{\gamma+p}=\frac{0.051}{0.03+0.096}=0.405$

注5：本稿の計算では、基本再生産数は

感染増大期間： $R_0=\frac{\beta S_0}{\gamma}=\frac{0.128}{0.04}=3.191$
感染減少期間： $R_0=\frac{\beta S_0}{\gamma}=\frac{0.022}{0.04}=0.545$

となります。そこで当該論文と同様の計算をしてみますと

感染頂上期間： $\tilde{R}_0=\frac{\beta S_0}{\gamma+p}=\frac{0.128}{0.04+0.088}=1$
感染減少期間： $\tilde{R}_0=\frac{(1-q)\beta S_0}{\gamma+p}=\frac{0.022}{0.04+0.088}=0.545$

となります。これから削減率 $q=1-0.022/0.128=0.829$ 、隔離率 $p=0.088$ 、 $1/p=11.4$ が分かります。これらを用いると次の計算結果を得ます。

謝辞
最後に、本サイトでまとめた感染シミュレーションは、３月末からわたって検討してきました。その間多くの知人に補遺を除く本サイトを見ていただき、コメントをいただき、その都度見直すことができて、大変感謝申し上げております。本来ならばお名前を挙げて謝意を表すべきところですが、まだ最終原稿となっていない今の段階では、文責を筆者のみに帰すため、遠慮させていただきます。

5/17/2020 梶原宏之　kajiwara@cacsd2.sakura.ne.jp

追記
有料サイトのため、全文を読むことができませんが、次の発表がありました。 $\gamma$ の値は0.04（感染性期間25日）は長すぎるかもしれません。

「コロナ感染は発症前から発症5日が最多、6日以後はほとんど感染しない！」
台湾では2020年5月14日現在、新型コロナウイルス感染症（COVID-19）患者440人、死者7人で感染制御に成功しています。その台湾衛生福利部疾病管制署（台湾CDC）からJAMA Internal Medicine（2020年5月1日オンライン版）に「台湾のCOVID-19感染ダイナミクス」の論文が出ました。100人の感染確定者と2,761人の接触者について追跡（contact tracing）を行い、時期による感染力の違いが明らかになりました。最重要点は次の3つです。接触者追跡は発症4日前までさかのぼれ、COVID-19は発症前から発症5日までに感染リスクがあり、6日以降はない！COVID-19の二次感染までの期間（serial interval）は４～５日、SARSは8.4日。
西伊豆健育会病院病院長　仲田和正
2020年05月14日 16:50