2021年6月26日土曜日

変異株の置き換わりについて

感染力の異なる複数の株があるときの感染者数の時間変化について考えたい。 特段の人為的な介入が行われず、なすがままに感染が広がっているとするとき、感受性のある (susceptible) 者の有限性を考える必要がない程度に感染が少数なら、株 $i$ の感染者数 $x_i\;(> 0)$ の時間変化は株に固有の定数 $k_i\;(> 0)$ を用いて、

$$\frac{\mathrm{d}x_i}{\mathrm{d}t} = k_i x_i \quad\tag{1}$$

で表されるだろう。 すなわち、この株の感染者数の時間に関する動きは指数関数 $x_i = A_i \exp(k_i t)$ で表される(ただし $A_i > 0$ は $t = 0$ における感染者数)。 $k_i$ が大きな株ほど広がりは速く $k_i$ は感染力を表すとみなせる。

広く他人との接触を抑制するような社会的な介入政策一般や(日本ではあまり有効に行われていない)防疫目的の検査による感染者数の時間変化に対する感染抑制の効果は、$k_i$ と同様に概ね感染者数に比例して作用するとみなせる。 その大きさは時間に関して変動するものの、株の種類には概ね依存しないと考えられる。 この影響を $c(t)$ と書けば、微分方程式 (1) は、

$$\frac{\mathrm{d}x_i}{\mathrm{d}t} = (k_i\!-\!c(t))\,x_i \quad\tag{2}$$

と変更される。 ワクチンの効果は株により抑制効果にある程度の違いがあるかもしれないが、その効果に時間による変動はないと考えられるので、(定数分の違いは無視して)同様に $c(t)$ に繰り込まれるものとする。 よって、このとき感染者数は、

$$x_i = \frac{A_i \exp(k_i t)}{\exp(C(t))} \quad\tag{3}$$

のように表される。ただし、$C(t) = \int_0^t c(\tau) \mathrm{d}\tau$ とした。 $C(t)$ は株に依存しない時間の関数と考えられたのであったから、複数の株 $i, j$ を考えるときその感染者数の比 $x_j/x_i$ を考えれば打ち消し合って、比は単純に人為的介入の効果に影響しない指数関数、

$$\frac{x_j}{x_i} = \exp((k_j\!-\!k_i)\,t) \quad\tag{4}$$

で表せることになると期待される。 ゲノム解析や特定の変異を検出する PCR 検査が無作為と見なせるように行われ市中の感染者数を反映しているなら、(サンプル誤差を除いて)(4) はそれら各株の感染確認者数においても成立する。

実際、2 変異株の比が概ね (4) の単純な関係に従っていることは図 1 でも実証される。 図 1 は東京都健康安全研究センターが週ごとに公表している PCR による変異株の変異検出数から作成したもので[1]、2021年3月以降に起こった第 4 波時の E484K 変異をもつ変異株から N501Y 変異をもつ変異株(後者は主に英ケントで最初に流行が確認された α 変異株)の比の時間変動を片対数グラフで示している。 比は、緊急事態宣言のような社会的介入の有無や感染確認者数の大小に影響を受けたようすはなく、(4) 式が示すように概ね直線で表される。 $k_{\,\mathrm{N501Y}}\!-\!k_{\,\mathrm{E484K}}$ に相当する傾きは 0.065 日−1 で、これら変異株それぞれの感染力の違いのみから定まる変異株間の移行の速さを表している。

図 1. 東京における E484K 変異をもつ変異株と N501Y 変異をもつ変異株の比(東京都健康安全研究センター調査[1]から)

一方、イギリスでは多くのゲノム解析を行って各都市の変異株の感染拡大状況を詳しく調べている。ロンドンにおける B.1.177+ 株(類似する株の集団)から α 株、および α 株から δ 株への置き換わりのようすを同様にプロットしてみると図 2 のようになる。

図 2. ロンドンにおける B.1.177+ 株(類似する株の集団)と α 株 (B.1.1.7)、および α 株と δ 株 (B.1.617.2) の比の時間変化(Wellcome Sanger Institute による lineage データ[2]から)

図 2 では、比の対数の時間変化は概ね直線に従うものの、途中からは反れてより緩やかな拡大となる傾向があるように見てとれる。 その理由は明らかではないが、地域的な人々のネットワークによる局所的な susceptible の有限性が効いているのかもしれない。 ただし、イギリスの数学者の Alex Selby 氏は Twitter 上で COG-UK (COVID-19 Genomics UK Consortium) のデータを用いたイギリスにおけるこの比が見事に指数関数に従うことを示している[3]。 Selby 氏は、イギリスにおける α 株から δ 株への置き換わりの傾き ($k_\delta\!-\!k_\alpha$) が安定して 0.11 日−1 であるとしている。 これは、以下で見る東京における δ 株への置き換わりの速さとおよそコンシステントである。

現在、首都圏および大阪などで進行中の α 変異株(B.1.1.7, イギリス株, ケント株)からより感染力の強い δ 変異株 (B.1.617.2, インド株)への置き換わりについて上の知見を具体的に当てはめてみたい。 しかし、それにはいくつかの困難がある。 日本においては、新型コロナのゲノム解析は感染研で行われているが、時間がかかる上にごく少数であり、報告は特定の都市別には公表されていない。 どのような株が存在するか絞られている場合には、代わりにそれぞれの変異株に特徴的な特定の変異(α 株の場合 N501Y, δ 株では L452R)を検出する自治体等による PCR 検査が参考となる。 しかし、一般には市中の変異株の割合を知るために重要な無作為抽出を考慮していない結果が報告されているようであり、これは特にクラスターが発生しそれらが変異検出の検査にかけられた場合、結果を歪ませることになる。 さらに検出された変異をもつ検体の数だけが報告され、いくつの検体を検査した結果なのか、サンプルサイズがわからない報告を行っている自治体も多いようだ。 また、最近の厚労省の指示により、主に α 株の数を知るために行われていた N501Y 変異検出のための検査を多くの自治体がやめてしまった。

ここでは、少なくとも検査数を毎日公表している東京都の調査に基づいた分析のみを取り上げる。 東京都も N501Y の検査をやめ L452R の検査結果のみを公表しているので、δ 株(正確には L452R 変異を持つ株, 以下 δ 株と同一視する)以外の株はすべて α 株(同様に正確には N501Y 変異のある株)なのだと仮定する。 こととき求めるべき比は、δ 株の存在比 $p$ のオッズ $p/(1\!-\!p)$ ということになる。 東京都新型コロナウイルス感染症対策本部の本部報[4]を 2021 年 6 月 1 日から 25 日まで集計した結果は、図 3 左パネルのようになる。

図 3. 東京における L452R 変異をもつ株(主に δ 株)のオッズの対数(ロジット)とそれへの直線の当てはめ(左). 同じものを割合で表したもの(右).(東京都新型コロナウイルス感染症対策本部報[4]より)

図 3 でも概ね置き換わりが (4) で示される指数関数に従っていることがわかる。 ただし、用いているデータは本来、無作為抽出を意図したものではないため、偏りが存在すると思われる。 特に、クラスターが発生し、クラスターの集団に含まれる検体が複数変異の検査を受けた場合、結果を過大に評価することになる。 市中での δ 株の割合を知りたいというここでの目的においてはこれは不適切なため、図 3 では報道を元にクラスターが含まれると思われる結果は修正している。しかし、想定より上に飛び出しているものには依然こうした互いに関連した検体の結果が含まれているかもしれない。

なお技術的なこととして、図 3 は、扱っている範囲において依然置き換わりの初期であって δ 変異株の割合は小さい。 またサンプルサイズも小さな場合があるため、サンプル誤差を与える二項分布の歪が大きく、対数をとって最小二乗法で回帰直線を求めることはあまり適切ではないと思われる。 このため、二項分布の歪はなるべく保ったままモデルの直線を与えるパラメーター空間を網羅的に走査してそれらの尤度を求めている。 ここでは、これらの手続きの数学的詳細は省く。

この結果を用いると、変異株の割合の今後の変化だけでなく、感染者数のある程度の予測を行うことができる。 これについては項を改めて論じる。

  • 2021-06-26: 図 1 を差し替えました。
  • 2021-06-28: 図 2 を差し替え、関連する記述を一部修正しました。
  • 2021-06-29: 一部推敲。

2020年5月26日火曜日

正単体のデカルト座標

シミュレーションのため、$n$ 次元ユークリッド空間 $\mathbf{R}^n$ において互いに等距離離れた $n+1$ 個の点の座標が欲しかったので、備忘のため記録。 自己流。調べればもっとよいものがあるかもしれない。

上のような点を結べば、正三角形、正四面体の拡張である $n$ 次元正単体(正 $n$ 単体)をなす。 $n$ 次元に 1 次元増やした $\mathbf{R}^{n+1}$ を考えると、その標準基底 $$\begin{array}{c} (1, 0, \cdots, 0),\\ (0, 1, \cdots, 0),\\ \vdots\\ (0, 0, \cdots, 1)\, \end{array}$$ が、こうした正 $n$ 単体をなしており、どの点も互いに $\sqrt{2}$ だけ離れていることはすぐにわかる。

よって、法線 $(1, 1, \cdots, 1)$ に直交する $n$ 次元部分空間への射影が、求める(頂点間の距離が $\sqrt{2}$ の)正 $n$ 単体のデカルト座標を与える。 以下の $(x_{ij})$ $(i = 0, ..., n,\ j = 0, ..., n-1)$ の行はそうしたものの例となる(添字を $0$ からとっていることに注意)。 $$(x_{ij}) = \left( \begin{array}{ccccccc} \displaystyle \frac{-1}{\sqrt{2}}, & \displaystyle \frac{-1}{\sqrt{6}}, & \displaystyle \frac{-1}{2 \sqrt{3}}, & \cdots, & \displaystyle \frac{-1}{\sqrt{(j+1)(j+2)}}, & \cdots, & \displaystyle \frac{-1}{\sqrt{n (n+1)}} \\ \displaystyle \frac{1}{\sqrt{2}}, & \displaystyle \frac{-1}{\sqrt{6}}, & \displaystyle \frac{-1}{2 \sqrt{3}}, & \cdots, & \displaystyle \frac{-1}{\sqrt{(j+1)(j+2)}}, & \cdots, & \displaystyle \frac{-1}{\sqrt{n (n+1)}} \\ \displaystyle 0, & \displaystyle \sqrt{\frac{2}{3}}, & \displaystyle \frac{-1}{2 \sqrt{3}}, & \cdots, & \displaystyle \frac{-1}{\sqrt{(j+1)(j+2)}}, & \cdots, & \displaystyle \frac{-1}{\sqrt{n (n+1)}} \\ \displaystyle 0, & \displaystyle 0, & \displaystyle \frac{\sqrt{3}}{2}, & \cdots, & \displaystyle \frac{-1}{\sqrt{(j+1)(j+2)}}, & \cdots, & \displaystyle \frac{-1}{\sqrt{n (n+1)}} \\ \vdots & \vdots & \vdots & \ddots & \vdots & & \vdots & \\ \displaystyle 0, & \displaystyle 0, & \displaystyle 0, & \cdots, & \displaystyle \sqrt{\frac{j+1}{j+2}}, & \cdots, & \displaystyle \frac{-1}{\sqrt{n (n+1)}} \\ \vdots & \vdots & \vdots & & \vdots & \ddots & \vdots & \\ \displaystyle 0, & \displaystyle 0, & \displaystyle 0, & \cdots, & \displaystyle 0, & \cdots, & \displaystyle \sqrt{\frac{n}{n+1}} & \end{array} \right)$$ すなわち、 $$x_{ij} = \left\{ \begin{array}{cl} \displaystyle \frac{-1}{\sqrt{(j+1)(j+2)}} & \hspace{3em} (i \leq j)\\ \displaystyle \sqrt{\frac{j+1}{j+2}} & \hspace{3em} (i = j + 1)\\ \displaystyle 0 & \hspace{3em} (i > j + 1)\\ \end{array} \right.$$ 正単体の重心(各座標値の平均)は原点となる。

法線ベクトル $(1, 1, ..., 1)$ の列を加えた $n+1$ 次元正方行列 $$\left( \begin{array}{rrrcrr} -1 & -1 & -1 & \cdots & -1 & 1 \\ 1 & -1 & -1 & \cdots & -1 & 1 \\ 0 & 2 & -1 & \cdots & -1 & 1 \\ 0 & 0 & 3 & \cdots & -1 & 1 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & n & 1 \\ \end{array} \right)$$ の列が互いに直交することはすぐに確かめられ、この各列を正規化して直交行列とすれば、元の標準基底のベクトルを行ベクトルとして右から施し移したものはその行そのものとなるので(法線方向成分を除いて)上が得られる。 ただしこの取り方では行列式は $1$ とは限らず(特殊直交行列とは限らず)反転があるかもしれない。

2020年4月26日日曜日

観測された感染者数に関する覚書き 3

前記事で用いた $I_E(t\,;t_P)$ および $c_E(t\,;t_P)$ を、韓国・ドイツ・日本の 3 か国およびアメリカ・ニューヨーク州について示し、それぞれについて検討する。 データにはジョンズ・ホプキンズ大学 CSSE が公開している 2020 年 4 月 23 日までの時系列データを用いた(差分を取るため 4 月 22 日が $t_P$ となる)[1]

韓国

韓国においては 2020 年 2 月 18 日に大邱(テグ)市で宗教団体に関係する感染者が見出だされ、翌 19 日に集団感染が発生していることが把握された。 その後、宗教団体に関係する者に関しては全数検査が行われた。 3 月 8 日の時点において、この集団感染による感染者が 4500 件近くあり、韓国の感染者の 6 割以上を占めていることが韓国における感染の特徴と言える[2]。 長期的に見れば、3 月後半以降、現在までに 2500 人程度の新たな感染者が見つかっており[1]、それもよく防がれているように見える。 しかし全体としてみれば発見の量は初期の動きに比べ小さい。 よって、これらはここに示した感染拡大初期に関して $I(t)$ および $c(t)$ の比較的よい近似を与えているものと思われる。

図に示されていない範囲で起こったと思われる大規模な集団感染時には、$I(t)$ の名目上の倍加時間は短かったといえるかもしれない。 しかし、それは飽くまで集団感染に起因するものである。 大邱での大規模な検査が始まった時点での観測された感染者数 $H(t)$(上図灰色の実線)の急激な上昇は多数の検査を行った帰結であり、隠れた感染者 $I(t)$ はその時点ですでに数千人いたとするのがもっともらしい。 仮にこの時期に感染自体の倍加時間 $D$ が 2-3 日のように短かい期間で拡大したのだとしても、上図のように検査と隔離の効果で増加は抑えられており、後に徐々に減少したのではないかと思われる。 実際には、初期の隠れた感染者 $I_0$ がもっと多く、$D$ がもっと長かったとして矛盾はない。 初期の韓国においては限定された集団を無症状者も含めて網羅的に検査できたことが、その後の $I(t)$ の減少に大きく寄与したと考える。 徹底した早期の検査がその後の感染拡大をうまく喰い止めた例は、規模は異なるが日本の和歌山県における病院での集団発生後の状況も事例として挙げられるだろう。

検査率を上から抑える $c_E(t\,;t_P)$ のグラフは概ね高い値を示し、倍加時間が長い場合も、初期の検査で隠れた感染者を 1 日あたり数 % のオーダーで見つけ出していたものと思われる。 $I_0$ の見積もりと合わせれば、最初の記事で述べた $A = c I_0$ はオーダーとして初期に 100 人/日程度であったろう。 例として、大邱の集団感染発覚当初からしばらくたった状況を考え、$H^* = 1000$ 人とすれば、2 月 25-26 日にこの基準を超える。 一方、2 月 23-28 日の見かけの累積確定例時定数は約 3.7 日(倍加時間約 2.6 日)の依然として短いものであった。 これより仮に $A/H^* \approx 0.1$ として、$d(\log \tilde{H})/dt = \lambda + A/H^*$ の関係を用いれば、実際の時定数は $1/\lambda \approx 6$ 日(倍加時間約 4 日)となり、検査によるの見かけの効果が実際に働いていたであろうことが裏付けられる。

ドイツ

ドイツで本格的に検査が拡大したのは 2020 年 2 月 26 日頃だが、政府が拘束力のあるロックダウンの措置を行ったのは 3 月 23 日だった[3]。 検査数は 3 月半ば頃より拡大され、1 日あたり 5 万件の大規模な検査が行われてきた[4]。 3 月下旬からは無症状者への検査も行われた[3]。 観測された感染者数(累積確定例)は、3 月中旬ごろまでおおむね倍加時間 3 日ほどであった。

極端に短くない倍加時間 ($D \geq 4$) を考えたとき、グラフから、ドイツにおいて観測された累積の感染者数 $H(t)$(灰色)が隠れた感染者数の下限 $I_E(t\,;t_P)$(実線)を超えたのは、ロックダウンが発令された 3 月下旬に入ってからであることがわかる。 すなわち、そのころには少なくともまだ把握された感染者と同等以上の隠れた感染者がいたことになる。 仮に、3 月 18 日時点で、観測された感染者の 5 倍の隠れた感染者がいたとしたとき、一定の倍加時間のもとそれ以前の感染者数がどのようであったかを破線で表している。 前出の $I_R$ に相当するこの追加の感染者は単純に指数関数的に増大するので、この程度の大きさであった場合、下限 $I_E$ との初期の比はずっと小さくなる。 よって、これらの時点で韓国ほどには初期換算感染者を削りきれていないであろうドイツにおいても、検査初期の 2 月 25 日に数千程度の隠れた感染者がいたとすれば、見かけと異なって倍加時間 4 日以上で拡大していたとして矛盾はない。

検査率の最大値 $c_E$ は、3 月初めごろまでは韓国より小さな値を示しており、長い倍加時間を想定した場合、1 日 1 % に満たないものだった。 よって、$A$ の効果もこのころには韓国より 1 桁ほど小さかったかもしれない。 その後、$c_E$ は急激に上昇しており、実際の検査率も同様の傾向があったと思われる。 このことは、ドイツの初期の倍加時間が韓国ほど短くなかったことを説明するかもしれない。

初期の倍加時間を 6 日とした場合、ロックダウン宣言の 3 月 23 日までに削った初期換算感染者数 $J(t)$(起点 2 月 26 日)は約 3100 人、倍加時間 3 日とした場合、約 490 人である(前記事の図参照)。 これは、広範な検査と隔離とが行われなかった場合、3 月 23 日までの 26 日間で、それぞれ約 6 万 2 千人と約 20 万人に拡大したはずのものとなり、判明した感染者の隔離が適切に行われているなら大規模な検査には隠れた感染者の削減に大きな効果があったことがわかる。

ニューヨーク州

大規模な感染拡大がみられたニューヨークでは、最初の確定例が見つかった 2020 年 3 月 1 日(JHU CSSE データでカウントされたのは翌 2 日)以降、検査による累積確定例の急激な上昇が見られた。 3 月 5-19 日の見かけの倍加時間は 1.8 日であった。 一方、2020 年 4 月 23 日のニューヨーク・タイムズは、3 月 1 日の時点で実際には(ニューヨーク市のみで)1 万人を超える隠れた感染者がいたのだとするノースイースタン大学の研究グループの研究を紹介している[5]。 4 月 25 日現在、研究の詳細は明らかでないが、この値は倍加時間がある程度長いとした場合の隠れた感染者の最小値 $I_E$ と矛盾しない。 図では $I_E$ が 1 万より小さい倍加時間 6 日以下の場合に関して、3 月 1 日の時点で 1 万人の隠れた感染者がいたとした場合のその後の隠れた感染者数を破線で示している。

日本

観測された感染者数 $H(t)$ や、それによる初期換算感染者数 $J(t)$ で特異的な少なさを示している日本においては、$I_E$ および $c_E$ のグラフも特異的なものとなっている。 パラメーターの倍加時間にあまり依存せず、隠れた感染者の最小値 $I_E$ は $H(t)$ に沿って増加し、検査率の最大値 $c_E$ は時間に依存しない大きな一定の値を示している。 検査の強い抑制によって、これらは実際の値と大きく離れたものでありうる。 しかし、検査の強い抑制が一定の基準で行われ続けていると考えるなら、$c(t)$ がほぼ定数となることと、$H(t)$ のほぼ一定の時定数が隠れた感染者 $I(t)$ のそれを反映していることとは想定しうるだろう。

2 月 15 日から 30 日間の倍加時間はおよそ 7 日であり、これが初期の $I(t)$ の倍加時間も反映していると思われる。 このときの $I_0$ の最小値、すなわち検査が削った初期換算感染者数は 240 人程であるが、4 月下旬も隠れた感染者が少なからず存在していることを考えれば、起点とした 2 月 14 日の感染者はこれよりはるかに多かった。 図では、$I_0$ が 1000 人とした場合と 10000 人とした場合のその後の動きを破線で示している。 隠れた感染者は検査からほとんど影響を受けることなく、ほぼ指数関数的に増大する。 倍加時間 7 日が続いていたとした場合、これらの $I_0$ による 4 月 1 日(47 日後)の隠れた感染者数は、それぞれ約 8 万 6 千人と、約 103 万人となる。

*   *   *

結局、観測された感染者数の動きを隠れた感染者数と正しく区別して論じるならば、日本の観測された感染者数の少なさを、日本の隠れた感染者数がヨーロッパ諸国で想定される値よりもずっと小さいとする根拠とすることはできないことがわかる。

「新型コロナウイルス感染症対策専門家会議」は2020 年 4 月 1 日の資料においていわゆる「オーバーシュート」に関し「我が国では、今のところ諸外国のような、オーバーシュート(爆発的患者急増)は見られていない」として、次のような注釈を加えた[6]

オーバーシュート: 欧米で見られるように、爆発的な患者数の増加のことを指すが、2~3 日で累積患者数が倍増する程度のスピードが継続して認められるものを指す。異常なスピードでの患者数増加が見込まれるため、一定期間の不要不急の外出自粛や移動の制限(いわゆるロックダウンに類する措置)を含む速やかな対策を必要とする。(後略)
さらに太字とアンダーラインで強調し、
いわゆる「医療崩壊」は、オーバーシュートが生じてから起こるものと解される向きもある。しかし、新規感染者数が急増し、クラスター感染が頻繁に報告されている現状を考えれば、爆発的感染が起こる前に医療供給体制の限度を超える負担がかかり医療現場が機能不全に陥ることが予想される。
とした。

すなわち専門家会議によれば「オーバーシュート」は「欧米」で見られたような単に倍加時間 2-3 日の継続的拡大を意味するとし、このような短い倍加時間による増大がなくとも、医療のリソース不足による機能不全が迫っていることに警鐘を鳴らした。 しかし、この小論で明らかとしたように、ヨーロッパ諸国やアメリカの都市で初期に見られたこのような短い倍加時間の「爆発的急増」は大規模な検査を行うことによってすでにいた隠れた感染者を発掘した見かけの効果だろう。 日本が検査数を絞っているが故に、日本においてはそのような増加は見えず、逆に日本でも諸外国でもそれよりはやや長い一定の倍加時間での指数関数的拡大を継続的に続けている。 この点で日本の状況は何ら諸外国と変わる点を見出だせない。 違いは、日本の倍加時間が幸いにもヨーロッパなどよりやや長いと思われることと、逆に日本がほとんど検査を行っていないために感染者を適切に削っていくことができなかったことである。 指数関数的拡大を続ける限り、倍加時間の多少の長短にかかわりなく医療崩壊のような同等の帰結がもたらされるのは当然である。

「オーバーシュート」(overshoot) の語は、本来、なにがしかの量に対して、他の要因から想定したなにがしかの基準が先にあり、実際の量がそれを超えていくことを表すのに相応しい語である。 そうした他の要因が見えない「爆発的急増」もしくは「倍加時間 2-3 日」を単に表すとした用法は甚だ奇妙にみえる。 ここからは想像であるが、むしろその基準は上の言と逆に専門家会議こそが医療崩壊を想定していたのではないか。 倍加時間 2-3 日という事態が実在のものと思い込み、そのような事態になれば否応なくすぐさま医療崩壊に至ると考えたとするならば、「オーバーシュート」という語を用いた意味が理解できる。 実際には、そのような増加がなくとも隠れた感染者数は長めの倍加時間での指数関数的増大を続け、医療崩壊が想定される状況まで至ったために、医療崩壊と関係しない「オーバーシュート」=「爆発的増加」という本来の語の意味と合わない置き換えが成立したのではないだろうか。 もしそうだとすれば、専門家会議は自らが立てた不適切な道標によって自ら勝手に道を誤っていることとなる。