はじめに
新型コロナウイルス感染症(COVID-19)の感染拡大の予測や分析において、SEIRモデルあるいはSIRモデルが良く使われますが、最近スケールフリーネットワークに基づくG.M. Szaboによる論文やY. Ohsawaによる論文が発表されて、注目されています。
スケールフリーネットワークとは、簡単に言えば、下図のように「一部の頂点が他のたくさんの頂点と辺で繋がっており、大きな次数を持っている一方で、その他の大部分はわずかな頂点としか繋がっておらず、次数は小さいという性質」(複雑ネットワーク)を持つネットワーク構造ですが、この性質がクラスター班の調査で見出されている少数の人の基本再生産数が大きい一方で、大多数の人の基本再生産数は小さいという性質(新型コロナクラスター対策専門家（２）)によく似ており、親和性が高いようです。
（G.M. Szaboによる論文から引用）
そこで、本記事では、この研究にインスパイアされて、**個々の基本再生産数が裾の広い分布に従って決まるとき全体としての感染拡大にどう影響するか？**をSEIRモデルを拡張して検証してみたいと思います。
用語の定義
SEIRモデルに基づくシミュレーションを行います。まず、集団と個人を定義します。
- N:対象集団の人数
- j:対象集団の個人j
次に、個人jに関する確率変数を定義します。
- $s_j$:個人jが"感染症に対して免疫を持たない者(Susceptible)"である確率
- $e_j$:個人jが"感染症が潜伏期間中の者(Exposed)"である確率
- $i_j$:個人jが"発症者(Infected)"である確率
- $r_j$:個人jが"感染症から回復し免疫を獲得した者(Recovered)"である確率
確率なので当然以下を満たします。
- $\forall j, 0 \leq s_j, e_j, i_j, r_j \leq 1$
- $\forall j, s_j + e_j + i_j + r_j = 1$
次に全体としての変数を定義します。
- S:感染症に対して免疫を持たない者(Susceptible)の集団としての割合
- E:感染症が潜伏期間中の者(Exposed)の集団としての割合
- I:発症者(Infected)の集団としての割合
- R:感染症から回復し免疫を獲得した者(Recovered)の集団としての割合
また、パラメータとして、次を導入します。
- $R_{0j}$:個人jの基本再生産数
- lp:潜伏期間（5日を採用）
- ip:感染期間（8日を採用）
- $\beta_j$:個人jの感染率($\beta_j = R_{0j}/(N * ip)$)
- $e_r$:感染源が換気の悪い環境にいた確率
基本再生産数と過分散
基本再生産数$R_0$は「（感染初期において）一人の感染者がその感染期間に生み出す２次感染者数」として定義されます。感染初期としているのは、時間の経過によって、SEIRの割合が変化することなどにより、その値が変動するためです。時間に依存させる場合は$R_t$と表記します。厚生労働省のクラスター対策班の調査によれば、基本再生産数は次のグラフにようにサンプル集計されています。
先ほどのパラメータ$e_r$はこのグラフ上で換気の悪い環境にいた確率を表します。このグラフから示されるように、少数の人の基本再生産数が大きい一方で、大多数の人の基本再生産数は小さいという性質が確認できます。この基本再生産数の分布は以下の近似コードで実現します。
def COVID19R0(er):
if np.random.rand() < er:
# good environment
if np.random.rand() < 0.8:
R0 = 0
else:
R0 = np.random.randint(1,4)*1.0
else:
# bad environment
R0 = np.random.randint(0,12)*1.0
return R0
また、基本再生産数は定義から、時間の次元と感染者数の次元の両方に意味を持つことに注意が必要です。例えば、
- 一人の人が8日間かけて、一日1人ずつ感染させ続けた。
- 一人の人が1日に8人感染させ、その後２次感染を起こさなかった。
いずれの場合も$R_0=8$ですが、1はどちらかというと個人の持つ感染力の強さを暗示しており、2はどちらかというと人々が３密など同じ場面にいたことの感受性の強さを暗示していると言えるでしょう。
さて、基本再生産数の分布に基づいて、SEIRモデルを計算する数理モデルを立てる訳ですが、次の２種類のモデルを考えます。
- 感染力の過分散を想定するモデル。過分散を発症者側の要因とみなす。
- 感受性の過分散を想定するモデル。過分散を感受者側の要因とみなす。
数理モデル１（感染力の過分散）
感染力の過分散を想定するモデルは次のように書けます。
\begin{eqnarray}
\beta_j &=& \frac{R_{0j}}{ip N} \\
\frac{ds_j}{dt} &=& - s_j \cdot \sum_{k \neq j} \beta_k i_k \\
\frac{de_j}{dt} &=& -\frac{ds_j}{dt} - \frac{1}{lp} e_j \\
\frac{di_j}{dt} &=& \frac{1}{lp} e_j - \frac{1}{ip} i_j \\
\frac{dr_j}{dt} &=& \frac{1}{ip} i_j
\end{eqnarray}
特に２行目が肝心ですが、この式を$j$に関して和とると、
\begin{eqnarray}
\frac{dS}{dt} &=& - S \cdot BI + \sum_j s_j i_j \beta_j \\
\frac{dS}{dt} &=& \sum_j \frac{ds_j}{dt} \\
BI &=& \sum_k \beta_k i_k
\end{eqnarray}
となります。$\forall k, \beta_k = \beta$の場合、
\begin{eqnarray}
\frac{dS}{dt} &=& - \beta S \cdot I + \sum_j s_j i_j \beta_j \\
I &=& \sum_j i_j \\
\end{eqnarray}
となり、通常のSEIRモデルと比べ第２項だけ違いますが、この項は自分自身が自分を２次感染者にさせることはないことを表しており、SEIRモデルの推移を考えるとほぼ$s_j i_j=0$と言えるので無視できます。つまり、$\beta_j$を均一化すると通常のSEIRモデルに帰着されるモデルです。
数理モデル２（感受性の過分散）
感受性の過分散を想定するモデルは次のように書けます。
\begin{eqnarray}
\beta_j &=& \frac{R_{0j}}{ip N} \\
\frac{ds_j}{dt} &=& - s_j \beta_j \cdot \sum_{k \neq j} i_k \\
\frac{de_j}{dt} &=& -\frac{ds_j}{dt} - \frac{1}{lp} e_j \\
\frac{di_j}{dt} &=& \frac{1}{lp} e_j - \frac{1}{ip} i_j \\
\frac{dr_j}{dt} &=& \frac{1}{ip} i_j
\end{eqnarray}
２行目に着目すると、先ほどの数理モデル１に対して、$\beta_k$が$\beta_j$に移動していることが分かります。このモデルも先ほどと同様、$\beta_j$を均一化すると通常のSEIRモデルに帰着されるモデルです。
Pythonで計算してみる
では、Pythonを使って、数理モデル１と数理モデル２を計算してみましょう。
まず、計算の前提を下記とします。
- $N=200$の集団を考える。
- 初期状態は、各個人が確率1%でI（感染者）、確率99%でS（感受性者）である状態。
- 基本再生産数の平均値は$E[R_0]=2.5$付近。
ライブラリをインポートします。
import numpy as np
import matplotlib.pyplot as plt
基本再生産数を計算して分布を生成する関数です。今回は、平均値が$E[R_0]=2.5$付近となるように$er=0.6$を選択しました。また、$R_{0j}$が$j$について降順になるようにソートしています。
def COVID19R0(er):
if np.random.rand() < er:
# good environment
if np.random.rand() < 0.8:
R0 = 0
else:
R0 = np.random.randint(1,4)*1.0
else:
# bad environment
R0 = np.random.randint(0,12)*1.0
return R0
N = 200
er = 0.6
R0 = np.array([COVID19R0(er) for i in range(N)])
R0 = np.sort(R0)[::-1]
plt.hist(R0,bins=10)
plt.title("Ro distribution with E[Ro]: {}".format(np.average(R0)))
plt.xlabel('R0')
plt.ylabel('num. of people')
plt.show()
常微分方程式を解く関数です。
def my_odeint(deq, ini_state, tseq, keys):
sim = None
v = np.array(ini_state).astype(np.float64)
dt = (tseq[1] - tseq[0])*1.0
c = 0
for t in tseq:
dv = deq(v,t, keys)
v = v + np.array(dv) * dt
if sim is None:
sim = v
else:
sim = np.vstack((sim, v))
pg = 100.* t / tseq[-1]
if pg >= c:
c = c + 10
print("progress: {}%".format(pg))
return sim
数理モデル１の計算用関数です。
#define differencial equation of seir model
def seir_eq10_1(v, t, keys):
N = keys['N']
b = keys['b']
lp = keys['lp']
ip = keys['ip']
#
ds = np.zeros(N)
de = np.zeros(N)
di = np.zeros(N)
dr = np.zeros(N)
#
BI = np.sum([ b[k]*v[2][k] for k in range(N)])
#
for j in range(N):
s = v[0][j];
e = v[1][j];
i = v[2][j];
r = v[3][j];
#
ds[j] = - s * (BI - b[j]*i) # 自己感染はしない
de[j] = - ds[j] - (1/lp) * e
di[j] = (1/lp)*e - (1/ip) * i
dr[j] = (1/ip)*i
return [ds, de, di, dr]
数理モデル２の計算用関数です。
#define differencial equation of seir model
def seir_eq10_2(v, t, keys):
N = keys['N']
b = keys['b']
lp = keys['lp']
ip = keys['ip']
#
ds = np.zeros(N)
de = np.zeros(N)
di = np.zeros(N)
dr = np.zeros(N)
#
CI = np.sum([ v[2][k] for k in range(N)])
#
for j in range(N):
s = v[0][j];
e = v[1][j];
i = v[2][j];
r = v[3][j];
#
ds[j] = - s * b[j] * (CI - i) # 自己感染はしない
de[j] = - ds[j] - (1/lp) * e
di[j] = (1/lp)*e - (1/ip) * i
dr[j] = (1/ip)*i
return [ds, de, di, dr]
表示用の関数です。
def showSim_01(n, sim, keys):
N = keys['N']
plt.rcParams["font.size"] = 12
fig, ax = plt.subplots(figsize=(10,5))
x = range(N)
ax.plot(x,sim[n+0]) # extract S
ax.plot(x,sim[n+1]) # extract E
ax.plot(x,sim[n+2]) # extract I
ax.plot(x,sim[n+3]) # extract R
ax.grid(which='both')
ax.legend([ 'Susceptible','Exposed', 'Infected', 'Recoverd'])
ax.set_xlabel('person no.')
ax.set_ylim(0,1)
plt.show()
def showSim_02(t, sim, keys):
N = keys['N']
plt.rcParams["font.size"] = 12
fig, ax = plt.subplots(figsize=(10,5))
n = len(t)
sdat = [np.average(sim[4*i+0]) for i in range(n)]
edat = [np.average(sim[4*i+1]) for i in range(n)]
idat = [np.average(sim[4*i+2]) for i in range(n)]
rdat = [np.average(sim[4*i+3]) for i in range(n)]
ax.plot(t,sdat) # extract S
ax.plot(t,edat) # extract E
ax.plot(t,idat) # extract I
ax.plot(t,rdat) # extract R
ax.grid(which='both')
ax.legend([ 'Susceptible','Exposed', 'Infected', 'Recoverd'])
ax.set_xlabel('date')
ax.set_ylim(0,)
plt.show()
シミュレーション用の関数です。初期状態は、各個人で均一です。
def calcsim( R0, keys):
# solve seir model
N = keys['N']
# S E I R
ini_state=[np.ones(N), np.zeros(N), np.zeros(N), np.zeros(N)]
ini_state=[np.ones(N)*0.99, np.zeros(N), np.ones(N)*0.01, np.zeros(N)]
t_max=180
dt=0.01
t=np.arange(0,t_max,dt)
keys['b'] = R0/(N*keys['ip'])
#
sim = my_odeint(seir_eq10, ini_state, t, keys)
#
return sim,t
数理モデル１を実行して表示するコードです。非常に重いです(-_-;)。
seir_eq10 = lambda v, t, keys: seir_eq10_1(v, t, keys)
keys = {'N':N, 'lp':5, 'ip':8, 'R0':2.5 }
sim,t = calcsim(R0, keys)
showSim_01(np.int(180/0.01*0.95), sim, keys)
showSim_02(t, sim, keys)
数理モデル２を実行して表示するコードです。
seir_eq10 = lambda v, t, keys: seir_eq10_2(v, t, keys)
keys = {'N':N, 'lp':5, 'ip':8, 'R0':2.5 }
sim,t = calcsim(R0, keys)
showSim_01(np.int(180/0.01*0.95), sim, keys)
showSim_02(t, sim, keys)
シミュレーション結果
それでは、計算結果を見てみましょう。
基本再生産数の分布
基本再生産数の分布は次のようになっています。基本再生産数の平均値は$E[R_0]=2.49$でした。
数理モデル１（感染力の過分散）
まずは、数理モデル１で横軸を各個人、縦軸を95%時間経過後の(s,e,i,r)の各確率としたグラフです。
何と、**感染力の過分散にもよらず個人差は殆ど現れない！**結果となりました。
次に、通常のSEIRモデルのように、集団としてのSEIRの推移を見てみます。
最終的に、およそ90%が感染し、集団免疫により終息する結果となります。
数理モデル２（感受性の過分散）
次に、数理モデル２で横軸を各個人、縦軸を95%時間経過後の(s,e,i,r)の各確率としたグラフです。
先ほどの結果とは大きく異なり、**個人のR0分布の影響が非常に強く現れる！**結果となりました。$R_0$が大きい個人ほど感染しやすく、$R_0$が小さい個人には段々感染しにくくなる傾向が見られます。
次に、集団としてのSEIRの推移を見てみます。
最終的に、**およそ37%が感染し、集団免疫により終息する結果!!**となります。
数理モデル１で計算した結果が90%なので、実に感染者数が59%も抑制されました。
考察
以上から、基本再生産数の過分散を考慮したSEIRモデルによる計算に関して、次の傾向がシミュレーションから導けます。
- 基本再生産数の過分散を感染力の過分散としたモデル（数理モデル１）では、基本再生産数が均一とした結果とほぼ変わらない。
- 基本再生産数の過分散を感受性の過分散としたモデル（数理モデル２）では、基本再生産数が均一とした結果と大きく異なり、最終的な感染者が59%も抑制される。
- スケールフリーネットワークとのアナロジーで言えば、スーパースプレッダーが問題なのではなく、個人がどれだけ多くの人と接触機会を持ったかが重要である。
さらに言えば・・・
- 日本で感染者数が比較的抑制された要因は、早期の３密機会回避徹底、マスク・手指衛生の習慣による個人の感受性の低減、ではないでしょうか。
- 逆に言えば、欧米で非常に感染が広がってしまった要因は、感受性の低減策が有効でなかったと言えます。ハグをする習慣、マスクを嫌う風習、宗教施設で定期的に集会を行う習慣、など様々な要因が考えられます。
- 計算があまりに重いので$N$を増やしすぎるのはお勧めしません。専ら速いと評判のJuliaに移行した方がいいかも…。
- 下記のスケールフリーネットワークの個人ごとの感染可能なノード数分布と、数理モデル２のR分布が非常に似ているのは偶然でしょうか？何か繋がりがあるように思います。
（G.M. Szaboによる論文から引用）
参考リンク
下記のページを参考にさせて頂きました。