Rei Frontier Tech Blog

人工知能を活用した位置情報分析プラットフォーム「SilentLog Analytics」を運営する、レイ・フロンティア株式会社のエンジニアメンバーで運営する技術ブログです。

感染症モデルとR0:その2

レイ・フロンティア株式会社のデータアナリストの齋藤です。前回の記事の続きを書いていこうと思います。

年齢構造をもつ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) )$$を得ます。以上のことから、次の定理が得られます:

定理3. (閾値原理)水平感染において\(\mathcal{R_0} = r(T)\)とするとき、
\(\mathcal{R}_0>1\)ならば伝染病は侵入できる。
\(\mathcal{R}_0<1\)ならば伝染病は侵入できない。

非自明平衡解と\(\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\)は大域的に漸近安定になっているとわかります。
したがって、次のような定理が得られます:

定理4.\(q = \theta = 0\)と仮定する。また、次世代作用素\(T\)はコンパクトかつノンサポーティングであると仮定する。
\(\mathcal{R}_0<1\)ならば定常解は自明なものに限り、大域的に漸近安定である。
\(\mathcal{R}_0>1\)ならば自明な定常解は不安定であり、少なくとも1つのエンデミックな定常解が存在する。
こうして、年齢構造を導入した場合においても、基本再生産数\(\mathcal{R}_0\)は非自明定常解の存在と自明定常解の安定性の両方をきめる閾値となっていることがわかりました。

参考文献

[1] 稲葉寿 (2002), "数理人口学", 東京大学出版会

感染症モデルとR0:その1

レイ・フロンティア株式会社のデータアナリストの齋藤です。すっかり寒くなってきましたね。

前々回の記事で、人口動態を記述する数理モデルにおいて、基本再生産数\(\mathcal{R}_0\)という量が重要な役割をもつことを述べました。次は、今回と次回の記事の2回にわたって、感染症の数理モデルについても\(\mathcal{R}_0\)という概念が有用であることを伝えたいと思います。

おさらい:安定人口モデルにおける\(\mathcal{R}_0\)

\(p(t,a)\)を時刻\(t\)における年齢\(a\)の年齢密度関数としたとき、\(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), \quad t>0,\ a>0 \\ & p(0,t) = \int_0^\omega \beta(a)p(t,a)\, da,\quad t>0 \\ & p(0,a) = p_0(a),\quad a \in [0,\omega] \end{align} \right.$$ここで\(\mu(a),\beta(a)\)はそれぞれ年齢別死亡率、年齢別出生率です。このモデルを変形することにより、再生方程式$$B(t) = G(t) + \int_0^t\Psi(a)B(t-a)\, da$$を得ます。ここで、$$\begin{align}&G(t) = \int_t^\omega \beta(a) \frac{\ell(a)}{\ell(a-t)}p_0(a-t)\, da \\ &\Psi(a) = \beta(a)\ell(a)\end{align}$$です。\(G(t)\)は初期時刻において生存している人口が時刻\(t\)で子を産むことに、\(\Psi(a)\)は\(a\)年前に産まれた人が\(a\)歳まで生存して子を産むことに、そして\(B(t)\)は時刻\(t\)で産まれた子の人数にそれぞれ対応しています。そこで、$$\mathcal{R}_0 = \int_0^\infty \Psi(a)\, da$$と定義すると、\(\mathcal{R}_0\)は一人の人が生涯に産む子の平均人数を意味します。
\(\mathcal{R}_0\)を基本再生産数と呼びます。

安定人口モデルにおいて、以下の定理が成り立ちます:

定理1. (Sharpe-Lotka-Fellerの定理)安定人口モデルにおける再生方程式に関して、以下が成り立つ:$$B(t) = q_0 e^{\lambda_0 t} (1 + O(e^{-\eta t})).$$ここで\(\eta\)はある正数、\(\lambda_0\)はLotkaの特性方程式$$\int_0^\omega e^{-\lambda a}\Psi(a)\, da = 1$$の唯一つの実根であり、$$q_0 = \frac{\int_0^\infty e^{-\lambda_0 t}G(t)\, dt}{\int_0^\infty a e^{\lambda_0 a}\Psi(a)\, da}$$である。
\(\lambda_0\)は内的成長率と呼ばれる量であり、基本再生産数との間に以下の符号関係が成り立ちます:$$\mathrm{sgn}\mathcal({R}_0-1) = \mathrm{sgn}\lambda_0$$すなわち、\(\mathcal{R}_0>1\)ならば人口は指数関数的に増大し、\(\mathcal{R}_0<1\)ならば減少します。「一人あたり一人より多く産めば人口は増えるし、一人より少なく産めば人口は減る」という直感的な推論が、この人口モデルにおいては確かに成立しているわけです。

感染症モデルにおける\(\mathcal{R}_0\)

次に、感染症の数理モデルを考えてみましょう。感染症モデルにおいても、様々なモデルについてそれに対応する\(\mathcal{R}_0\)が計算されています。人口の場合とは違い、感染症モデルにおける\(\mathcal{R}_0\)とは病気に感染した一人の人間が生涯にわたって病気をうつす人数のことです。しかしこの場合も人口モデルと同様に、「一人あたり一人より多くうつせば感染人口は増えるし、一人より少なくうつせば感染人口は減る」という法則が成立すると期待できます。
本記事では、年齢を考慮しない常微分方程式モデルを考えます。感染者の年齢を考慮した感染症モデルについては次回に述べます。

Kermack-McKendrickモデル

ある感染症について、\(S(t)\)を時刻\(t\)における感受性人口(susceptibles, 感染する可能性のある人口)、\(I(t)\)を感染人口(infected/infectious, 感染していて、かつ感染させる能力のある人口)、\(R(t)\)を隔離された人口(recovered.removed, 病気からの回復による免疫保持者、あるは隔離者、死亡者)とし、全人口を3つに分別します。今回は、病気の流行期間が非常に短く、感染プロセスにおいて出生や死亡などの人口動態は無視できると考えます。このとき、KermackとMcKendrickが提起したモデルは以下のような常微分方程式系です:
$$\left\{ \begin{align} & \frac{dS(t)}{dt} = -\beta S(t) I(t) \\ & \frac{dI(t)}{dt} = \beta S(t) I(t) - \gamma I(t) \\ & \frac{dR(t)}{dt} = \gamma I(t) \end{align} \right.$$ここで、\(\beta\)は感染率、\(\gamma\)は隔離率を表す正の定数です。一般に、このように人口を感受性・感染・隔離の3種類に類別した伝染病モデルをSIRモデルといいます。このモデルにおいて、総人口\(S(t)+I(t)+R(t)\)は$$\frac{d(S(t)+I(t)+R(t))}{dt} = 0$$より一定です。これを\(N\)とおきます。上の微分方程式系の定常解のひとつとして、$$(S(t),I(t),R(t)) = (N,0,0)$$があります。これは感染者がまだ一人もいない状態に対応するものです。このような状態を、感染者のいない定常状態(disease-free steady state, DFSS)と呼びます。感染症モデルの解析においては、disease-freeな人口集団に少数の感染者が発生した場合に流行(= 感染人口の増大)が発生する条件を調べるのが定石です。この条件は侵入条件と呼ばれます。\((N,0,0)\)の状態の集団にごく少数の感染者が発生した状態を考えます。\(S(t) \approx N, I(t) \approx 0\)より、感染者数のダイナミクスは平衡点\((N,0,0)\)での線形化方程式$$\frac{dI(t)}{dt} = (\beta N - \gamma)I(t)$$によって記述されます。したがって、流行が発生した直後における感染者人口は$$I(t) \approx I(0) e^{(\beta N - \gamma )t}$$というMalthus法則にしたがっています。このことから、\(\beta N - \gamma > 0\)ならば感染人口は増大し、\(\beta N - \gamma < 0\)ならば再びゼロに近づいていきます。条件\(\beta N - \gamma > 0\)は$$\mathcal{R}_0 := \frac{\beta N}{\gamma} > 1$$と同値ですが、この\(\mathcal{R}_0\)が基本再生産数と呼ばれています。\(\beta N\)はサイズ\(N\)の感受性人口集団において、一人の感染者が単位時間あたりに生産する2次感染者数を表します。そして、\(1/\gamma\)は感染者の感染状態にある時間の平均を表します。したがって、\(\mathcal{R}_0\)は感受性人口に侵入した感染者が一人あたりに生産する2次感染者数と解釈できます。

修正Kermack-McKendrickモデル

前節では人口動態のないSIRモデルを考えましたが、次は出生率を\(b\)、自然死亡率を\(\mu\)として、ホスト人口に人口動態を導入してみます。修正されたSIRモデルは以下のとおりです:
$$\left\{ \begin{align} & \frac{dS(t)}{dt} = b - \mu S(t) - \beta S(t) I(t) \\ & \frac{dI(t)}{dt} = \beta S(t) I(t) - (\mu + \gamma) I(t) \\ & \frac{dR(t)}{dt} = -\mu R(t) + \gamma I(t) \end{align} \right.$$このとき、総人口\(N(t) = S(t) + I(t) + R(t)\)は安定な平衡値\(b/\mu\)をもつので、はじめから全人口は一定値\(N=b/\mu\)をもつと仮定しておきます。このとき、上の常微分方程式系は以下の2次元の系と同値です:$$\left\{\begin{align}&\frac{dS(t)}{dt} = \mu N - \mu S(t) - \beta S(t) I(t) \\& \frac{dI(t)}{dt} = \beta S(t) I(t) - (\mu + \gamma)I(t)\end{align}\right.$$前節と同様に、初期侵入の状況(\(S(t) \approx N, I(t) \approx 0\))を考えると、線形化方程式は$$\frac{dI(t)}{dt} = (\gamma + \mu)\left[ \frac{\beta N}{\gamma + \mu} - 1 \right] I(t)$$となり、基本再生産数は$$\mathcal{R}_0 = \frac{\beta N}{\gamma + \mu}$$と求まります。さらに平衡点を探せば、$$E_1 = \left( \frac{b}{\mu}, 0 \right), \ E_2 = \left( \frac{N}{\mathcal{R}_0}, \frac{\mu}{\beta}(\mathcal{R}_0 - 1) \right)$$という2つの定常解が存在することがわかります。\(E_1\)はつねに存在する自明な定常解であり、感染者のいない定常状態です。一方、\(E_2\)は閾値条件\(\mathcal{R}_0 > 1\)が満たされた場合にのみ正となり、生物学的に意味をもちます。この定常解は伝染病が人口に定着して共存している、エンデミックな定常状態(endemic steady state, ESS)といいます。このとき、次のような定常解の安定性に関する閾値定理が成り立ちます:

定理2.\(\mathcal{R}_0 > 1\)ならば、DFSSは大域的に漸近安定である。
\(\mathcal{R}_0 < 1\)ならば、DFSSは不安定であり、ESSは大域的に漸近安定である。

次回に続きます。

テンソルの定義について: 普遍性と多重配列

レイ・フロンティア株式会社のデータアナリストの齋藤です。最近、食事がうどんばかりです。

本記事では、機械学習において盛んに応用され、深層学習のライブラリ tensorflow の名前の由来にもなっているテンソルについて述べます。機械学習においては、たとえば映画館の顧客について、「顧客の好み」「映画の種類」「放映された季節」など複数の"次元"を持つデータを多重配列として表現し、その類似度を求めるということが行われます。実際の計算処理のなかではテンソルは単なる多重配列として表現されますし、そのような理解で十分とする解説記事も多くあります。その一方で、テンソルの定義が気になって調べてみると難しい数学の説明ばかりでわけが分からなかった……という経験をした方も多いと思われます。そこで、本記事ではテンソルの数学的な定義と多重配列との関連について、その"気持ち"の部分をわかりやすく解説することを目指しました。少しでも納得に近づければ幸いです。

テンソルを導入するモチベーション

テンソル(tensor)という言葉の由来は、弾性学における"引っ張る力"を意味するテンション(tension)といわれています。現在テンションの概念は、材料力学において、物体の内部に加えられる力の大きさと方向を表し、物体の変形・破壊などを調べるのに使われる応力テンソル(stress tensor)に継承されています。まずは、この応力テンソルによって、テンソルの概念が自然に現れる様子を見ていきます。

応力テンソル

物体(連続体)の点\(r\)で、単位法線ベクトル\({\bf n}\)と面積\(\delta S\)をもつ面を考えます。この面を通じて、面の表側(\(r\)が向いている側)から裏側へもたらされる力を\(T({\bf n})\delta S\)と書くことにします。\(T({\bf n})\)は点\(r\)においてその面に作用する、単位面積あたりの力を表します。これを応力といいます。
ここで、物体を微小な四面体の集まりとして近似することを考え、面と接している四面体にはたらく力を考えます。ここで、四面体は3つの座標軸に垂直な3つの面をもっているとします。

f:id:reifrontier:20171011111719p:plain

この四面体にはたらく力は慣性力・外力・面積力の3種類ありますが、そのうち慣性力と外力は体積力であり、四面体の長さのスケール\(\delta l\)の3乗のオーダーを持ちます。面積力は2乗のオーダーなので、これらの力は四面体が十分に小さければ無視しても構いません。このとき、四面体に作用する面積力がつり合うと考えて、つり合いの式$$T({\bf n})\delta S + \sum_{j=1}^3 T(-{\bf e_j})\delta S_j = 0$$を得ます。ここで\(\delta S, \delta S_j\)はそれぞれの面の微小面積であり、\({\bf e}_j\)は\(x_j\)軸方向の単位ベクトルをさします。
\(\delta S_j = ({\bf n} \cdot {\bf e}_j ) \delta S\)であることから、$$T({\bf n}) = \sum_{j = 1}^3 T({\bf e}_j)({\bf n} \cdot {\bf e}_j )$$ となるので、3つの座標軸に垂直な面に作用する応力\(T({\bf e}_j)\ (j = 1,2,3)\)を知ることができれば、任意の面に作用する応力\(T({\bf e})\)を知ることができます。
ベクトル\(T\)の第\(i\)成分を\(T_i\)と書いたうえで、$$\sigma _{ij} := T_i({\bf e}_j)$$と定義すれば、\({\bf n} \cdot {\bf e}_j = {\bf n}_j\)より、$$T_i({\bf n}) = \sigma_{ij} {\bf n}_j$$と書けます。この9個の量\(\sigma_{ij}\)を成分とする量を応力テンソルと呼びます。

応力テンソルは9次元ベクトル空間のベクトルですが、力や位置のように"まっすぐ縦に数を並べた"量ではなく、\(i\)と\(j\)という2つのパラメータによって行列のように"2次元的に並べた"量として表現されます。その意味で、応力テンソルは通常のベクトルとは区別されるべき存在といえます。このような、成分を2次元的に並べられる量は2階のテンソルと呼ばれます。同様に、3,4,...次元的に並べられる量は3,4,...階のテンソルと言うことができます。

応力テンソルなど、物理学に現れるテンソルは、上にみた例のように多次元的に並べた時の各成分\(\sigma_{ij}\)のことをさして言うのが通例です。一方、近年応用の目覚ましい機械学習におけるテンソルは多次元的に並べられたデータの集まりであり、数学におけるテンソル積は普遍性と呼ばれるあるベクトル空間の性質から抽象的に定義される量です。このように学問領域によって言葉の使い方にばらつきがあることが、テンソルの定義の理解の妨げになっているのではないかと邪推しています。

テンソル積の定義

ここでは、数学的なテンソル積の定義を紹介することになります。しかし困ったことに、数学におけるテンソル積の定義にも異なる流儀が存在するのです(もちろん、どの定義を用いても結局は同じものが作られます)。本記事では、ベクトル空間のテンソル積の3種類の定義を紹介します。定理の証明を詳しく追いたい方は、参考文献に挙げた線形代数の本を参照してください。とくに[5]は本記事の執筆におおいに参考にしています。

方法1: 普遍性による定義: 基底を用いる

\(V,W\)を線形空間とします。係数体は何でもいいですが、\(\mathbb{R}\)や\(\mathbb{C}\)だと思って差し支えません。直積\(V\times W\)上の双線形写像の集合\(\mathscr{L}(V,W;U)\)を考えます。ここで、像の空間\(U\)を色々変えてみることを考えます。そうすると、\(U\)によらない特別な線形空間\(U_0\)と特別な双線形写像\(\iota : V \times W \to U_0\)が存在することが知られています。どういう風に特別なのかというと、\(U\)と、双線形写像\(\Phi:V \times W \to U\)をどのようにとったとしても、\(\iota\)によって\(\Phi\)は"線形化"されてしまうのです! この"双線形写像の線形化"という側面は、線形空間におけるテンソル積では本質的です。
そのような\(\iota\)を実際に作ってみましょう。\(V,W\)をそれぞれ\(n,m\)次元線形空間とし、それらの基底を1つずつ選んでそれぞれ\(\langle e_1,e_2,\dots,e_n\rangle , \langle f_1,f_2,\dots,f_m\rangle\)とおきます。また、\(U_0\)を\(mn\)次元線形空間として、\(\langle g_{11},g_{12},\dots, g_{1m},g_{21},g_{22},\dots, g_{2m},\dots, g_{n1},g_{n2},\dots, g_{nm}\rangle\)をその基底とします(\(n\)行\(m\)列に並べています)。そして、写像\(\iota: V\times W \to U_0\)を$$\begin{align}&V \ni v = \sum_i \alpha_i e_i,\ W \ni w = \sum_j \beta_j f_j\ (\alpha_i, \beta_j \in K)\ に対して\\& \iota(v,w) = \sum_{i,j} \alpha_i \beta_j g_{ij} \end{align}$$ によって定義します。
このとき、\((U_0,\iota)\)に関して次が成り立ちます:

定理1.
(1)\((U_0,\iota)\)に関して以下が成立する:
\((\otimes)\) \(^\forall \Phi \in \mathscr{L}(V,W;U)\)に対して、\(\Phi=F\circ \iota\)となる線形写像\(F: U_0 \to U\)がただ一つ存在する。

(2) (1)の性質をもつ\((U_0,\iota)\)は次の意味で一意的である: 線形空間と双線形写像の組\((U_0,\iota), (U_0^\prime,\iota^\prime)\)がともに(1)の性質をもつならば、線形同相写像\(F_0: U_0\to U_0^\prime\)であって\(F_0\circ \iota = \iota^\prime\)となるものがただ一つ存在する。

写像\(\iota\)の定義は基底を1つ選ぶことによってなされていますが、(2)によって基底によらない定義を与えていることがわかります。
重要なのは性質\((\otimes)\)です。感覚的に述べると、任意の(任意の空間への!)双線形写像\(V\times W \to U\)は、写像\(\iota\)を通すことによって線形写像\(U_0 \to U\)へと"線形化"することができるのです。このような、\(V\)と\(W\)のみによって定まり、ある意味で任意の\(U\)と\(\mathscr{L}(V,W;U)\)を統制する性質のことは、普遍性(universality)とよばれています。"線形化"写像\(\iota\)を改めて\(\otimes\)、\(U_0\)を\(V\otimes W\)と書き、\(V\otimes W\)をベクトル空間\(V,W\)のテンソル積、\(\otimes\)をテンソル積の標準写像と呼びます。また、ベクトル\(v\in V,w \in W\)について\(\otimes(v,w)\)を\(v\otimes w\)と表記します。これで、ひとまず我々が求めていたテンソル積の定義が手に入りました。

テンソル積は以下の性質をもちます:

命題2.(\(\otimes\)の双線形性)
\(\alpha, \beta \in K,\ v,v_1,v_2 \in V,\ w,w_1,w_2 \in W\)として、$$\begin{align}(\alpha v_1 + \beta v_2)\otimes w = \alpha (v_1 \otimes w) + \beta (v_2 \otimes w)\\ v\otimes (\alpha w_1 + \beta w_2) = \alpha (v \otimes w_1) + \beta (v \otimes w_2)\end{align}$$
また、\(V,W\)の基底\(\langle e_1,e_2,\dots,e_n\rangle , \langle f_1,f_2,\dots,f_m\rangle\)に対して、\(V\otimes W\)の元の集まり$$\begin{align}&\langle e_1\otimes f_1,e_1\otimes f_2,\dots, e_1\otimes f_m,\\ &e_2\otimes f_1,e_2\otimes f_2,\dots, e_2\otimes f_m,\\& \dots \dots \dots \dots, \\ &e_n\otimes f_1,e_n\otimes f_2,\dots, e_n\otimes f_m\rangle \end{align}$$は\(V\otimes W\)の基底になっています。したがって、$$\dim (V\otimes W) = mn = \dim V \dim W$$です。

方法2: 普遍性による定義: 基底を用いない

方法1では具体的に\(V,W,U\)の基底を1つとることによって標準写像を定義しましたが、線形空間の双対空間を用いれば、基底をとることなしにテンソル積を構成することもできます。
\(V,W\)の双対空間をそれぞれ\(V^\ast, W^\ast\)とし、$$U_0 = \mathscr{L}(V^\ast, W^\ast;K)$$とおきます。与えられた\((u,w) \in V\times W\)に対して、対応$$V^\ast \times W^\ast \ni (\phi, \psi) \mapsto \phi(v)\psi(w) \in K$$は\(V^\ast \times W^\ast\)から\(K\)への双線形写像を成すのでこれを\(\iota(u,w)\)とすればOKです。あとは方法1と同様にテンソル積を定義して、定理1や命題2を確かめることができます。

方法3: 商空間を用いる定義

3つめの方法は、線型空間の商空間を用いるものです。この流儀の定義は一見すると分かりにくいですが、重要なのはあくまで普遍性と双線形性であって、構成そのものを気にする必要はありません。
直積集合\(V\times W\)について、以下のような線形空間が定義できます:$$\mathscr{V}(V\times W) = \left\{ \sum_{(v,w) \in V\times W} \alpha(v,w) e(v,w) \middle| \alpha(v,w)\not=0となる(v,w)は高々有限個 \right\}$$ここで、\(\alpha(v,w) \in K\)であり、\(e(v,w)\)は\( (v,w)\in V\times W\)に対応する"もの"です。多項式\(\sum_i \alpha_i X^i\)の不定元のようなものだと思っておけばよいでしょう。
次に、\(\mathscr{V}(V\times W)\)の部分空間であって、次の形の元全体によって生成されるものを\(X\)とおきます:$$\left\{\begin{align} &e(v_1+v_2,w) - e(v_1,w) - e(v_2,w)\\ &e(v,w_1+w_2) - e(v,w_1) - e(v,w_2)\\ &e(\alpha v,w) - \alpha e(v,w)\\ &e(v,\alpha w) - \alpha e(v,w)\quad (v_1,v_2,v \in V,\ w_1,w_2,w \in W,\ \alpha \in K)\end{align}\right.$$商線形空間\(\mathscr{V}(V\times W)/X\)を\(U_1\)とおきます。また、\(\mathscr{V}(V\times W)\)から\(U_1\)への自然な射影を\(\pi\)とします。\(X\)の定義から、$$\pi (e(v_1+v_2,w)) - \pi (e(v_1,w)) - \pi (e(v_2,w)) = 0$$などが成り立ちます。ここから双線形性の雰囲気が何となく感じ取れるのではないでしょうか。
そして、写像\(e: V\times W \to \mathscr{V} (V\times W)\)と\(\pi\)の合成写像を\(\iota_1\)とおくと、次が成り立ちます:

定理3.
\(U_1,\iota_1\)は(方法1,2で定義した)テンソル積\(V\otimes W\)に同型である。
つまり、\(U_1,\iota_1\)はテンソル積のひとつの定義を与えているとみなすことができます。

高階テンソル

上で定義したテンソル積は2つの線形空間のある種の積でした。よってこれを2階のテンソルと呼びます。
テンソル積と線形空間のテンソル積を考えることも可能です。3つの線形空間\(V_1,V_2,V_3\)について、$$(V_1\otimes V_2)\otimes V_3 \cong V_1\otimes (V_2\otimes V_3)$$が成り立ちます。よって括弧の順番を気にすることなく\(V_1\otimes V_2\otimes V_3\)と書いてOKです。これを3つの線形空間のテンソル積、すなわち3階のテンソル積と呼びます。同様に、\(n\)個の線形空間\(V_1,V_2,\dots, V_n\)のテンソル積\(V_1\otimes V_2\otimes \cdots \otimes V_n\)を\(n\)階のテンソルと呼びます。\(n\)階のテンソルに関しても2階の場合と同様、\(n\)重の多重線形写像を"線形化"するという意味での普遍性と、\(n\)重線形性をもつことが確認できます。
スカラーは0階のテンソル、普通のベクトルは1階のテンソルと呼ばれます。データを"\(n\)次元に並べられる"ときに\(n\)階テンソル、ということになります。

多重配列としてのテンソル

さて、もう一つの目的であったテンソル積と多重配列との関係ですが、いったん定義を与えてしまえば簡単です。線形空間のテンソル積を方法1によって与えることにして、\(n\)個の線形空間\(V_1,V_2,\dots, V_n\)の基底を\(\langle e_{i,1},e_{i,2},\dots e_{i,m_i} \rangle,\ i=1,2,\dots, n\)とおきます。\(n\)階のテンソル\(V_1\otimes V_2\otimes \cdots \otimes V_n\)の元を、基底を用いて$$\sum_{j_1,j_2,\dots, j_n}\alpha(j_1,j_2,\dots, j_n) e_{1,j_i}\otimes e_{2,j_2}\otimes \cdots \otimes e_{n,j_n},\quad \alpha(j_1,j_2,\dots, j_n) \in K$$と書いたとき、係数\(\alpha(j_1,j_2,\dots, j_n)\)を\(n\)次元配列の成分\(\alpha[j_1][j_2]\dots [j_n]\)に対応させればOKです。プログラミングの世界における多重配列としてのテンソルの背景には、このような"普遍性(双線形写像の線形化)・双線形性"という興味深い構造が隠されていたのでした。

参考文献

[1] 石黒勝彦, 林浩平 (2016), "関係データ学習", 講談社

[2] 河合佑太, "気象力学を学ぶ上での基礎知識", http://epa.scitec.kobe-u.ac.jp/~ykawai/project/meteo_basic/

[3] 斎藤毅 (2007), "線形代数の世界", 東京大学出版会

[4] 佐武一郎 (1958), "線型代数学", 裳華房

[5] 横沼健雄 (1977), "テンソル空間と外積代数", 岩波書店

安定人口モデルについて

レイ・フロンティア株式会社でデータアナリストをしている齋藤です。これから、考えていることや興味を持って勉強していること、業務のなかで学んだ技術などについて定期的に発信していくことになりました。私個人としては数学や応用数理に興味を持って勉強を続けているので、そのような記事が中心になるかと思います。
お付き合い頂けると幸いです。

さて、現実の現象を理解するための手段として、現象を数理モデルによって記述し、モデルの性質を数学の問題として調べ、得られた結果を現実の言葉に翻訳するという数理的アプローチは基本的なものであり、今日ではあらゆる分野において数理モデルを用いた研究が行われています。今回はその中から、人口現象を説明するモデルである安定人口モデルを紹介します。

モデルの導出

歴史的経緯

人口学的現象を数理モデルによって説明しようとする試み自体は古くから存在しています。最も古いものは、おそらく13世紀にフィボナッチ(Leonardo Fibonacci)が彼の著書 "Liber Abaci"の中で紹介したフィボナッチ数列でしょう。
いま、ウサギのひとつがいを単位として、最初に1つがいのウサギが生まれてから\( n \)ヶ月後のつがいの数を\(P_n\)とします。また、

  • ウサギは生まれてから2ヶ月後に1つがいのウサギを産む。
  • ウサギが死ぬことは考えない。

とします。すると、漸化式 $$\left\{\begin{align} & P_{n+2} = P_{n+1} + P_n \\ & P_0 = 1 \end{align}\right.$$ が成り立ちます。\(n \to \infty\)とすると、比\(P_{n+1} / P_n\)は黄金比\(\phi = 1.6180339887...\)に収束します。これは、ウサギのつがい数は十分に時間が経つと公比\(\phi\)の等比数列に近づいていくことを意味します。このことを、つがい数は漸近的にMalthus径数\(\phi\)の指数的成長を行うと言い表します。後にも触れますが、人口が漸近的に指数的成長をすると仮定した場合の成長率を調べるというのは、人口モデルの解析における常套手段です。

人間の個体数の議論に数理モデルが広く受容されるようになったのは、マルサス(T. R. Malthus)が著書『人口論』において差分方程式モデル $$ P_{n+1} = (1 + r) P_n $$ を用いたことがきっかけでした。マルサスは、人口は指数関数的に増大するが食料は一次関数的にしか増大しないことから、人類は将来深刻な生活資源不足に陥ると警告しました。
人口モデルに年齢という概念がはじめて導入されたのは、ロトカ(A. J. Lotka)とシャープ(F. R. Sharpe)の論文([4])で提起された再生方程式(renewal equation) $$ B(t) = \int_0^\infty B(t-a) \ell(a) \beta(a)\, da$$ でした。ここで、\(B(t)\)は時刻\(t\)における単位時間あたりの出生数、\(\ell(a)\)は\(a\)歳までの生残率、\(\beta(a)\)は\(a\)歳の人口の年齢別出生率です。ロトカらは、再生方程式の解が漸近的に指数的成長をすることと、その漸近的な成長率が\(\ell\)と\(\beta\)によって決定されることを(数学的な厳密さは欠いていたものの)示しました。すなわち、外部との人口の流入出がない封鎖人口において、年齢ごとの死亡率と出生率が時間によらず一定ならば、時間が経つにつれて人口構造(ピラミッド型・壺型の図でおなじみの、年齢ごとの人口の比率のことです)は一定となり、総人口は指数的に変化するようになることを発見したのです。この理論は安定人口理論(stable population theory)と呼ばれ、人口学における基本的理論のひとつとなっています。彼らの結果はフェラー(W. Feller)によって数学的に正当化されました([1])が、その証明はLaplace変換を利用した極めて技巧的なものであり、それ以上の発展は到底見込めませんでした。それから40年以上経って、ウェブ(G. F. Webb)によってより一般の非線形人口モデルについて、それに対応する非線形半群による解を構成し、その振る舞いを調べるという関数解析的手法が打ち立てられました([5])。その後、メッツ(J. A. J. Metz)やディークマン(O. Diekmann)らによって年齢だけでなく位置やサイズなど様々な構造をもつモデルが研究され(その成果は論文集 [3] にまとめられています)、構造化人口動態学(structured population dynamics)という応用数学の一分野が確立されるに至ります。本記事で紹介するのは、構造化人口動態学の枠組みで考えられる最も単純なモデルです。あくまで紹介記事ですので、定理の証明などを細かく追うことはしません。きちんとした数学的な議論を見たい方は、たとえば[2], [3], [5]を参照するとよいでしょう。

パラメータの定義

まず、モデルに登場するパラメータたちについて説明します。
年齢構造をもつ人口集団において、時刻\(t\)における年齢密度関数(age-density function)を\(p(t,a)\)と書きます。積分$$ \int_{a_1}^{a_2} p(t,a)\, da $$は時刻\(t\)における年齢\(a_1\)以上\(a_2\)以下の人口になります。本記事で紹介するモデルは、\(p(t,a)\)に関する初期値境界値問題として定式化されます。
関数\(\mu (a)\)を、年齢\(a\)における個体の死亡の瞬間的な発生率とし、死亡力(force of mortality)あるいは年齢別死亡率(age-specific death rate)と呼びます。積分 $$\int_{a_1}^{a_2} \mu(a) p(t,a) \, da$$は単位時間あたりの年齢\(a_1\)以上\(a_2\)以下の死亡数となります。確率変数\(X\)を個体の死亡年齢とし、\(X\)の分布関数を\(F(a)\), 密度関数を\(f(a)\)とすると、\(\mu (a)\)は $$ \begin{align}\mu (a) &= \lim_{h \to 0} \frac{\mathbb{P}(a \leq X < a + h \mid a \leq X)}{h} = \lim_{h \to 0} \frac{\mathbb{P}(a \leq X < a + h)}{\mathbb{P}(a \leq X) h} \\ &= \frac{f(a)}{1 - F(a)} = -\frac{d}{da} (1 - F(a)) \end{align}$$ です(ここで、\(\mathbb{P}(A)\)は事象\(A\)がおこる確率を表します)。 このとき、 $$ \ell (a) := 1 - F(a) $$は個体が\(a\)歳まで生存する確率を意味し、生残率(survival rate)と呼ばれます。\(\mu(a)\)の定義から、微分方程式 $$ \frac{d\ell (a)}{da} = -\mu(a)\ell(a) $$が成り立ちます。これと\(\ell(0)\)(0歳の時点では確率1で生き残っている)から、$$ \ell (a) = \exp \left( - \int_0^a \mu (b) \, db \right) $$を得ます。いま、各個体は有限の最大年齢\(\omega\)をもつと仮定しましょう。すると\(\ell(\omega) = 0\)でなければならないので、$$\int_0^\omega \mu(a)\, da = \infty $$が必要です。たいていの場合\(\mu\)は上の式を満たす局所可積分な非負関数と仮定されますが、実際のデータでは年齢が\(0, \omega\)に近いところで急激に増加するU字型のグラフを示します。
関数\(\beta(a)\)を、単位時間あたりに微小な年齢区間\([a, a + da]\)にいる一個体から産まれる新生児の数とし、これを年齢別出生率(age-specific birth rate)と呼びます。積分$$\int_{a_1}^{a_2} \beta(a) p(t,a) \, da$$は単位時間に年齢\(a_1\)以上\(a_2\)以下の個体から産まれる新生児の数です。単位時間に産まれる新生児の総数\(B(t)\)は、$$B(t) = \int_0^\omega \beta(a) p(t,a) \, da$$によって与えられます。再生方程式は、この\(B(t)\)についての積分方程式です。実際のモデル解析においては、\(\beta\)はある有界閉区間\([\beta_1,\beta_2]\)の外で0である非負有界可積分関数と仮定されるのが通例です。区間\([\beta_1,\beta_2]\)は子どもを産める年齢の区間を意味し、再生産年齢(reproduction age)と呼ばれます。

方程式の導出

時刻\(t\)における\(a\)歳以下の人口を\(F(t,a)\)とします。すなわち、$$F(t,a) := \int_0^a p(t,s)\, ds$$とします。すると、時間間隔\([t, t+h]\)における人口の増減を考えることにより、人口方程式
$$ \begin{align}F(t+h,a+h) - F(t,a) &= ( [t, t+h] で生まれた人数 ) - ( [t, t+h] で亡くなった人数 ) \\ &= \int_t^{t+h} B(s)\, ds - \int_0^h \!\!\! \int_0^{a+s} \mu(b) p(t+s,b)\, dbds \end{align}$$を得ます。両辺を\(h\)について微分して\(h=0\)とおくと、積分記号下での微分を形式的に用いて、連続方程式$$p(t,a) + \int_0^a \frac{\partial p(t,s)}{\partial t}\, ds = B(t) - \int_0^a \mu(t,s)p(t,s)\, ds$$が得られます。特に\(a=0\)とおけば、$$p(0,t) = B(t) = \int_0^\omega \beta(a)p(t,a)\, da$$を得ます。連続方程式を\(a\)について微分すれば、$$\frac{\partial p(t,a)}{\partial t} + \frac{\partial p(t,a)}{\partial a} = -\mu(a) p(t,a)$$となります。これをマッケンドリック方程式(McKendrick equation)と呼びます。初期条件を\(p(a,0) = p_0(a)\)と設定することにより、年齢密度関数\(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), \quad t>0,\ a>0 \\ & p(0,t) = \int_0^\omega \beta(a)p(t,a)\, da,\quad t>0 \\ & p(0,a) = p_0(a),\quad a \in [0,\omega] \end{align} \right.$$これが安定人口モデル(stable population model)と呼ばれるものです。

モデルの解析

解の存在と一意性、人口作用素の性質

以下では、1980年代半ばに大きく発展を遂げた安定人口モデルへの関数解析的アプローチを概観します。安定人口モデルを\(L^1\)における抽象的な微分方程式として捉え直せば、Hille-Yosidaの半群の生成定理を用いて解作用素を導くことができます。本記事では線形の安定人口モデルのみを扱っていますが、この手法が真に効力を発揮するのは、モデルの右辺がより複雑な非線形問題を考えるときです。実際、1980年代以降の構造化人口動態学の発展は半群の理論を駆使した関数解析的手法の開発によるところが大きいのです。
いま、安定人口モデルにおける死亡力と年齢別出生率は\(\mu, \beta \in L_+^\infty (0,\omega) \)をみたすと仮定します。すなわち、$$ \bar\mu := \sup_{a \in [0,\omega]} \mu(a) < \infty , \ \bar\beta := \sup_{a \in [0,\omega]} \beta(a) < \infty $$が成り立つとします。このとき、安定人口モデルは次のようなBanach空間\(L^1(0,\omega)\)上の抽象的Cauchy問題として考えることができます: $$\left\{ \begin{align} & \frac{dp(t)}{dt} = Ap(t) \\ & p(0) = p_0 \end{align} \right.$$ ここで、\(A\)は以下のように定義される人口作用素(population operator)です: $$ \left\{ \begin{align} & (A\phi)(a) = -\frac{d\phi (a)}{da} - \mu(a)\phi(a) \\ & \mathcal{D}(A) = \left\{ \phi \in L^1(0,\omega) \,\middle| \, \phi は絶対連続かつ \phi(0) = \int_0^\omega \beta(a) \phi(a) \, da \right\} \end{align} \right. $$ \(\phi \in \mathcal{D}(A)\)ならば、ほとんどいたるところ微分係数\(\phi^\prime\)が存在して\(\phi^\prime \in L^1(0,\omega)\)となります。
人口作用素は以下の性質をもちます:

定理1.
人口作用素\(A\)に関して、以下の(1)-(6)が成り立つ。ただし、$$\hat\Psi := \int_0^\omega e^{-\lambda a} \beta(a) \ell(a)\, da , \quad \underline{\mu} := \inf_{a \in [0,\omega] \mu(a)}$$と定義しておく。
(1) \( \sigma(A) = P_\sigma(A) = \{ \lambda \in \mathbb{C} \mid \hat\Psi(\lambda) = 1 \}\) である。ここで、\(\sigma(A), P_\sigma(A)\)はそれぞれ\(A\)のスペクトル集合、点スペクトルを意味する。
(2) レゾルベント \(R(\lambda, A)\)はコンパクト作用素であって、$$\begin{align}(R(\lambda,A)\phi)(a) &= \int_0^a e^{-\lambda (a-z)} \frac{\ell(a)}{\ell(z)}\phi(z)\, dz \\ & + \frac{e^{-\lambda a}\ell(a)}{1-\hat\Psi(\lambda)} \int_0^\omega \!\!\! \int_z^\omega e^{-\lambda (a-z)}\frac{\ell(a)}{\ell(z)} \beta(a)\, da \phi(z) \, dz \end{align}$$ (3)\(\lambda \in \sigma(A)\)であれば\(\lambda\)の代数的重複度は有限であって、\(\lambda\)に対応する固有空間は\(N(\lambda-A) = \{ce^{-\lambda a} \mid c \in \mathbb{C} \}\)で与えられる。
(4) \( \rho(A) \supset \{\lambda \in \mathbb{C} \mid \Re \lambda < \bar\beta - \underline{\mu}\}\)であり、\(\lambda \in \mathbb{R}\)で\(\lambda>\bar\beta-\underline{\mu}\)であれば\(R(\lambda,A)\)は正値作用素である。
(5) \(\lambda\in\sigma(A)\)ならば\(\bar\lambda\in\sigma(A)\)である。
(6) \( ^\forall \alpha > -\infty \)について、\(A\)は半平面\(\Re \lambda \geq \alpha\)に高々有限個の固有値しかもたない。

定理2.
人口作用素\(A\)は閉作用素であり、かつ\(\mathcal{D}(A)\)は\(L^1(0,\omega)\)で稠密である。

以上の準備のもとで、Hille-Yosidaの定理により次が成り立ちます:

定理3.
人口作用素\(A\)は\(C_0\)半群 \(T(t)\ (t\leq 0) \)を生成する。\(T(t)\)は以下をみたす:$$ \begin{align} & \|T(t)\| \leq e^{(\bar\beta - \underline{\mu}) t} \\ & T(t)(L_+^1) \subset L_+^1 \end{align}$$しかも、\(p_0 \in \mathcal{D}(A)\)ならば\(p(t)=T(t)p_0\)は抽象的Cauchy問題の古典解を与える。

人口作用素のスペクトルとSharpe-Lotka-Fellerの定理

\(C_0\)半群とその無限小生成作用素のスペクトルに関する一般的な性質について少し述べます。作用素\(A\)のスペクトル上限\(s(A)\)とスペクトル半径\(r(A)\)はそれぞれ $$\begin{align} & s(A) := \sup\{ \Re \lambda \mid \lambda \in \sigma(A) \} \\ & r(A) := \{ |\lambda| \: \mid \lambda \in \sigma(A) \} \end{align}$$ また、半群\(T(t)\)の成長上限(growth bound)は$$\omega_0 := \inf \{w \in \mathbb{R} \mid ^\forall t\leq 0,\ ^\exists M(w)\in \mathbb{R}_+,\ s.t.\ \|T(t)\| \leq M(w)e^{wt} \}$$ と定められます。このとき、$$\omega_0 = \lim_{t\to \infty} \frac{\log\|T(t)\|}{t}$$であり、\(-\infty \leq \omega_0 < \infty\)が成り立ちます。

定理4.
\(\omega_0\)を\(C_0\)半群\(T(t)(t\geq 0)\)の成長上限、\(A\)を無限小生成作用素とすれば、\(^\forall t\geq 0\)に対して以下が成り立つ:$$e^{s(A) t}\leq r(T(t)) = e^{\omega_0 t}$$ \(A\)が有界作用素であれば、$$e^{s(A) t} = r(T(t)) = e^{\omega_0 t}$$とくに、\(\omega_0 = s(A)\)である。

\(\omega_0 = s(A)\)であるならば、無限小生成作用素のスペクトルの実部の上限が負であれば、\(\lim_{t\to \infty}\|T(t)\| = 0\)となります。これは、連立微分方程式\(x^\prime (t) = Ax(t)\)の解の漸近挙動を調べるには、行列\(A\)の固有値の実部の正負に注目すればよい……という議論の無限次元バージョンとなっています。
\(t_0\)が存在して\(T(t_0)\)がコンパクト作用素になるとき、半群\(T(t)\)は終局的コンパクト(eventually compact)であるといいます。終局的コンパクトな半群に関して、以下の定理が成り立ちます:

定理5. (スペクトル写像定理) $$\sigma(T(t)) \setminus \{0\} = e^{t\sigma(A)} := \{e^{t\lambda} \mid \lambda \in \sigma(A)\} , \quad t\geq 0 .$$ 定理3によって存在が保証された半群は人口半群(population semigroup)と呼ばれますが、これは終局的コンパクトであることが知られています。よってスペクトル写像定理が成り立ち、そのとき\(\omega_0 = s(A)\)なので、抽象的Cauchy問題のゼロ解が大域的安定である(すなわち、\(t\to\infty\)とすると人口がゼロに近づいていく)ためには、人口作用素\(A\)のスペクトル上限が負であれば十分であることがわかります。一般に、そうなっているかどうかは\(\mu\)と\(\beta\)のバランスによって決定されます。
ゼロ解の大域安定性について、より精密な条件を与えるのがSharpe-Lotka-Fellerの定理です。この定理によって、特性方程式\(\hat\Psi(\lambda) = 1\)のただ一つの実根\(\lambda_0\)によって安定人口モデルの漸近挙動が支配されることが判明するのです。ここでは、関数解析的な言葉を使って記しておきます:

定理6.
\(\lambda_0\)を人口作用素の実固有値とする。このとき、1次元の射影子\(P_0\)と正数\(\varepsilon >0, M(\varepsilon)\geq 1\)が存在して、$$ \| e^{-\lambda_0 t} T(t) - P_0\| \leq M(\varepsilon)e^{-\epsilon t}$$となる。ここで、\(u_0,v_0\)を\(A\)の固有値\(\lambda_0\)に対応する固有ベクトルおよび共役固有ベクトルとすれば、$$P_0\psi = \frac{\langle v_0,\psi \rangle}{\langle v_0,u_0 \rangle}u_0$$と表される。

つまり、半群\(T(t)\)は初期値によって決まる定ベクトルの\(e^{\lambda_0 t}\)倍に近づいていきます。\(\lambda_0\)は自然成長率(intrinsic rate of natural increase)あるいは安定人口成長率、内的成長率などと呼ばれています。\(\lambda_0<0\)ならば人口はゼロ解に近づいていき、\(\lambda_0>0\)ならば増大していきます。

Lotkaの再生方程式と基本再生産数

安定人口モデルを特性線\(t-a=Const.\)に沿って積分することにより、以下の表示を得ます:$$p(t,a) = \left\{ \begin{array}[ll] \, B(t-a) \ell(a), & t-a>0 \\ p_0(a-t) \frac{\ell(a)}{\ell(a-t)}, & a-t>0 \end{array} \right.$$したがって、出生数\(B(t)\)が求まれば年齢密度関数も求まります。これを境界条件に代入することにより、以下の積分方程式を得ます:$$B(t) = G(t) + \int_0^t \Psi(a)B(t-a)\, da$$ ここで、\(G(t), \Psi(a)\)は以下によって与えられる既知関数です:$$G(t)= \int_t^\omega \beta(a) \frac{\ell(a)}{\ell(a-t)}p_0(a-t)\, da , \quad \Psi(a) = \beta(a)\ell(a)$$特に、\(G(t)\)は時刻\(0\)に生存している個体から時刻\(t\)に産まれる単位時間あたりの新生児数と解釈されます。これが、Lotkaの積分方程式、あるいは再生方程式と呼ばれているものです。定理6から、\(B(t)\)に関するSharpe-Lotka-Fellerの定理が帰結します:

定理7. (Sharpe-Lotka-Fellerの定理)
安定人口モデルにおける再生方程式の解\(B(t)\)に関して、以下が成り立つ:$$B(t) = q_0 e^{\lambda_0 t}(1+O(e^{\eta t}))$$ここで\(\eta\)はある正数であって、$$q_0 = \frac{\int_0^\infty e^{-\lambda_0 t}G(t)\, dt}{\int_0^\infty a e^{-\lambda_0 a}\Psi(a)\, da}$$である。

さて、関数\(B_n(t)\ (n=0,1,2,\dots)\)を次のように定義してみましょう:$$\left\{ \begin{align} &B_0(t) = G(t)\\ &B_n(t) = \int_0^t \Psi(t-a) B_{n-1}(a)\, da ,\quad n\geq 1 \end{align} \right.$$関数\(B_n(t)\)は、初期人口の第\(n\)世代の子孫の出生率を与えています。各世代にわたる和\(\sum_{n=0}^\infty B_n(t)\)はLotkaの積分方程式の解となっています。
各世代の総出生数を計算してみましょう。漸化式から、$$\begin{align} \int_0^\infty B_n(t)\, dt &= \int_0^\infty \, dt \int_0^t \Psi(t-a) B_{n-1}(a)\, da \\ &= \int_0^\infty B_{n-1}(a)\, da \int_a^\infty \Psi(t-a) \, dt = \mathcal{R}_0 \int_0^\infty B_{n-1}(t)\, dt \end{align}$$ を得ます。ここで、$$\mathcal{R}_0 = \int_0^\omega \Psi(a) \, da$$は基本再生産数(basic reproduction number)といい、一個体が生涯に産む個体数と解釈される量です。人口学的な意味から、\(\mathcal{R}_0>1\)ならば人口は増大し、\(\mathcal{R}_0 <1\)ならば人口は減少すると予想されますが、実際、自然成長率は特性方程式$$\int_0^\omega e^{-\lambda_0 a}\Psi(a)\, da = 1$$をみたす実数であったので、符号関係$$\mathrm{sgn}(\mathcal{R}_0-1) = \mathrm{sgn} \lambda_0$$が成立します。基本再生産数は安定人口モデルだけではなく、その拡張のモデルや感染症モデルなどの様々な個体群動態モデルに現れており、構造化人口動態学における最も重要な概念の一つとして知られています。

参考文献

[1] W. Feller (1941), "On the integral equation of renewal theory", Ann. Math. Stat. 12: 243-267.

[2] 稲葉 寿 (2002), "数理人口学", 東京大学出版会.

[3] J. A. J. Metz, O. Diekmann (1986), "The Dynamics of Physiologically Structured Populations", Lect. Notes in Biomath. 68, Springer, Berlin.

[4] F. R. Sharpe, A. J. Lotka (1911), "A problem in age-distribution", Philosophical Magazine, Series 6, Vol. 21: 435-438.

[5] G. F. Webb (1985), "Theory of Nonlinear Age-Dependent Population Dynamics", Marcel Dekker, New York and Basel.