感染症モデルに基づくゾンビ大量発生時の危機管理指針

新型インフルエンザが全世界を席巻しているが、エマージングウイルス（emerging virus：新興感染症）の感染爆発（outbreak: アウトブレイク）の危険性は、熱帯雨林の開発に伴う野生生物との接触拡大、交通機関の発達による感染経路の拡大により日に日に増大している。このあたりの話は山内一也東京大学名誉教授による連続講座人獣共通感染症が大変充実しており興味深い。

エボラ出血熱などの極めて高い致死率を有する感染症が全世界に広がらず人類が滅亡しないのは、その高い致死率にある。他人に感染させる前に宿主が死亡してしまうので、感染が局所的範囲に抑えられるのだ。万が一、発病までに長い潜伏期間をもつ致命的な感染症が流行すると人類の存続が危うくなる。新型インフルエンザがそうならない保証はない。

ゾンビ感染症

ヴードゥーの司祭により生み出されるゾンビ*1は、多くの小説や映画の題材となっているが、おおよその傾向として、ゾンビにやられた人間がゾンビ化する感染症としての性質を持っている。普通の感染症では宿主が死ねば他者への感染は止まるが、ゾンビ感染症の場合は、死亡した感染者が動き回り感染を広げるところが厄介である。ゾンビを止めるには首をはねるか、頭を潰すかするしかない。

仮に上記のような性質をもつゾンビ感染症が現実のものとなった場合、感染者数の推移はどのようになるのか、ゾンビ感染症を制圧にするにはどのようにすればよいのかを、数学モデルを構築して考察した論文が、"When Zombies Attack!: Mathematical Modelling of an Outbreak of Zombie"*2である*3。

人類の叡智はこの厄介なゾンビ感染症に打ち勝つことができるのか、早速基本モデルから紹介したい。

基本モデル（SZRモデル）

基本モデルでは次の3つのクラスから構成される。

S: 生者
Z: ゾンビ
R: 死者

出生率Πで生まれる生者はパラメタδでゾンビとは関係がない要因で死亡し、死者となる。死者はパラメタζによりゾンビとして生まれ変わる。生者は感染パラメタβによってゾンビと接触し新たなゾンビになる。ゾンビは生者によって首をはねるか頭を潰すかしてパラメタαによって殺され死者となる（ゾンビは他のゾンビを襲わない）。このモデルは次の数式で表される。

　　 $S' = \rm{\Pi} - \beta SZ-\delta S$
　　 $Z' = \beta SZ+\zeta R -\alpha SZ$
　　 $R' = \delta S + \alpha SZ - \zeta R$

これは感染症をモデル化するために通常利用されるSIRモデル*4よりも若干複雑になっている。ここで1匹のゾンビから単位時間あたりに感染して生まれる新たなゾンビはβSZ、一人の生者が単位時間あたりに殺すゾンビの数はαSZとなる。

この常微分方程式は
　　 $S'+Z'+R'=\rm{\Pi}$
を満たすため、
　　 $S+Z+R=\infty$
となる。これは、生者とゾンビと死者の数の合計が無限大に発散するということを意味するが、t→∞としても生者の数は無限大にはなり得ないので（人はやがて死ぬ）、ゾンビあるいは死者の数が爆発的に増えることになる。

ゾンビ感染症のアウトブレイクが極めて短期間に発生したと仮定すると、出生率およびゾンビと関係がない死亡率を無視することができる。そこで、Π=δ=0を代入し、常微分方程式の零解を求める。

　　 $- \beta SZ = 0$
　　 $\beta SZ+\zeta R -\alpha SZ = 0$
　　 $\alpha SZ - \zeta R = 0$

最初の数式からSないしはZが0になることがわかる。それに従って何れの場合でもRは0となる。S=0とすると、生者が絶滅し、"終末の日"の解を得る。

　　 $(\overline{S},\, \overline{Z},\, \overline{R}) = (0,\, \overline{Z},\, 0)$

一方、Z=0とすると、ゾンビに悩まされることのない"平和な日々"の解を得る。

　　 $(\overline{S},\, \overline{Z},\, \overline{R}) = (N,\, 0,\, 0)$

つまり、このモデルでは解の安定性に関わらず両者が共存することはないことが導き出される。相手を絶滅させるしか生き延びる道はないのだ。

オイラー法に基づいて微分方程式を解いた結果を次に示す。文末にMATLABコードを掲載しているので、実行環境のある人は試してみると良いだろう。

左のグラフはゾンビが発生しない状況を表しているが、残念ながらこの解はラウス・フルビッツの安定判別法より不安定であることが分かる（詳細割愛）。一方右のグラフは、生者を0とした場合のゾンビ感染症アウトブレイク発生時のシナリオを示している（α=0.005, β=0.0095, ζ=0.0001, δ=0.0001）。生者は早々に絶滅させられ、ゾンビが代わりに感染を広げていく。こちらの平衡点は安定しており、短期間のアウトブレイクではゾンビ感染は易々と広がることを示している。

潜伏感染を考慮したモデル（SIZRモデル）

既にエントリが肥大化してきたので、以降はモデルの解説に限定して紹介する。数式等詳細は論文を参照されたい。

さて、ゾンビ感染症が恐ろしいのは、発病までに潜伏期間（"The Zombie Survival Guide"によれば、おおよそ24時間）があり、感染者がしばらく動き回って他者に感染を広げる点にある。基本モデルではこの観点が考慮されていなかったので、新たに感染者クラス(I)を追加する。生者はまず感染者になって、しばらくの期間そのクラスに留まる。感染者はそのまま自然死をすることもあるが、さもなければゾンビとして生まれ変わる。このSIZRモデルは次のように図示される。

常微分方程式の零解を求めると、このモデルにおいても基本モデルと同様に、生者が絶滅し"終末の日"が訪れるか、ゾンビに悩まされることのない"平和な日々"の2つに1つしかないことが分かる。"終末の日"の平衡解が安定で、"平和な日々"が不安定なのも同様だ。

次にSIZRモデルにおける"終末の日"シナリオの推移グラフを示す。各パラメタは先の例と同一でやはりゾンビが地上を席巻するが、最終的なゾンビの数はおおよそ2倍にまで増加している。

隔離対策を行った場合（SIZRQモデル）

人類としても座して滅びを待つわけにはいかない。初期の新型インフルエンザへの対応にも見られたように、感染症への対策のひとつに強制隔離がある。感染者を隔離し感染が広がらないように被害を局所的に封じ込めるのだ。そこで隔離者クラス(Q)を新たに追加する。

隔離エリアには感染者及びゾンビの一部がそれぞれパラメタκ, σに従って隔離される。中には隔離エリアからの逃亡を企てるものも出てくると考えられるが、それらは成功前に銃殺される（パラメタγ）。殺された人々は死者クラス(R)に移動するが、一部はゾンビとして復活する。このSIZRQモデルは次のように図示される。

短期間のアウトブレイク(Π=δ=0)を想定すると、次の2つの平衡解を得る。
　　 $(\overline{S},\,\overline{I},\,\overline{Z},\,\overline{R},\,\overline{Q}) = (N,\,0,\,0,\,0,\,0),\,(0,\,0,\,\overline{Z},\,\overline{R},\,\overline{Q})$

やはり食うか食われるかの2つに1つだが、人類が勝利する解が安定するためには次に導かれる感染率R₀ (一次感染者から二次感染する患者数の平均値) が1未満である必要がある。ちなみに28日に発表されたCDCの報告によれば、新型インフルエンザの感染率は1.5と見積もられている*5。

　　 $R_{0} = \frac{\beta N \rho}{(\rho+\kappa)(\alpha N+\sigma)}$

感染率R₀を下げるためには、感染者やゾンビを隔離する割合を表すκかσを大きくする必要がある。人口が大きい場合R₀の近似値は次のようになる。

　　 $R_{0} \approx \frac{\beta \rho}{(\rho+\kappa)\alpha}$

β>α（人がゾンビを殺すよりも、ゾンビ感染症が人に感染する方が速い）の場合、ゾンビを撲滅できるかどうかは感染拡大の初期段階における隔離にかかっている。しかし実際のところ、感染者を見分けることが困難である場合には、徹底した隔離をすることは困難だろう。大量の感染者を隔離するための施設面を考えても、大規模隔離は非現実的である。そのため、結局κかσは小さくならざるを得ず、R₀は1よりも大きくなる。やはり今までと同様、人類が勝利する解は安定しない。

次にSIZRQモデルにおける"終末の日"シナリオの推移グラフを示す。隔離対策によって若干人類絶滅の日が後ろに延びるが結果は同じだ。

治療法が確立した場合（SIZR with cureモデル）

今まで苦杯をなめてきた人類だが、人類の叡智を結集させればゾンビ感染症を治療するワクチンなど、早期に治療法を確立できるかもしれない。どのようにゾンビ化したかに関わらずゾンビ感染症を治療し、再び生者とする治療法が確立したと仮定する。ただしこの治療法は免疫を与えるものではなく、一旦ゾンビ感染症から回復しても再び罹患する可能性は残るとする。治療法の確立をもって非実用的な隔離をやめると、モデルは次のように図示される。

短期間のアウトブレイク(Π=δ=0)を想定し、Z=0とすると次の2つの平衡解を得る。

　　 $(\overline{S},\,\overline{I},\,\overline{Z},\,\overline{R}) = (N,\,0,\,0,\,0)$
　　 $(\overline{S},\,\overline{I},\,\overline{Z},\,\overline{R}) = (\frac{c}{\beta} ,\,\frac{c}{\rho}\overline{Z},\,\overline{Z},\,\frac{\alpha c}{\zeta \beta}\overline{Z})$

この2つ目の人類とゾンビが共存する解は安定的である。次に治療法が確立した場合に人類が生き残るシナリオの推移グラフを示す。確かに人類は一定数生存することができるが、自らの何倍もの数のゾンビが跳梁跋扈する世界で、ゾンビに怯えながら生活することを余儀なくされそうだ。

攻性掃討作戦

最後に攻性対応が実現する場合を想定する。人類は準備ができた段階で戦略的にゾンビの掃討作戦を決行する。そのようなリソースの準備や調整は困難を極めることから、作戦は複数回に分けて継続的に遂行される。

ここで再び基本モデルに戻り、掃討作戦の効果を加える。

　　 $S' = \rm{\Pi} - \beta SZ-\delta S \hspace{50} t \neq t_{n}$
　　 $Z' = \beta SZ+\zeta R -\alpha SZ \hspace{25} t \neq t_{n}$
　　 $R' = \delta S + \alpha SZ - \zeta R \hspace{38}t \neq t_{n}$
　　 $\Delta Z = knZ \hspace{101} t = t_{n}$

ここで、kは殺傷率であり、nは根絶までに繰り返される掃討作戦の回数を表す。ゾンビ数の推移は次のグラフに示される（縦軸のオーダーは他のグラフの1/1000であることに注意）。

α=0.0075, β=0.0055, ζ=0.09, δ=0.0001, k=0.25としている。2.5日後に25%のゾンビが破壊され、5日後に50%、7.5日後に75%、そして10日後に100%のゾンビが破壊される。

まとめ

ゾンビ感染症のアウトブレイクは、アンデッドに対して攻性の作戦行動が伴わない限り、絶望的な状況へと陥る可能性が高い。積極的な隔離対策は感染症を撲滅できるかもしれないが、現実的に取り得るオプションとは考えにくい。治療法の確立は一部の人類の存続を可能とするが、地球はゾンビの手に落ちる。早い段階から大規模な掃討作戦を戦略的・継続的に繰り返すことだけが、ゾンビの撲滅を可能とする。

これらの結果は短期間のアウトブレイクを想定しているが、長期にわたると結果は破滅的になる。文明は破壊され、すべての人間は感染するか死亡することになる。生まれる人間も死ぬ人間も、どちらも新たなゾンビの供給源となるからだ。人類存亡の危機に際し、我々が行わなければならないことは、強い覚悟をもって絶滅させられる前に相手を絶滅させることだ。

ゾンビの大量発生が現実に起こるとは考えにくいが、対応において重要なのは、危機管理における基本である、迅速な状況把握、早期の政治決断、戦略的な作戦行動であることは現実も変わらないようだ。

Appendix

MATLABコードを以下に示す。

function [ ] = zombies(a,b,ze,d,T,dt)
% This function will solve the system of ODE’s for the basic model used in
% the Zombie Dynamics project for MAT 5187. It will then plot the curve of
% the zombie population based on time.
% Function Inputs: a - alpha value in model: "zombie destruction" rate
% b - beta value in model: "new zombie" rate
% ze - zeta value in model: zombie resurrection rate
% d - delta value in model: background death rate
% T - Stopping time
% dt - time step for numerical solutions
% Created by Philip Munz, November 12, 2008
%Initial set up of solution vectors and an initial condition
N = 500; %N is the population
n = T/dt;
t = zeros(1,n+1);
s = zeros(1,n+1);
z = zeros(1,n+1);
r = zeros(1,n+1);
s(1) = N;
z(1) = 0;
r(1) = 0;
t = 0:dt:T;
% Define the ODE’s of the model and solve numerically by Euler’s method:
for i = 1:n
s(i+1) = s(i) + dt*(-b*s(i)*z(i)); %here we assume birth rate
= background deathrate, so only term is -b term
z(i+1) = z(i) + dt*(b*s(i)*z(i) -a*s(i)*z(i) +ze*r(i));
r(i+1) = r(i) + dt*(a*s(i)*z(i) +d*s(i) - ze*r(i));
if s(i)<0 || s(i) >N
break
end
if z(i) > N || z(i) < 0
break
end
if r(i) <0 || r(i) >N
break
end
end
hold on
plot(t,s,’b’);
plot(t,z,’r’);
legend(’Suscepties’,’Zombies’)
------------
function [z] = eradode(a,b,ze,d,Ti,dt,s1,z1,r1)
% This function will take as inputs, the initial value of the 3 classes.
% It will then apply Eulers method to the problem and churn out a vector of
% solutions over a predetermined period of time (the other input).
% Function Inputs: s1, z1, r1 - initial value of each ODE, either the
% actual initial value or the value after the
% impulse.
% Ti - Amount of time between inpulses and dt is time step
% Created by Philip Munz, November 21, 2008
k = Ti/dt;
%s = zeros(1,n+1);
%z = zeros(1,n+1);
%r = zeros(1,n+1);
%t = 0:dt:Ti;
s(1) = s1;
z(1) = z1;
r(1) = r1;
for i=1:k
s(i+1) = s(i) + dt*(-b*s(i)*z(i)); %here we assume birth rate
= background deathrate, so only term is -b term
z(i+1) = z(i) + dt*(b*s(i)*z(i) -a*s(i)*z(i) +ze*r(i));
r(i+1) = r(i) + dt*(a*s(i)*z(i) +d*s(i) - ze*r(i));
end
%plot(t,z)
------------
function [] = erad(a,b,ze,d,k,T,dt)
% This is the main function in our numerical impulse analysis, used in
% conjunction with eradode.m, which will simulate the eradication of
% zombies. The impulses represent a coordinated attack against zombiekind
% at specified times.
% Function Inputs: a - alpha value in model: "zombie destruction" rate
% b - beta value in model: "new zombie" rate
% ze - zeta value in model: zombie resurrection rate
% d - delta value in model: background death rate
% k - "kill" rate, used in the impulse
% T - Stopping time
% dt - time step for numerical solutions
% Created by Philip Munz, November 21, 2008
N = 1000;
Ti = T/4; %We plan to break the solution into 4 parts with 4 impulses
n = Ti/dt;
m = T/dt;
s = zeros(1,n+1);
z = zeros(1,n+1);
r = zeros(1,n+1);
sol = zeros(1,m+1); %The solution vector for all zombie impulses and such
t = zeros(1,m+1);
s1 = N;
z1 = 0;
r1 = 0;
%i=0; %i is the intensity factor for the current impulse
%for j=1:n:T/dt
% i = i+1;
% t(j:j+n) = Ti*(i-1):dt:i*Ti;
% sol(j:j+n) = eradode(a,b,ze,d,Ti,dt,s1,z1,r1);
% sol(j+n) = sol(j+n)-i*k*sol(j+n);
% s1 = N-sol(j+n);
% z1 = sol(j+n+1);
% r1 = 0;
%end
sol1 = eradode(a,b,ze,d,Ti,dt,s1,z1,r1);
sol1(n+1) = sol1(n+1)-1*k*sol1(n+1); %347.7975;
s1 = N-sol1(n+1);
z1 = sol1(n+1);
r1 = 0;
sol2 = eradode(a,b,ze,d,Ti,dt,s1,z1,r1);
sol2(n+1) = sol2(n+1)-2*k*sol2(n+1);
s1 = N-sol2(n+1);
z1 = sol2(n+1);
r1 = 0;
sol3 = eradode(a,b,ze,d,Ti,dt,s1,z1,r1);
sol3(n+1) = sol3(n+1)-3*k*sol3(n+1);
s1 = N-sol3(n+1);
z1 = sol3(n+1);
r1 = 0;
sol4 = eradode(a,b,ze,d,Ti,dt,s1,z1,r1);
sol4(n+1) = sol4(n+1)-4*k*sol4(n+1);
s1 = N-sol4(n+1);
z1 = sol4(n+1);
r1 = 0;
sol=[sol1(1:n),sol2(1:n),sol3(1:n),sol4];
t = 0:dt:T;
t1 = 0:dt:Ti;
t2 = Ti:dt:2*Ti;
t3 = 2*Ti:dt:3*Ti;
t4 = 3*Ti:dt:4*Ti;
%plot(t,sol)
hold on
plot(t1(1:n),sol1(1:n),’k’)
plot(t2(1:n),sol2(1:n),’k’)
plot(t3(1:n),sol3(1:n),’k’)
plot(t4,sol4,’k’)
hold off

*1:写真はJoel Friesenによる。Wikipediaからの引用。

*2:J.M. Tchuenche and C. Chiyaka: When Zombies Attack!: Mathematical Modelling of an Outbreak of Zombie, Infectious Disease Modelling Research Progress, Nova Science Publishers, pp. 133-150 (2009).

*3:論文の序章ではゾンビの起源から古今東西のゾンビ小説・映画の紹介が続く。本論文ほど多くの小説や映画を参考文献としてあげた論文は久し振りに読んだ。本論文は痛いニュースでゾンビがいたらどうなるの？→カナダの科学者「世界が終わります」として取り上げられ話題になっていたが、誰も元論文を読んでいないっぽい。英語圏なら相当数が元論文を参照すると期待されるのに、言語の壁はやはり大きいようだ（せっかく面白いのに……）。本論文で取り上げられているゾンビ感染症モデルは、ゾンビが再び蘇るという点で他の感染症モデルと決定的に異なり、現実的に発生しそうもないが、政党支持者のモデル化や潜伏感染を有する疾病のモデル化等への応用は可能かもしれないと筆者らは述べている。これは最初のゾンビ感染症のモデル化であり利用される数学モデルの柔軟性を示したと言えるだろう。

*4:Brauer, F. Compartmental Models in Epidemiology. In: Brauer, F., van den Driessche, P., Wu, J. (eds). Mathematical Epidemiology. Springer Berlin 2008.

*5:新型インフルエンザの感染率の1.5という数字は今のところそれほど高くない（季節性インフルエンザの感染率が1.3程度）。1918年のスペインかぜのパンデミック時には感染率は2.5程度だったとされる。感染者6億人、死者4000〜5000万人を記録した大災害である。