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

投稿日時: 投稿者:

ここでは次の文献にもとづいて「感染ダイナミックス」について述べます。
[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検査を増やすことの効果を確認できたと言えます。ただし、偽陽性と偽陰性の可能性を勘案して、陽性の出る可能性の高い場合に検査を行う必要があることは述べるまでもありません。