著者：畝見達夫，作成：令和3年8月18日，編集：10月10日
この文書は個体ベース感染シミュレータ SimEpidemic に使われているモデル ver. 1.9.0 について記述したものである。 過去のモデル仕様についてはつぎのリンク先を参照されたい。 1.8.6， 1.8.5， 1.8.4， 1.8.3， 1.8.2， 1.8， 1.7.3
シミュレーションの場となる世界は，連続値直行座標で位置を表現する正方形の2次元平面である。 初期状態では人口に対応する数の個体がランダムに配置される。 初期状態における個体の健康状態は，ランダムに選ばれた数人の未発症感染者を除き， 未感染であるものとする。 辺の長さと人口を調整することで人口密度を調整する。
シミュレーション過程は，まず，現在の状態での各個体の相互作用を計算し， その結果と個体自身の状態に応じて次の状態を算出するというステップを繰り返すことで実現する。
各個体に拠点位置を与え人口の空間分布を指定することで，人口密度の地域的偏りを表現する。 設定の選択肢としてつぎの3種類を用意する。
上記 3. の場合の各個体の自宅位置 \(H\) は， 区間 \([-1,1]\) の2つの独立な一様乱数 \(x, y\) から下式により決定する。
\[ H = \left(\frac{px+1}{2}S,\, \frac{py+1}{2}S\right)\] \[ p = \frac{q}{1-(1-q)\cdot\max (|x|, |y|)},\; q=\frac{a}{1-a}\]
ただし，\(S\) は世界の一辺の長さ，\(a\) は密度の偏りを表す区間 \((0,0.5]\) の係数である。 \(a=0.5\) のとき \(q=p=1\) となり一様分布になる。 \(a \rightarrow 0\) では中央の点に集中し， 極限では \(q=p=0\) で \(x, y\) の値にかかわらず \(H\) は中央位置 \((S/2, S/2)\) となる。
上記 4. は，与えられた画像データを 512x512 画素の正方形に比例変換し， 各画素の明度に比例した密度となるよう個体の拠点位置を定める。 アプリケーション版のシミュレータには，複数箇所の密集地点をランダムに配置した分布パターンを生成する機能ががあり，抽象的なモデルとしてこのパターンを利用することも有用であろう。 詳細は 人口空間分布の与え方と分布の生成 を参照のこと。
年齢は個人の特性として広く使われており，感染症に関連する統計指標も様々なものが利用可能である。 利用可能な実測データとの比較をもとにシミュレーションモデルのパラメータ設定を調整することは， 現実との乖離を小さくするために有効であろう。 対象地域の年齢構成を国勢調査などのデータから抽出し，各年代の比率に合うよう各個体に年齢を割り付ける。 年齢による行動の違いや発症機序の違いをモデル化する。 （年齢と行動の相関については未実装）
各個体は各ステップにおいて，周囲から受ける力ベクトルに従い移動する。 昼間は長距離移動や集会参加をするが，夜間には自宅位置へ戻るような行動をとる。 力には，周囲の個体との距離の２乗に反比例する斥力と周囲のもっとも好む個体への引力を合成して力のベクトルを計算し， 質量と摩擦係数を想定した積分の近似計算1により次の位置を算出する。 世界の境界に隣接する場合は境界線からの斥力も受ける。 個体 \(i\) が受ける力 \(F_i\) は以下の式で表される。 \[F_i = \alpha\cdot\sum_j\frac{p_i - p_j}{|p_i - p_j|^3}+\beta\cdot(p_b-p_i) +\gamma\cdot\sum_{j=1}^4\frac{p_i-w_j}{|p_i-w_j|^3}+B_i\] \(p_i\) は個体 \(i\) の位置，\(b\) は周囲で最も好みの個体，\(w_j\) は境界線 \(j\) 上の最も近い点，\(\alpha,\:\beta,\:\gamma\) は係数である2。 \(B_i\) は後に述べる集会の中心または自宅へ向かうベクトルである。
発症していない個体は設定された確率で世界の領域内の別の位置に移動する。 対策としての移動制限は，この確率と距離を調整することで実現する。 ただし，夜間に移動が発生した場合は，移動先を自宅位置とし，帰宅するものとする。
昼間の学校，職場，イベント，集会，さらに家庭など人が集まる現象を再現するモデルとして， 集会 を確率的に発生される。個々の集会は円形の領域で．大きさ，持続時間，および，強さが割り当てらる。 集会の円周付近の発症していない個体は，集会の中心へ向かう力を受ける。力の大きさは円周上で最も強く， そこからの距離に反比例して小さくなる。 右図に個体が受ける力の分布を示す。 上の式における \(B_i\) は以下に定義される集会の中心へ向かうベクトル \(G_i\) に 係数 \(\zeta\) を掛けたものになる。
\[ G_i=\begin{cases} 0 & \mbox{if } |p_i - G_p| > G_r + C_r\\ 1 - (|p_i - G_p|-G_r)\:/\:C_r & \mbox{else if }|p_i - G_p| > G_r \\ |p_i - G_p|\:/\:G_r & \mbox{otherwise}\\ \end{cases} \] ただし \(G_p\) はその個体から最も近い集会の中心，\(G_r\) はその半径，\(C_r\) は円周から影響する領域の境界までの距離を表す定数である。
さらに，集会に留まるよう中心から \(G_r\) 以内の領域では摩擦係数が増すものとし，中心で通常の倍の摩擦になるよう距離に比例した補正を適用する。
自宅位置がある場合には，ランダムに選択した個体の自宅位置を集会の中心位置とし， ない場合は，世界の範囲で一様乱数により集会の中心位置を決定する。 また，自宅位置に偏りがある場合は，人口密度に反比例するよう集会の大きさを変える。
夜間は個体は集会の影響を受けず，代わりに自宅へ向かう力を受ける。 このとき \(B_i\) は \[ B_i=\begin{cases} 0 & \mbox{if } |p_i - H_i| >\theta_H\\ \psi\cdot\theta_h & \mbox{else if }|p_i - H_i| > \theta_h \\ \psi \cdot|p_i - H_i| & \mbox{otherwise}\\ \end{cases} \] ただし，\(H_i\) は個体 \(i\) の自宅位置，\(\theta_H\) と \(\theta_h\) は距離の閾値， \(\psi\) は係数である。つまり，自宅から \(\theta_H\) 以上離れた個体は帰宅せず， 自宅からの距離が \(\theta_h\) から \(\theta_H\) の間である個体は自宅へ向かう一定の力を受け， \(\theta_h\) 以内にいる個体は距離に比例した力を受ける。 また，自宅から \(\theta_H\) 以上離れた個体が夜間に長距離移動する場合は， 行き先が任意の位置ではなく自宅位置になる。
各個体の健康状態を 1. 未感染，2. 未発症，3. 発症，4. 快復，5. 死亡の５種類に分類する。
未発症または発症中の個体 \(a\) と未感染の個体 \(b\) が距離 \(D_i\) 以内に近づくと， 設定された確率で \(a\) が感染しているウイルスが \(b\) に感染する。 確率は，距離，感染者の感染性，被感染者の免疫力に依存し， \(a\) から \(b\) への感染確率 \(P_{a\rightarrow b}\) は以下の式で表される。
\[ P_{a\rightarrow b} = P_i \cdot R_v \cdot C_a \cdot D_{ab} \cdot (1-I_b) \]
\(P_i\) は全体に適用されるパラメータであり，一般的な感染確率を表す。 マスク着用や手洗いの励行，季節性要因などを総合した値と仮定する。 \(R_v\) は変異ウイルスの種類に依存した感染率である。 感染者 \(a\) から他個体に感染する感染性の強さ \(C_a\) は，感染から感染性の遅れ期間を経るまでは 0 で， その後ピーク到達期間と潜伏期間の短い方までは，時間経過に比例して強まる。 \(D_{ab}\) は \(a\) と \(b\) の間の距離の効果であり，以下の式で定義される。
\[ D_{ab} = \min \left(1, \left(\frac{D_i \sqrt{R_v} - | p_a - p_b |}{2}\right)^2\right) \]
ここで \(D_i\) は感染可能な標準的な距離であり，\(D_i\sqrt{R_v}\) より離れた個体間では感染は起きない。 \(I_b\) は被感染者 \(b\) の免疫による予防効果である。 免疫は感染あるいはワクチン接種による獲得免疫を想定する。 未感染者がワクチンを接種済みの場合，あるいは快復者の場合は， 接種あるいは感染による獲得免疫効果 \(I_b\) に従って感染確率が低くなる。
感染した個体は個体に割り振られた速度で症状が進み，ある閾値を越えると発症し， さらに次の閾値を越えると死亡する。 感染後，個体に割り振られた快復開始日数に達した個体は，症状が一定速度で快方に向かい， 感染前の健康状態に戻った時点で完治したとみなす。 快復した個体は免疫を獲得したものとみなされるが， 設定された期間を過ぎると免疫を失ったとみなし，未感染の状態にもどる。
加齢による重症化リスクの増加は，快復開始期間の増加によるものと考え，
個体の年齢と快復開始期間の分布の間に相関があるものとして，
個体の年齢から快復開始期間を求める関数を用いて計算する。
関数は 指定した分布に従う乱数の生成アルゴリズム に基づき，
期間の最小値，最大値，最頻値からベル型の確率密度関数に従う値を出力する。
このときの最頻値 \(M_r(A_i)\) は
\[
M_r(A_i) = \zeta\; M_f \exp \left( \frac{A_i - 105}{\eta} \right)
\]
とし，最小値，最大値は最頻値に係数をかけた値とする。ここで，\(A_i\) は個体 \(i\) の年齢，
\(M_f\) は感染から死去までの期間の分布の最頻値，
\(\zeta\) および \(\eta\) はパラメータとして与えられる正の係数である。
右のグラフに \(\zeta = 1.5, \eta = 50, M_f = 21\) のときの，
\(M_r(A)\) を縦軸，\(A\) を横軸にして描いた曲線，および，
最大値の比を 3倍，最小値の比を 0.4 倍としたときの範囲を示す。
快復開始期間が潜伏期間より短ければ発症せずに感染は消滅し，
快復開始期間が死亡までの期間より長ければ死に至る。
また，治療法の開発により期間が短縮されるものと考え，治療効果の係数を乗じた値を用いる。
一度感染した個体は，感染により免疫を獲得するが， その効果の強さについては罹患時の最大重症度に依存するものと考え， 下記の図に示すように，快復前の最大重症度に応じた有効期間と再感染確率を個体に割り当てる。
重症化の比率が年齢と相関があることは早くから指摘されており，また，長距離移動や集会への参加などの行動の活発さも年齢に依存すると考えられる。 それが感染拡大の状況に及ぼす影響については不明であるが，終息を難しくしている1つの原因として，無症状のまま活発に行動する個体の存在も仮説として考えられる。 このような様相をシミュレーションに反映されるため， 発症機序における快復開始期間， 発症から死亡までの期間などの個体の生理的事象に関係する属性については年齢との相関性を考え， 質量，長距離移動，集会への参加など，社会的な行動に関係する属性については， 「活発さ」との相関性を考える。 相関性のある分布のに従う値の割り当てのについては， 範囲，中間値，尖度，他のパラメータとの相関性などを指定した乱数の発生アルゴリズムを用いる。 分布の形式をとるパラメータは，最小値，最大値，最頻値の３つの数値の組み合わせにより表現し， その条件を満足するベル型の確率密度関数に従った擬似乱数発生により個々の場合の値を決める3。 アルゴリズムの詳細はこちらを参照されたい。
社会的距離，移動の制限，隔離，接触者追跡を次節に述べる検査とともに設定する。
行動モデルの節で述べた他の個体との衝突回避の斥力の係数を変化させることで社会的距離を表現する。 パラメータとして，斥力の係数の増加率と要請への協力率を設定する。
行動モデルの節で述べた他の場所への移動の距離と頻度を調整することで，移動制限を表現する。
次節で述べる検査の結果が陽性であった個体を隔離する。 隔離された個体は世界とは別に設けられた病院領域へ移動する。 隔離された個体が快復すると元の位置へ戻る。
接触者追跡は，各ステップで感染可能な距離以内に近づいた別の個体を接触者として，
個体が持つリストに追記することで実現する。既にリストにある個体と再度接近した場合は，
過去の記録を消去し，新しい記録で書き換える。ただし，接触から14日を過ぎた記録は消去する。
接近したすべての個体を記録するのではなく，設定された確率により記録するかどうかが確率的に決まる。
検査の結果が陽性であった個体の接触者リストを参照し，1. 検査，2. ワクチン接種，3. 両方のいづれかの処置を施す。2. および 3. の場合は 追跡接種 が行われる。
ただし，ワクチン接種の待ち行列中の順位が前方に移動するだけであり，
実際に接種されるかどうかは，接種速度，および，順番が来た時点での感染の有無に依存する。
検査の対象となる個体の選択方法と確率，検査の感度と特異度を設定する。 対象とする個体は，症状の出た個体，疑いのある個体，陽性の検査結果が判明した個体との接触者 の中から選ぶ。接触者以外の対象は設定された確率により選択する。 接触者の検査確率は，接触者の記録の確率に吸収できるので，別にパラメータを設けない。
ワクチンを接種した個体の感染抑制効果と， ワクチン接種を実施する速度および優先順の設定について述べる。
ワクチンを接種すると，徐々に獲得免疫の効果により感染リスクが下がり，また，ある期間を過ぎると効果が薄まり， 感染リスクが徐々に元に戻るものと考える。ワクチン接種の時期と感染リスクの変化の関係を下の図のようにモデル化する。 Pfizer-BioNTech の治験情報4,5,6 に基づき，2回接種で1回目から3週間後に2回目を接種する場合などを設定できるようにする。 Moderna7 や AstraZeneca など他社製ワクチンについても，同様のモデルを用いる。 接種による獲得免疫効果の経時的減衰については， 他の国々に比べ接種が先行したイスラエルでの接種実績8 および 2021年5月のアルファ株に対する効果と 2021年6月のデルタ株に対する効果9 などから推測する。
接種を受けた個体の感染リスクの推移のモデルを下の図に示す。 同じ接種を受けても個人差は生じると考えられるが，ここでは感染は確率的に起きるものと想定するので， 下図は多数の個人についての平均とみなしてもシミュレーション上の支障はない。
ワクチンの直接的な効果は感染に対する免疫系の即応性である。 実際にはまったく感染しない確率が上昇するのではなく， 少量の感染でも素早い免疫応答により，検出できない程度にウイルスの排出量を低く抑えられると考えられる。 同様に，免疫応答よりも早くウイルスが体内に侵入した場合でも， 強化された免疫機能により症状の進行が遅くなり，快復も早くなると考えられ， 発症および重症化の防止が観測される効果として現れる。 これを再現するため，接種済みの個体が感染した場合には， 発症機序 としてモデル化された快復開始期間が短縮されるものとする。
接種速度と対象個体の優先順はワクチンの種類毎に指定する。
接種の実施数は，１日当たりの接種人数の人口に対する割合を示すパラメータ値により指定する。 接種により獲得された免疫の有効期間は個体間で共通のパラメータとして設定される。 つまり，個体差はないものと仮定する。
接種の優先順の選択肢として， 1. ランダム，2. 高齢者優先，3. 中央優先，4. 人口密度の高い地域を優先 の4種類を用意する。人口密度は拠点位置の分布によるため， 分布が一様な場合はランダムと，中央に集中した分布の場合は中央を優先と同じである。 ラインダム以外の順序については， 順位の遵守の程度を表顕するパラメータとして，順位から外れた個体に接種する確率を用意する。
追跡接種を行う場合は，陽性判明者の接触者を優先的に接種する。 いずれの場合も，未感染者と検査されていない無症状感染者を接種対象とする。 追跡接種に用いるワクチンの種類は１つを指定するものとする。
何らかの理由によりワクチンを接種できないまたはしない個体のモデルとして， その年齢層別に人口に対する割合を設定し，接種の順番が回ってきてもその個体には接種しないこととする。 ただし，そのような個体がいる場合でも1日あたりの接種数は変わらないものとする。 つまり，ワクチンの効果が消えなければ，接種しない個体以外の個体がワクチン接種を終えた時点で接種は終了する。
接種しないクラスタの影響を見るため，接種しない個体の地理的分布を指定できるよう， クラスタ数および配置を以下のようなパラメータにより指定する。 クラスタ内個体数とクラスタ外個体数の比、 クラスタ数とクラスタ内個体数の比、 クラスタの配置。詳細は 非接種クラスタ分布の生成 を参照のこと。
ウイルスの遺伝子は増殖するたびに絶えず変異を繰り返しているが， 進化生態系の基本的な動力学に従い増殖力の強い変異は保存され， より弱い遺伝子は駆逐される。ここでは，変異の機構をモデル化せずに， 既知の変異株が外部から移入されるものと仮定する。 また，変異株の種類による発症機序の差異は，基本的に感染者体内での増殖速度によるものと仮定する。
増殖速度に比例して，潜伏期間および死に至るまでの期間が早まるものと仮定する。 増殖力の上昇は感染者から排出されるウイルス量の増加につながる。 それによる感染力と毒性の変化も増殖速度に沿うものとし， 感染確率の式 における係数 \(R_v\) に影響する。 また，\(R_v\) が大きいほど症状の亢進が速くなるものと仮定し， 発症機序における潜伏期間，快復開始期間，感染性の遅れ，感染性ピーク到達期間のそれぞれについて， \(\sqrt[3]{R_v}\) で除した値を用いる。
個体の感染の有無を判定する検査においても， ウイルスの増殖が速ければ検体の含まれるウイルス量が増え， 偽陰性の確率は減る，つまり，感度は上昇すると考えられる。 ここでは，増殖力が 0 で感度 0%，無限大のときに感度が 100% となると考え， 以下の式で感度 \(Se\) を補正する。
\[ S_e = 1 - (1-S_1)^{R_v} \] ここで \(S_1\) は \(R_v=1\) のときの標準となる感度である。
接種したワクチンの種類，あるいは既に感染したウイルスの種類と変異株の種類の組み合わせにより予防効果が異なることは原理的にも経験的にも明らかである。 それぞれの組み合わせについて，標準的な効果との比を与えることで効果の違いを表現する。 ここでは獲得免疫の効果を，発症機序における快復開始期間の短縮に単純化し， 感染，発症，重症化それぞれに対する予防効果はその結果として観測されるものとする。
© Tatsuo Unemi, 2020, 2021. All rights reserved.
オイラー法を用いる。 ↩
2個体間の組み合わせをすべて計算せずに，あらかじめ定めた距離内にいる個体だけについて計算する。 ↩
ボックス＝ミュラー法により生成した正規分布乱数を元に，上下限と歪度を設定する関数により変更を加えた数値を用いる。 ↩
Centers for Disease Control and Prevention, USA, "Information about the Pfizer-BioNTech COVID-19 Vaccine," 2021, Jan. 25. https://www.cdc.gov/coronavirus/2019-ncov/vaccines/different-vaccines/Pfizer-BioNTech.html ↩
Pfizer and BioNTech, FDA Briefing Document Pfizer-BioNTech COVID-19 Vaccine, Vaccines and Related Biological Products Advisory Committee Meeting December 10, 2020, https://www.fda.gov/media/144245/download ↩
Sara E. Oliver et al, The Advisory Committee on Immunization Practices’ Interim Recommendation for Use of Pfizer-BioNTech COVID-19 Vaccine -- United States, December 2020, Morbidity and Mortality Weekly Report, Centers for Disease Control and Prevention, USA. https://www.cdc.gov/mmwr/volumes/69/wr/mm6950e2.htm?scid=mm6950e2w ↩
Centers for Disease Control and Prevention, USA, "Grading of Recommendations, Assessment, Development, and Evaluation (GRADE): Moderna COVID-19 Vaccine," https://www.cdc.gov/vaccines/acip/recs/grade/covid-19-moderna-vaccine.html ↩
Our World in Data: Statistics and Research - Coronavirus (COVID-19) Vaccinations, https://ourworldindata.org/covid-vaccinations?country=ISR ↩
The Wall Street Journal: Pfizer Vaccine Less Effective Against Delta Infections but Prevents Severe Illness, Israeli Data Show, July 6, 2021, https://www.wsj.com/articles/pfizers-covid-19-vaccine-is-less-effective-against-delta-variant-israeli-data-show-11625572796 ↩