Pythonプログラミング（ステップ６・SIRモデル）

伝染病の流行の微分方程式によるモデル化

以下のような状況で感染が拡大するような疾病（インフルエンザ等）を考える。

1. 新型の病気のため、人々は免疫を持っていない。
2. 人口が密集しており、人々の行き来が活発で、「平均的に」不特定多数の人と接触している。
3. 人口の流入・流出はない（ロックダウンされた国や都市を想定）。
4. 接触によって、一定の割合で感染する。
5. 感染者は快復後に免疫が生じ、罹患者と再び接触しても再発はしないし、未感染者に罹患させることもない。

そして、未感染（Susceptable)者、感染(Infected)者、および快復した（Recovered）人の３種類が「一様に」混在し、 以下のようなプロセスで感染が伝搬・拡大するような状況を想定する。

ここでは、一人ひとりの状態を考えるのではなく、S, I, Rそれぞれの人数の変化として、流行の様子を捉えることにする（ゆらぎの効果は無視する）。 総人口を$N$、未感染者(S)が $x$人、感染者(I)が $y$人、快復した人(R)が $z$ 人と置く。 人口が十分に大きければ、これらを連続量とみなしても構わないだろう。 また、経過時間を変数 $t$（日）で表すことにする。

集団の中で、一人あたり毎日平均 $m$ 人と接触するとする。 そのうち感染している人の割合は $y(t)/N$ であるから、 感染者と非感染者の接触回数は1日あたり $m \frac{y(t)}{N} x(t)$ となる。 それぞれの接触毎に感染が生じる確率を $p$ とすれば、 $\Delta t$ 日の間に生じる新たな感染者は、全体として $$ m \, \frac{y(t)}{N} \, x(t) \, p \, \Delta t $$ となるだろう。 ここで、$\beta = m p / N$ と置くと、$t$日目から $t+\Delta t$日目の間に、新たな感染によって非感染者の数が $$ x(t+\Delta t) - x(t) = - \beta \, x(t) \, y(t) \, \Delta t $$ だけ変化することになる。これを $$ \frac{x(t+ \Delta t) - x(t)}{\Delta t} = - \beta \, x(t) \, y(t) $$ と書き直して、$\Delta t \to 0$の極限を取ると $$ \frac{d x(t)}{d t} = - \beta \, x(t) \, y(t) \tag{1} $$ を得る。

一方、感染者は１日あたり一定の率$\gamma$で快復してゆくとすると、快復者の増加分は $$ \frac{d z(t)}{d t} = \gamma \, y(t) \tag{2} $$ と書ける（$\gamma$は一群の感染者の数が快復によって$1/e$になるまでに要する日数の逆数*）。

人口の流入・流出は考えないので、各人口の総和は一定 $$ x(t) + y(t) + z(t) = N $$ すなわち $$ y(t) = N - x(t) - z(t) $$ であるから、未感染者の人数の変化が $$ \frac{dy(t)}{dt} = - \frac{dx(t)}{dt} - \frac{dz(t)}{dt} = \beta \, x(t) \, y(t) - \gamma \, y(t) \tag{3} $$ となる。

上記はSIRモデル^†と呼ばれ、古典的な伝染病の流行モデルのひとつとして知られており、いくつかのバリエーションも含め、その振る舞いが研究されている。

補足

　練習：方程式の無次元化

各変数を $\tilde{x} = x/N, \, \tilde{y} = y/N, \, \tilde{z} = z/N, \, \tilde{t} = m \,p\,t$ のように変数変換（無次元化）することで、 SIRモデルはパラメータをひとつしか含まないかたち $$ \begin{eqnarray} \frac{d \tilde{x}}{d \tilde{t}} & = & - \tilde{x} \, \tilde{y} \\ \frac{d \tilde{y}}{d \tilde{t}} & = & \tilde{x} \, \tilde{y} - \frac{1}{R_0} \tilde{y}\\ \frac{d \tilde{z}}{d \tilde{t}} & = & \frac{1}{R_0} \tilde{y} \end{eqnarray} $$ に変形できることを確かめなさい。

SIRモデルのシミュレーション

二次のルンゲクッタ法を使ってSIRモデルのシミュレーションを行うコードと結果の例を以下に示す。このコードでの状況設定は

1. 人口100万人($N=10^6$)で、そのうち10人が初期に感染
2. 毎日1人あたり10人に接触。1日・1人の接触あたりの感染率は0.02（$m p=0.2$)。
3. 感染者の快復日数が14日程度（$\gamma = 1/14 \approx 0.07$）。

である（これらの値そのものに強い根拠は無い）。

# coding: utf-8
import math
import matplotlib.pyplot as plt

N = 1000000
x = N-10
y = 10
z = 0

taxis=[ ]
xaxis=[ ]
yaxis=[ ]
zaxis=[ ]

beta=0.2/N
gamma=0.07
dt=0.001

t = 0
cnt=0
while t<180:
    if cnt%100==0:
        taxis.append(t)
        xaxis.append(x)
        yaxis.append(y)
        zaxis.append(z)
# step 1
    kx1 = - beta*x*y
    ky1 = beta*x*y - gamma*y
# step 2
    t2 = t+dt
    x2 = x + kx1*dt
    y2 = y + ky1*dt
    kx2 = - beta*x2*y2
    ky2 = beta*x2*y2 - gamma*y2
# update
    x = x + (kx1+kx2)*dt/2
    y = y + (ky1+ky2)*dt/2
    z = N - x - y
    t = t + dt
    cnt = cnt + 1

plt.title("SIR MODEL")
plt.plot(taxis,xaxis, color=(0.0,1,0.0), linewidth=1.0, label='S')
plt.plot(taxis,yaxis, color=(1.0,0,0.0), linewidth=1.0, label='I')
plt.plot(taxis,zaxis, color=(0.0,0,1.0), linewidth=1.0, label='R')
plt.xlim(0,180)
plt.legend()
plt.xlabel('DAY')
plt.ylabel('X, Y, Z')
plt.grid(True)
plt.show()

　練習：シミュレーション：爆発的な感染拡大の抑止

上記のコードで、初期の感染者数, beta, gamma を変えてシミュレーションし、感染の拡大・収束にはどのようなシナリオがあり得るかを予想し、 流行のピークを下げるのに効果的と考えられる対策について考察しなさい。

　ヒント：パラメータの目安と状況設定

　練習：４次のルンゲクッタ法

上記のコードを４次のルンゲクッタ法に書き直し、計算結果の違いをみてみなさい。

　解説: SciPyの利用

# coding: utf-8
import math
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

def sir_model(v, t, beta, gamma):
    x = v[0]
    y = v[1]
    z = v[2]
    dxdt = - beta*x*y
    dydt = beta*x*y - gamma*y
    dzdt = gamma*y
    return [dxdt, dydt, dzdt]

N = 1000000
beta=0.2/N
gamma=0.07
dt=0.1

init_vals=[N-10, 10, 0]
times = np.arange(0.0, 180, dt)
args=(beta, gamma)
results = odeint(sir_model, init_vals, times, args=args)

plt.title("SIR MODEL")
plt.plot(times[::10],results[::10,0], color=(0.0,1,0.0), linewidth=1.0, label='S')
plt.plot(times[::10],results[::10,1], color=(1.0,0,0.0), linewidth=1.0, label='I')
plt.plot(times[::10],results[::10,2], color=(0.0,0,1.0), linewidth=1.0, label='R')

plt.xlim(0,180)
plt.legend()
plt.xlabel('DAY')
plt.ylabel('X, Y, Z')
plt.grid(True)
plt.show()