レイ・フロンティア株式会社のデータアナリストの齋藤です。前回の記事の続きを書いていこうと思います。
年齢構造をもつSIRモデルにおける\(\mathcal{R}_0\)
前記事では、感染症の常微分方程式モデルにおいて、感染者が生産する2次感染者の人数を表す基本再生産数\(\mathcal{R}_0\)が、方程式系の漸近挙動や非自明な平衡解の存在をきめる閾値になっていることを見てきました。次に、年齢構造を導入した場合の感染症モデルにおいても基本再生産数が重要な指標となっていることを、かなり駆け足になってしまいますがご紹介します。
モデルの構成
現実の伝染病の流行仮定を考える際には、年齢によって感染率が大きく異なるという事情から、年齢構造を導入することでより現実味のあるモデルを考えることができると期待できます。さらに、疫学的な問題として、ワクチン接種をどの年齢層に行うべきかを考えるうえでも年齢構造の考察は欠かせません。そこで、前節で考えたSIR型感染症モデルに年齢構造を入れた以下のモデルを考えてみます:
$$\left\{ \begin{align} &\frac{\partial S(t,a)}{\partial t} + \frac{\partial S(t,a)}{\partial a} = -(\mu(a) + \theta(a) + \lambda(t,a))S(t,a) \\ & \frac{\partial I(t,a)}{\partial t} + \frac{\partial I(t,a)}{\partial a} = \lambda(t,a)S(t,a) - (\mu(a) + \gamma(a)) I(t,a) \\ &\frac{\partial R(t,a)}{\partial t} + \frac{\partial R(t,a)}{\partial a} = \theta(a)S(t,a) + \gamma(a)I(t,a) - \mu(a)R(t,a) \\& S(t,0) = \int_0^\omega m(a) (S(t,a) + (1-q)I(t,a) + R(t,a))\, da \\& I(t,0) = q\int_0^\omega m(a)I(t,a)\, da \\ &R(t,0) = 0 \\ &S(0,a) = S_0(a),\ I(0,a) = I_0(a),\ R(0,a) = R_0(a) \end{align} \right.$$ここで、\(S(t,a),I(t,a),R(t,a)\)はそれぞれ感受性人口、感染人口、隔離人口の年齢密度関数とし、\(m(a)\)は年齢別出生率、\(\mu(a)\)は自然死亡率、\(\gamma(a)\)は隔離率、\(q\)は感染人口から産まれた新生児が感染している(これを垂直感染と呼びます)確率、\(\theta(a)\)はワクチン接種による免疫化率を表します。また、\(\lambda(t,a)\)は\(a\)歳の感受性人口に作用する感染力(force of infection)であり、以下のように与えられます:$$\lambda(t,a) = \frac{1}{N(t)}\int_0^\omega \beta(a,\sigma) I(t,\sigma)\, d\sigma .$$ただし、\(P(t,a) := S(t,a) + I(t,a) + R(t,a)\)とするとき、$$N(t) := \int_0^\omega P(t,a), da$$です。つまり\(N(t)\)は時刻\(t\)における総人口です。\(\beta(a,\sigma)\)は\(a\)歳の感受性個体と\(\sigma\)歳の感染個体の間の感染率であり、\(\omega\)は最大年齢、\(S_0,I_0,R_0\)は初期データです。
このシステムにおいて、全人口の年齢分布\(P(t,a)\)は以下の安定人口モデルにしたがうことがわかります:$$\left\{ \begin{align} &\frac{\partial P(t,a)}{\partial t} + \frac{\partial P(t,a)}{\partial a} = -\mu(a) P(t,a) \\& P(t,0) = \int_0^\omega m(a) P(t,a)\, da\\& P(0,a) = P_0(a) \end{align} \right.$$ここで、\(P_0(a) := S_0(a) + I_0(a) + R_0(a)\)です。したがって、全人口のサイズと年齢分布の時間発展は、感染の状況とは無関係に決定されます。安定人口モデルの一般論により、人口の年齢分布は時間に依らない分布へ収束します:$$\lim_{t \to \infty} \frac{P(t,a)}{\int_0^\omega P(t,a)\, da} = c(a) := \frac{e^{-\lambda_0 a} m(a) \ell(a)}{\int_0^\omega e^{-\lambda_0 a} m(a) \ell(a)\, da} .$$ここで\(\lambda_0\)は全人口の自然成長率、\(\ell(a)\)は生残率\(\ell(a) = \exp (-\int_0^a \mu(\sigma)\, d\sigma)\)です。\(c(a)\)は安定年齢分布といいます。そこで以下では、ホスト人口は既に安定年齢分布に到達しているものと仮定します。
全体人口が定常状態にあると仮定します。また、簡単のため、垂直感染が存在せず、ワクチン接種も行わないと仮定します。すなわち、\(q=0, \theta(a) \equiv 0\)であり、$$\int_0^\omega m(a)\ell(a)\, da = 1$$が成り立つとします。また、全人口が定常状態であるとすれば、任意の\(t>0\)について、$$P(t,a) = P_0(a) := B\ell(a) .$$ここで\(B\)は粗出生率と呼ばれる定数です。このとき\(c(a) = b\ell(a),\ b=1/\int_0^\omega \ell(a)\, da\)となります。\(b\)は規格化された定常人口分布の粗出生率です。
漸近挙動と\(\mathcal{R}_0\)
このとき、全人口が感受性であるような定常状態に少数の感染者が発生した場合を考えます。このとき、線形化された感染者のダイナミクスを考えれば、
$$\left\{ \begin{align} &\frac{\partial I(t,a)}{\partial t} + \frac{\partial I(t,a)}{\partial a} = c_2(a) \int_0^\omega \beta(a, \sigma) I(t,\sigma)\, d\sigma - (\mu(a) + \gamma(a))I(t,a)\\ &
I(t,0) = q\int_0^\omega m(a) I(t,a)\, da \end{align} \right.$$を得ます。この線形方程式を、Laplace変換を用いて形式に解くことを考えます。IのLaplace変換を$$\hat I(a,\alpha) = \int_0^\infty e^{-\alpha t} I(t,a)\, da$$とし、システムの上の方程式の両辺について時間\(t\)に関するLaplace変換をとれば、以下を得ます:$$\frac{d\hat I(t,a)}{da} = I_0(a) + c(a) \int_0^\omega \beta(a,\sigma)\hat I(\sigma,\alpha)\, d\sigma - (\alpha + \mu(a) + \gamma(a))\hat I(a,\alpha) .$$これを\(\hat I(\cdot, \alpha)\)に関する微分方程式とみて解けば、$$\begin{align}&\hat I(a,\alpha) = q \langle m,\hat I(\cdot, \alpha) \rangle e^{-\alpha a}\ell(a) \Gamma(a) \\ &\quad + \int_0^a e^{-\alpha (a-\sigma)}\frac{\ell(a)\Gamma(a)}{\ell(\sigma)\Gamma(\sigma)} [I_0(\sigma) + c(\sigma) \langle \beta(\sigma,\cdot), \hat I(\cdot, \alpha) \rangle ]\, d\sigma \end{align}$$を得ます。ここで、\(\Gamma(a) :=\exp(-\int_0^a \gamma(\sigma)\, d\sigma)\)です。上の等式の両辺に\(m(a), \beta(a,\sigma)\)をそれぞれ掛けてから\([0,\omega]\)上で積分すれば、\( \langle m,\hat I(\cdot, \alpha) \rangle, \langle \beta(a,\cdot), \hat I(\cdot, \alpha) \rangle\)に関する2つの等式が得られ、それらを連立して解けば\(\hat I(a,\alpha)\)が求まります。とくに、\(I\)によらない作用素\(T(\alpha)\)と関数\(f(a,\alpha)\)があって、$$\langle \beta(a,\cdot),\hat I(\cdot,\alpha) \rangle = (1 - T(\alpha))^{-1}f(\cdot,\alpha)(a)$$となります。安定人口モデルとのアナロジーから、\(\hat I(a,\alpha)\)の特異点、すなわち作用素\( (1-T(\alpha) )\)が有界な逆をもたないような点の集合\( \Lambda := \{ \lambda \in \mathbb{C} \mid 1 \in \sigma ( T(\alpha) ) \} \)によって\(I(t,a)\)の漸近挙動が決まると予想されます。実際、\(T(\alpha)\)がコンパクトかつノンサポーティングという性質をもち、かつ垂直感染がないならば、\(T(\alpha)\)は全ての\(\alpha \in \mathbb{R}\)で正作用素であり、方程式$$r(T(\alpha))=1$$(\(r\)はスペクトル半径を表す)は唯一つの実根\(\lambda_0 \in \Lambda\)をもち、しかも\(\lambda_0\)は\(\Lambda\)の任意の他の元の実部よりも大きいことが示されます。そして、感染症の流行の拡大がおこるためには、\(r( T(0) ) > 1\)となることが必要十分であることが知られています。
いま、作用素\(T\)を$$(Tf)(a) := c(a) \int_0^\omega \beta(a,x) \int_0^x \frac{\ell(x)\Gamma(x)}{\ell(\sigma)\Gamma(\sigma)}f(\sigma)\, d\sigma dx$$と定義すると、\(Tf\)はちょうど感染したばかりの感染者の年齢分布\(f\)から水平感染によって生産された2次感染者の年齢分布を表します。その意味で\(T\)は次世代作用素(next generation operator, NGO)と呼ばれます。次世代作用素のスペクトル半径は$$r(T) = \lim_{n \to \infty} \sqrt[n]{\| T^n \|}$$であり、漸近的には連続する世代の感染人口のサイズの比と解釈されるので、システムの基本再生産数\(\mathcal{R}_0\)であると考えられます。このとき、掛け算作用素\(L:f\mapsto cf\)について\(T = L^{-1}T(0)L\)が成り立つことから、$$r(T) = r( T(0) )$$を得ます。以上のことから、次の定理が得られます:
非自明平衡解と\(\mathcal{R}_0\)
次に、定常人口下におけるシステムの定常解の存在を議論してみましょう。\( (S^\ast, I^\ast) \)を定常解として、定常解に対応する感染力を\(\lambda^\ast\)とすれば、$$\left\{ \begin{align} &\frac{dS^\ast(a)}{da} = -(\mu(a) + \theta(a) + \lambda^\ast(a))S^\ast(a) \\ &\frac{dI^\ast(a)}{da} = \lambda^\ast(a)S^\ast(a) - (\mu(a) + \gamma(a))I^\ast(a) \\ &S^\ast(0) = B - I^\ast(0) \\ &I^\ast(0) = q \int_0^\omega m(a) I^\ast(a)\, da \\ &\lambda^\ast(a) = \left( B\int_0^\omega \ell(a)\, da \right)^{-1} \int_0^\omega \beta(a,\sigma)I^\ast(\sigma)\, d\sigma . \end{align}\right.$$定数変化法を用いれば、以下の表現が得られます:$$\begin{align} &S^\ast(a) = S^\ast(0)\ell(a)\Theta(a)e^{-\int_0^a \lambda^\ast(\sigma)\, d\sigma},\\ &I^\ast(a) = I^\ast(0) \ell(a)\Gamma(a) + \int_0^a \frac{\Gamma(a)\ell(a)}{\Gamma(\sigma)\ell(\sigma)} \lambda^\ast(\sigma) S^\ast(\sigma)\, d\sigma \end{align}.$$ ここで、$$\Theta(a) := \exp \left( -\int_0^a \theta(\sigma)\, d\sigma \right)$$です。
いま、\(\lambda^\ast\)の非線形汎関数\(F(\lambda^\ast)\)を、$$\begin{align} F(\lambda^\ast) &:= \frac{qb^{-1}}{1 - q\int_0^\omega m(a)\ell(a)\Gamma(a)\,da} \\ &\times \int_0^\omega m(a) \int_0^a \frac{\Gamma(a)\ell(a)}{\Gamma(\sigma)\ell(\sigma)}\Theta(\sigma)c(\sigma)\lambda^\ast(\sigma) e^{-\int_0^\sigma \lambda^\ast(z)\, dz}\, d\sigma \, da \end{align}$$と定義すれば、\(I^\ast(0)\)は\(F(\lambda^\ast)\)を用いて次のように表されます:$$I^\ast(0) = \frac{BF(\lambda^\ast)}{1 + F(\lambda^\ast)}$$また、\(S^\ast(a), I^\ast(a)\)の表現を\(\lambda^\ast\)の式に代入すると\(S^\ast, I^\ast\)に関する文字を消去することができて、\(\lambda^\ast\)についての関係式$$\begin{align} \lambda^\ast(a) &= \frac{bF(\lambda^\ast)}{1 + F(\lambda^\ast)} \int_0^\omega \beta(a,\sigma)\ell(\sigma)\Gamma(\sigma)\, d\sigma \\ &+ \frac{1}{1 + F(\lambda^\ast)} \int_0^\omega \beta(a,\sigma) \int_0^\sigma \frac{\Gamma(\sigma)\ell(\sigma)}{\Gamma(z)\ell(z)} c(z)\Theta(z)\lambda^\ast(z) e^{-\int_0^z \lambda^\ast(s)\, ds}\, dz\, d\sigma \end{align}$$を得ます。この等式の右辺を\(\Phi(\lambda^\ast)\)と書きましょう。すると、非線形作用素\(\Phi\)が正の不動点を持てば、それが非自明な定常解に対応する感染力であり、感染力が求まればエンデミックな定常解が計算されることになります。一般に\(\Phi\)の不動点の個数は簡単には分かりませんが、ワクチンによる免疫化が無い場合を考えれば、侵入条件が満たされればエンデミックな定常解が存在します。
このことは、次のような流れで証明されます: \(\Phi(x)\)の\(x=0\)におけるFréchet微分を\(\Phi^\prime[0]\)とすれば、\(u \in L^1\)に対して、$$\begin{align} (\Phi^\prime[0]u)(a) &= \int_0^\omega \beta(a,\sigma) \int_0^\sigma \frac{\Gamma(\sigma)\ell(\sigma)}{\Gamma(z)\ell(z)} c(z)\Theta(z) u(z) \, dz\, d\sigma \\ &+ \frac{q\int_0^\omega \beta(a,\sigma)\ell(\sigma)\Gamma(\sigma)\, d\sigma}{1 - q\int_0^\omega m(a)\ell(a)\Gamma(a)\, da}\int_0^\omega m(a) \int_0^a \frac{\Gamma(a)\ell(a)}{\Gamma(\sigma)\ell(\sigma)} c(\sigma)\Theta(\sigma) u(\sigma)\, d\sigma\, da . \end{align}$$正値錐\(L_+^1\)の\(\Phi\)による像は正有界であるので、Krasnoselskiiの不動点定理により、\(\Phi^\prime[0]\)がコンパクト、ノンサポーティングかつ\(r(\Phi^\prime[0]) > 1\)ならば、\(\Phi\)は少なくとも1つの正の不動点をもちます。とくに\(\theta \equiv 0\)であれば、\(\Phi^\prime[0] = T(0)\)なので、\(T(0)\)がコンパクト、ノンサポーティングかつ侵入条件\(r(T) = r(T(0)) > 1\)が満たされていれば、少なくとも1つのエンデミックな定常解が存在します。
このときさらに垂直感染がないと仮定すると、\(q = F = 0\)として、\(\Phi(x) \leq \Phi^\prime[0]x\)となります。この場合、\(r(T)\leq 1\)であれば自明な定常解しかないことが知られていますが、この自明な定常解は大域安定となります。実際、\(S(t,a)\)が与えられていると考えると、\(I(t,a)\)は非自律的線形問題$$\left\{\begin{align} &\frac{\partial I(t,a)}{\partial t} + \frac{\partial I(t,a)}{\partial a} = \frac{S(t,a)}{B\int_0^\omega \ell(a)\, da}\int_0^\omega \beta(a,\sigma)I(t,\sigma)\, d\sigma - (\mu(a) + \gamma(a))I(t,a) \\ &I(t,0) = q\int_0^\omega m(a) I(t,a)\, da \end{align}\right.$$の正値解になっています。ここで、$$S(t,a) = \left\{ \begin{align}& B\ell(a)\exp\left( -\int_0^t \lambda(t-a+\sigma,\sigma)\, d\sigma \right) \quad (t - a > 0)\\ & S_0(a-t) \frac{\ell(a)}{\ell(a-t)} \exp \left( -\int_0^a \lambda(\sigma, a-t+\sigma)\, d\sigma \right) \quad (a - t > 0) \end{align}\right.$$であるから、\(t\)をじゅうぶん大きくとれば、$$\frac{S(t,a)}{B\int_0^\omega \ell(a)\, da} \leq \frac{\ell(a)}{\int_0^\omega \ell(a)\, da} = c(a)$$となります。ゆえに比較定理により、\(J(t,a)\)を初期条件\(J(0,a) = I_0(a)\)のもとでの線形化問題(前セクションで登場した線形化方程式のことです)の解とすれば、$$I(t,a) \leq J(t,a)$$となります。定理3により、\(q=0\)の場合、\(r(T)<1\)のときは$$\lim_{t \to \infty}J(t,a) = 0$$なので、\(I=0\)は大域的に漸近安定になっているとわかります。
したがって、次のような定理が得られます:
参考文献
[1] 稲葉寿 (2002), "数理人口学", 東京大学出版会