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: 一部推敲。