2021.10.21
『感染症の数理モデル』 稲葉寿編著（培風館）
これは日本語で読める唯一の教科書である。2008年初版の後、しばらく絶版であった。去年 COVID-19 の流行で必要となって、その部分を追加して2020年12月に増補版になった。この分野の日本での創始者、稲葉寿氏が編集して、第１章で数理モデルの一般論（解析的議論）を書いて、それを引き継いでいる格好の西浦博氏がデータ解析の側面を書いている。これらに引き続いて、個別の解析手法をいろんな人が書いている。５９００円もする本なので、図書館で借りて勉強中である。なかなか大変なので、まずは、基礎となる最初の２章分のメモを書いておく。

● 第１章：基本再生産数 R0 と閾値原理（稲葉寿）

・・・第１章で気づいたのは、SIRモデルやSEIRモデルでの変数は人口密度とした方が自然だということである。今では殆どのモデルで人口の絶対値、あるいは全人口に対する比率を変数としているので、感染係数を地域間で比較する時に注意を要する。福岡と広島の比較においても、再生産数への人々の活動度の影響を見る場合に、人口密度が高ければそれだけ影響も大きいと考えなくてはならないだろう。ただ、歴史上、いろんな感染症で基本再生産数を比較する時にそのことはあまり気に留められていないのは、感染が起こる場面における人口密度というのはそんなに違わないということなのだろう。具体的には、都市生活と家庭内や施設内である。また、集団免疫が問題になるほどまでに感染が拡がっていなければ、そもそも感染者数（密度）に定数を掛けても時間変化は変わらないという事もある。

・・・状態間遷移の感染時間（感染からの経過時間）や滞在時間（その状態に滞在している時間）依存性も含めた議論になっている処はさすがと思った。それらの依存性が無い場合として、待機時間依存性が指数関数となり、方程式が常微分方程式になる、としている。その背景の元で、1.3.5では未発症感染モデルを議論している。これは AIDS を念頭に置いている。

1.1.1 Kermack-McKendrick モデル I・・感染率 β、回復率 γ

・・・よく知られているが、一応方程式を書いておく。

      dS(t)/dt = -βI(t)
     ｄI(t)/dt = βI(t) - γI(t)
      dR(t)/dt = γI(t)

dI/dt = 0 となる条件は、感受性人口密度 S(t)=γ/β である。これを Scr と表す。
初期感染者数 I(0) と初期感受性人口 S(0)=N-I(0) から始まって、Scr で I が最大となり、やがて I=0 になる。これが SIRモデルの相図 図1.1 である。
 ζ=βI(0)/γ と置き、初期の感受性人口から最終的に除去される人口の割合を p とすると、

      1 - p = exp(-R0p-ζ)

が成り立ち、ζ→0 の極限では

      1 - p = exp(-R0p)

となる。この解 p が最終的に感染する人口割合であるから、ここから逆に R0 が判る。
図1.2 p と R0 の関係、図1.3 1978年インフルエンザの流行データ、図1.4 ペストの死者数とモデル比較

1.2.1 Kermack-McKendrick モデル II

・・・Kermack-McKendrick は SIR モデルの発案者としてよく知られているが、現在一般的に使われているのはその内でもっとも単純なケース I であって、引き続く論文で、より一般的な方程式を提案している。感染してからの経過時間 τ 別に感染者数を変数 i(t,τ) として、更に感染率や治癒率も τ に依存する β(τ)、γ(τ) としている。1927年に提案されたものの、その数学的構造が解明されたのは1970年代になってからなので、あまり知られていない。感染初期近似として線形化すれば、稲葉氏の求めた一般解が得られる。

      t > τ では、i(t,τ)=Γ(τ)v(t-τ)、
      t < τ では、{Γ(τ)/Γ(τ-t)}i0(τ-t)
 ここで i0 は初期条件である。

  Γ(τ) は感染齢 τ において感染状態に留まっている割合である。具体的には、

      Γ(τ)=exp{-∫γ(σ)dσ}：σ 積分は 0 から τ まで。
 v(t) は i(t,0) つまり単位時間あたりに発生する新規感染者数である。

      v(t)=g(t) + ∫Ψ(τ)v(t-τ)dτ：積分は 0 から t まで、
      Ψ(τ)=S0β(τ)Γ(τ)：感染齢 τ での感染係数×人口。
      g(t)=S0∫β(τ){Γ(τ)/Γ(τ-t)}i0(τ-t)dτ： 積分は t から ∞ まで。
g(t) は初期感染者 i0 から単位時間あたりに発生する二次感染者数である。

感染齢が長くなればいずれ感染性が無くなる（β→0）から、その最大時間を ω とする。
t > ω となれば、g= 0 になってしまい、初期条件が影響しなくなる。
v(t)は同次方程式に従うので、指数関数型の解を持つ。

      v(t)=kexp(λt)=∫Ψ(τ)kexp(λ(t-τ)dτ  より、
      1=∫Ψ(τ)exp(-λτ)dτ： 積分は 0 から∞。
これをオイラー・ロトカの特性方程式という。

λ=0 の時の右辺が基本再生産数 R0 である。

      Φ(τ)=Ψ(τ)/R0
は二次感染発生までの待機時間（世代時間）の確率密度関数である。
Φ(τ) が平均 τ0、分散 σ^2 の正規分布で近似できるとすれば、

      (σ^2/2) λ^2 - τ0 λ + log(R0) = 0
という二次方程式の解として λ が与えられる。
人口学のダブリン・ロトカの公式が使えて、

      λ0 ≒ {τ0 - (τ0^2 - 2σ^2log(R0))^(1/2)}/σ^2
λ0 から R0 が求められる。

      R0 ≒ exp(λ0τ0 - σ^2λ0^2/2)

（注：僕の最初の論文 も、大枠で言えば、この方程式の枠組みの中に入っている。
β(τ) を固定された待機時間（2日）に始まって隔離時間に終わる指数関数として、
γ(τ) をその待機時間に始まる定数とすればよい。更に、検査・隔離の項が時間遅れで γ(t) に追加されるというモデルである。
遅延微分方程式とは言っても、滞在時間に応じて変数を設定すれば、普通の偏微分方程式になる、ということである。）

年齢・性別だとかも含めて、一般的に人口を分割して論じることもできる。これらは人口学の分野で開発されてきたものの借用である。

1.3.5 未発症感染モデル

・・・SEIR モデルにおいて、待機状態 E においても２次感染を起こすとしたモデルである。だから、これは待機状態というよりも感染性潜伏期間ということである。
待機状態 E から感染状態 I への遷移の速度係数を ε とする。感染が両方から起きるので、それぞれの感染係数を β1、β2 としている。感染すると E が増える。
誕生と自然死亡速度、b と μ が考慮されているが、ややこしいので無視すると、未発症感染者数密度、発症感染者数密度、x(t)、y(t) に対して、

      dx/dt = -εx + β1x + β2y
      dy/dt = -γy + εx
という方程式になる。２つの感染状態間の次世代行列は

       K = (R1  R2)
               ( 1   0)
      R1 = β1/ε、R2 = β2/γ、となる。
K の固有値から

      R0 = {R1 + √(R1^2+4R2)}/2
実態としての基本再生産数は T1 であって、

      T1 = R1 + R2
となる。R0 も T1 も 1 を境界として感染の広がりが決まる。
R1 < 1 の時は、発症者による２次感染によって感染が拡がり得る。
この場合には、発症者から発症者への（タイプ別）再生産数は

      T2 = R2/(1-R1)
となる。もしも、発症者を e の比率で隔離したとすると、

      (1-e)T2 < 1
の時に感染が終息する。閾値 e* は

      e* = 1 - 1/T2 = (1 - 1/T1)/(1 - R1/T1)

となる。(1 - 1/T1) は集団ワクチン接種での臨界免疫割合、
(1 - R1/T1) は全２次感染者中で発症感染者の割合である。
・・・ややこしいので、すべて R1 と R2 で表すと、

      ε* = (R1 + R2 - 1)/R2
となり、

      R1 + (1 - ε*)R2 = 1
である。R1 は未発症感染者の寄与で、R2 は発症感染者の寄与であり、
発症者の ε* の割合を隔離することで、再生産数 R1 + R2 を 1 まで落とす、という意味。
TOST の図で言えば、時間軸負の部分が R1/R0、正の部分が R2/R0 である。
特に新規なことを言っているわけではないが、僕の２番目の論文 へのモデル的な見方を与えている。

1.5.2 リスク構造モデル

・・・HIV/AIDS 等性的感染症の流行モデルで、性的接触時の感染確率を β とするが、感染確率はこれに一人当たりの性的接触頻度 ζ がかかる。この接触頻度は分布を持つから、その平均と分散を m、σ^2 とする。μ、γ をそれぞれ自然死亡率、隔離率とすると、基本再生産数が

      R0 = (βm/(μ+γ)){1+(σ/m)^2}
となる。接触頻度の平均が同じでも、その分散が大きいと R0 が大きくなり、分布がベキ法則にしたがうような場合（スケールフリーネットワーク）には R0 が∞になりうる。実際 AIDS では ３乗則と言われており、その場合には、区間[a,b]であったとすると、分散が log(b/a) に比例するから、実際上感染侵入を防げない。

（注：COVID-19 においても人々の活動度には大きな分散があるので、こういう事情で R0 が大きくなっていて、その活動度の大きな人たちが感染して免疫を獲得してしまえば、自然収束すると言えるのかもしれない。第６章と第１０章に関係してくる。）

● 第２章：感染症数理モデルのデータサイエンス（西浦博）
この章には下記の付録が付いている。
http://plaza.umin.ac.jp/~infepi/appendix.pdf
・・・伝染病データは通常の医療データとは異なり、患者同士の相互作用が大きいので、リスクが状況依存している。極端な例が集団免疫である。

2.2 内的動態の推定問題
・・・数理モデルにおけるパラメータを感染のデータから推定する方法である。遷移パラメータ（状態間の遷移速度が質量作用の法則にしたがうという仮定の時の係数）の逆数は時間の次元を持ち、その感染症を特徴づける。感染待ち時間（latent period）、感染性期間（infectious period）、潜伏期間（incubation period）、症候性期間（symptomatic period）、発症間隔（serial interval）、世代時間（generation time）。これらの図示が 図2.2 である。

2.2.1 潜伏期間の推定
・・・しばしば対数正規分布が仮定される。

2.2.2 病気齢に対する２次感染相対強度の推定
・・・伝搬ネットワークのデータと潜伏期間の分布を畳み込みで組み合わせて、逆計算する。ただし、接触頻度が病気齢に対して独立でなければ、とりわけ発症後については過小評価となる。

2.3 家庭内感染
・・・家庭内感染では接触が濃厚であって、しかも機関が明確な場合が多いので、解析が容易となる。発症率（attack rate）は接触した場合に感染する割合である。以下、実際上役に立ちそうな解析方法である。

2.3.1 家庭内感染確率の推定
・・・
      確率 B：感受性宿主が一つの流行を通じで（家庭外）コミュニティから感染しない確率
      確率 Q：感受性宿主が一人の家庭内感染性宿主から感染しない確率
      Pr(k|s)： s 人の家庭において、k 人が感染する確率。
      Pr(1|1)= 1-B、Pr(1|2)=2(1-B)BQ= 2Pr(1|1)BQ、
      Pr(2|2)=(1-b)^2 + 2(1-B)BQ= 1 - Pr(0|2) - Pr(1|2)
一般化すると、
      Pr(k|s) = sCk Pr(k|k) B^(s-k) Q^(k(s-k))  : k=0,1,,,s-1
      Pr(k|k) = 1 - ∑ Pr(k|s) : 和は k=0,,,s-1

実際の家庭内感染データをこの式に当てはめることで、 B や Q の推定値が得られる。
家庭内濃厚接触において一日あたり感染する確率を β とすれば、
感性性の持続日数 m を使って、

      Q = (1 - β)^m
となるので、感染性についての情報が得られ、B と Q の比較から、家庭外での接触状況が推定できる。

2.3.2 家庭内の観察データにおける発症間隔の利用

・・・家庭内で二人の感染者が出た場合を集めて、それらの発症間隔の分布を s(t) とすると、

      s(t) ＝ ph(t) + (1-p)l(t)
と表せる。h(t) は家庭外からの感染者同士の発症間隔であり、 l(t) は家庭内での感染による。
家庭外からの感染タイミングはランダムと考えれば、h(Δt) は
潜伏期間分布 f(t) において、時間差 Δt の持つ分布である（∝∑ f(t）f(t+Δt) ：和は t について）。
l(t) は本来の意味での発症間隔であるから、q を病気齢に対する二次感染の確率分布とすれば、

      l(t) = ∫q(σ)f(t-σ)dσ （σ 積分は t -xからt）
となる。つまり、家庭内感染データの分析と潜伏期間の情報から二次感染のタイミング分布が推定できる。

2.4 基本再生産数の推定

・・・R0 は感染に対する介入が無い場合の初期再生産数であるが、その推定には仮定が必要であり、どのような仮定に基づいているかを吟味すべきものである。

2.4.1 流行（Epidemic）データによる推定

・・・指数関数的感染者数の振る舞い（マルサス的）が見られるので、SIRモデルでは、その増殖率に感染期間(1/γ)を掛ければ R0-1 が得られる。

      I(t) = I(0)exp{Λt} = I(0)exp{(R0-1)γt}

不顕性感染者が多い場合は推定が難しい。
スーパースプレッダーなど個体レベルの異質性や集団レベルの異質性には対応できない。
I(t) は発症者数ではなくて、二次感染させるような感染者数であることに注意すべきである。

・・・SEIRモデルでは、
      R0 = 1 + VΛ + f(1-f)(VΛ)^2

Λ は実測データから得られた指数関数の係数、V は世代時間＝非感染期間+ 1/γ であり、f=1/γV である。
この式は SARS の流行期に提案された。第３項は SIR モデルによる推定に追加されている。

2.4.2 流行（Epidemic）データによる推定（その２）

・・・SIRモデルが成り立つとすれば、感染初期ではなく、感染が終了した後の最終規模（どれくらいの人口比が感染したか）によっても R0 の推定ができるが、これは感染対策をしていないことが前提となるので、人間社会の感染症に対してはあまり有用ではない。

・・・ここでは確率モデル（計数プロセス）を議論している。人口 n+1 の中の一人が初期感染者である。βの定義が少し違う。

      dI(t)/dt = -γI(t)+β(S(t)/n)I(t)<br>K(t)=S(0)-S(t)= I(t)+R(t)-1 と R(t) を変数として考える。

確率過程を考えて、これらが増加する確率を与える。

     P(dK=1,dR=0|)=β(S(t)/n)I(t)dt
     P(dK=0,dR=1|)=γI(t)dt
     P(dK=1,dR=0|)={1-β(S(t)/n)I(t)-γI(t)}dt

u を過去の時刻とすると、u における増加確率を積分した量を考えることができるが、
これは必ずしも時刻 t での実現量ではない。しかし、それらの期待値は等しい。
これがマルチンゲールプロセスである。つまり、次の確率変数は平均値が 0 となる。

      M1(t) = K(t) - ∫β(S(u)/n)I(u)du
      M2(t) = R(t) - ∫γI(u)du    :積分は 0 から t まで

ここから先はよく判らない。マルチンゲール理論というのがあるらしい。
それを使うと、β/γ の推定値 Θ について、
その期待値が -log(1-p)/R(T) となるが、
確率過程なので、その分散も評価できるのが特徴であり、感染者数が少ない時に効力を発揮する。
R0 の期待値は n×Θの期待値、R0 の標準偏差は Θ の標準偏差×n となる。
Θ の標準偏差は {1/(S(T)+0.5)-1/(S(0)+0.5)-Θの期待値の自乗×R(T)}^(1/2)/R(T)

2.4.4 単純分岐過程の基礎

・・・なかなか難しい。確率過程で考えると平均的感染状況だけでなく、そのバラツキも議論出来て、感染の完全終息確率等が判る。

・・・感染の世代を考え、第 n 世代の感染者数を Cn とする。
各感染者からの二次感染者発生数を ρi とする。その分布は感染者に依存しないとする。

      pk = Pr(ρi=k)
      R0 = E(ρi) = ∑kpk :和は k の 1以上
      Var(ρi) = ∑(k-R0)^2
Galton-Watson 型分岐過程という。

第 n 世代の感染者数の期待値 E(Cn)=Mn、分散 Var(Cn)=Vn が問題となる。
確率母関数 gρ で議論できる。

      gρ(s) = E(s^ρ) = ∑ pk s^k ：和は k=0 から∞
R0 は s=1 における s についての一回微分、
感染者がこれ以上二次感染者を生み出さない確率は g(0)： 絶滅確率。
第 0 世代は感染者が１人なので、特別に gC0=s となる。
それ以降の g は下記漸化式から得られる。

      g_n+1(s) = E(s^C_n+1) 
                     = ∑E(s^(ρ1+ +ρi)|Cn=i) Pr(Cn=i) ：和は i=0 から∞
                     = ∑g_(ρ1+ +ρi)(s) Pr(Cn=i) ：和は i=0 から∞
                     = ∑g_ρ1(s)・・・g_ρi(s) Pr(Cn=i) ：和は i=0 から∞
                     = ∑g(s)^i Pr(Cn=i) 
                     = g_n(g(s))
                     = g(g_n(s))

pk がポワッソン分布とし、平均再生産数 v の時の２次感染者数の密度関数を f_v(u) とすると、

      g(s) = ∫exp{-u(1-s)}f_v(u)du ：積分は 0 から∞

流行初期とすれば、S は一定として、v=R0 と置いて、

      g(s) = exp{-R0(1-s)}

ポワッソン分布であることを利用すると、観察データ Cn となる確率（尤度）が得られる。

      Π(R0C_k-1)^Ckexp(-R0C_k-1) ：積は k=1 から n まで

尤度最大とすれば、R0 の推定ができる。
R0 の推定量は (Ck/C0)^(1/k)
g を微分すれば期待値であるから、

      Mn = R0M_n-1 =  = R0^nM0
分散は
      R0=1 の時、Vn = nVar(ρ)
      R0>1 の時、Vn = R0^(n-1) Var(ρ)(R0^n - 1)/(R0 - 1)

・・・pk がばらつきパラメータ k のガンマ分布の場合には、
      g(s) = {1 + (R/k)(1-s)}^(-k)
となり、負の二項分布となる（スーパースプレッダー）。SARS、COVID-19。

2.4.5 単純分岐過程を利用した絶滅の分析

・・・確率的な推移モデルからは、R0>1 における感染終焉（絶滅）の確率 q が得られる。
      q=g(q)
の最小解である。但し、初期感染者が a 名の場合は q^a となる。

2.5 実効再生産数の推定
・・・２種類ある。
症例再生産数 Rc(t) において t は感染時刻である。つまり時刻 t に感染した感染者が生み出す２次感染者数の期待値である。感染対策（接触回避）との相関で言えば、あまり有用でない。瞬間再生産数 Rt(t) において、t は２次感染者の感染時刻であるから、感染対策の効果を評価するのに適している。

・・・
      Rt(t) = i(t)/{∫i(t-σ)s(σ)dσ} ：積分は 0 から∞
となる。s(t) は世代時間、近似的には発症間隔の分布である。
i(t) は時刻 t において新たに感染した感染者数（単位時間あたり）である。
これはデータから直接得られる情報ではないので、推定する必要がある。
c(t) を時刻 t において新たに発症する感染者数（単位時間あたり）とすれば、
潜伏期間分布を f として、
      c(t) = ∫f(τ)i(t-τ)dτ ：積分は 0 から∞
なので、逆計算で推定される。

・・・この辺は馴染みがある。ただ、去年の５月頃までは感染症の専門家の間でしか知られていなかった。