ibaibabaibaiのサイエンスブログ

サイエンス中心の予定ですが，何を書くかわかりません．統計とかの話はこっちに書くつもり． https://sites.google.com/site/iwanamidatascience/memberspages/ibayukito　 ツイッターは@ibaibabaibai

はしか物語（３）【改訂版】 当たり年・外れ年・カオス


最初に戻って，麻疹の流行を非線形力学系の視点から眺める．ここに書いたようなことはいろいろな伝染病にあてはまる可能性があるが，（１）感染力が強く，（２）感染して症状の出ない人がほとんどおらず，（３）多くの人が比較的すみやかに回復し，（４）回復後は相当の期間にわたって強い免疫が保たれる，といった性質が，麻疹を特に興味深いものにしている．

麻疹のダイナミクス再び

最初の回では患者の数の変動をあらわす確率モデルの例を考えた．伝染病が「自滅」しては再生する様子を示すにはよいが，実際の麻疹のモデルと考えると，いろいろ問題がある．

まず，ロンドンのような大きな都市では，麻疹は毎年流行するので，流行の周期は１年に固定されていることになる．大きい都市での変動は，最初の回で青い線で書いたグラフに相当するが，モデルでは流行の周期は都市の大きさ（あるいは新しく出生して入れ替わる数）によって連続的に変わるので，明らかに事実とは合っていない．

実際には，大きな都市で周期的な流行を作るのは，さまざまなパラメータの季節変動だと考えられる．特に新入学などによる「患者からまだ感染していない人への伝達効率」の季節変動が重要だとされる．あとで紹介するモデルでも，この部分に季節変動を入れている．

一方，最初の回で赤い線で書いたグラフに相当するのはもっと小さい都市である．この場合には，絶滅と再侵入によって周期的な流行が起きる．実際の麻疹でも，小さな都市での流行の周期は都市のサイズによって変化することが知られており，その意味ではいまの確率モデルはかなりいい線をいっている．バートレットなどの確率モデルの研究は主にこちらの場合を想定していると思われる．

最初の回に紹介した文献に実際の麻疹の時系列データがいくつか出ているので紹介しておく．

こちらのリンクから読める英文文献のFig.1には，London(300万人），Plymouth(21万人)，Teignmouth（1万人）の各都市での1940年代から1960年代の麻疹の流行が示されている．この3つはそれぞれ，(1) ほぼ毎年流行し非流行期にも患者数がゼロにはならない (2) ほぼ毎年流行するが非流行期には患者数がゼロになる (3) 独自の周期でやや不規則に流行する，の典型例とされる（ゼロになるかどうかは地図上のどの範囲をユニットとして考えるかにもよるように思えるが）
A stochastic model for extinction and recurrence of epidemics: estimation and inference for measles outbreaks

日本語のこれ（PDF)　の図５（ワクチンが無い時代の麻疹の流行周期）には都市や島の人口による周期の変化を示すデータが出ている．大きな都市ほど周期が短いことがわかる．
http://www.eiken.co.jp/modern_media/backnumber/pdf/MM1007_03.pdf

最初の回に紹介したモデルは，それ以外にも問題がある．実は，確率的要素のないモデル（この回の最後に出てくる微分方程式モデル）は季節変動を組み込まないと振動がしだいに減衰して安定点に収束してしまうのである．最初の回のモデルで突然の絶滅まで振動が続くのは，乱数を入れているため以外に，時間を離散化したためである可能性もある．この場合，ある種の生態学のモデルなどでは，離散時間のモデル（世代の重なりのないモデル）が現実的な場合もあるが，麻疹ではちょっと疑問である．

当たり年と外れ年

以下では，大きな都市での流行について考える．季節変動を取り入れたモデルを作ること自体は難しくないが，単に新入学とか休み明けの影響で流行が決まって，それで１年周期になっている，というのはあまり面白くないように思われる．この場合には，最初の回で見たような，麻疹自体がもつ「振動する性質」は全然出てこないのだろうか．

実は面白いことがある．毎年流行するといっても，実際のデータをよく見ると当たり年と外れ年がおおむね交互にある．つまり2年周期である．

これを論じた論文としてよく引用されるのは

Oscillatory fluctuations in the incidence of infectious disease and the impact of vaccination: time series analysis. - PubMed - NCBI

である（フリーで落とせる）．Robert Mayはカオスの研究で有名だが，この論文にはカオスは出てこない．オタフク風邪や百日咳では3年周期が見られるとも書いてある．

直観的にいうと，感染のしやすさなどの季節変動が僅かであれば，それに直接に対応して1年周期で患者数や免疫を持つ人の数が増減するだけである．

ところが，もっと季節変動が大きくなると話が違ってくる．たとえば，春先だけ感染のしやすさが増えるとすると，春だけ大量の患者が出て，免疫を持つ人の数もそこでどかっと増える．すると翌年にまで影響が及んで，翌年の春には患者がそれほど発生しない．すると，その次の年にはまたどかっと発生し・・のようになるわけだ．

外部からの駆動による周期のほかに，その整数倍の周期があらわれることがあるのは，一般に「非線形」のシステムの特徴である．おまけ１でちょっとだけ一般論を説明したので，興味のある方は読んでみてほしい．ナトリウム灯の例はちょっと面白いかもしれない．

麻疹のシミュレーション（周期外力あり/乱数なし）

上では2年周期の原因を言葉で説明したが，実際にそうなるかは，そんなに簡単に断言できない．できれば，ちゃんとした数式で2年周期が出したい．

数式としてはいろいろ考えられるが，ここでは以下の論文のものを取り上げよう（麻疹の微分方程式モデル自体は1920年代のSoperとか，さらに前のHamerとかに遡る）．

www.researchgate.net

具体的な数式とその解説はおまけ２に書いた．以下では，RのパッケージdeSolveを使って【上のモデル＋麻疹用のパラメータ】を解いた結果を見せる（コードはおまけ３を参照）．

まず，季節変動の強さをあらわすパラメータ\beta_1の値が0.05の場合．これは純粋の1年周期である．

f:id:ibaibabaibai_h:20160814220701p:plain

同じモデルで\beta_1の値を0.2に増やすと当たり年・外れ年のある2年周期に変化する．どちらも最初の過渡的な部分を捨てていることに注意．

f:id:ibaibabaibai_h:20160814220734p:plain

とりあえず，当たり年・外れ年の交代をシミュレートすることができた．

これまであげた文献のほか，LondonとYorkの1973年の論文がデータ解析とシミュレーションの両方でよく知られているらしい．読んでいないが，メモとしてリンクしておく（PDF)．
http://yorke.umd.edu/Yorke_papers_most_cited_and_post2000/1973_03_London_I_AMJEpi.pdf
http://yorke.umd.edu/Yorke_papers_most_cited_and_post2000/1973_04_London_II_AMJEpi.pdf
上のほうにあるデータだとオタフク風邪や水疱瘡では，１年以上の周期ははっきりしないようだ．

それ以外の周期も出る

ところが，このモデルは，2年周期だけでなく，他にもいろいろな挙動を示す．

たとえば\beta_1=0.1とすると，1年の整数倍でないような振動が見られる．実はこれは文献に載ってないので，計算ミスと思ったのだが，たぶん正しいようだ（読者が再現できたらツイッターなどで教えてほしい）．
f:id:ibaibabaibai_h:20160814221028p:plain

また，季節変動をあらわすパラメータを大きくして行くと，4年周期，8年周期，16年周期・・も出現する．たとえば，8年周期の例はこうなる（\beta_1=0.268）．
f:id:ibaibabaibai_h:20160814221240p:plain
2年周期がベースにあって，丁寧にみると8年周期になっているという感じだ．

カオスも出る

だんだん長い周期が出てきて，最後はどうなってしまうのか？　実はある限界があって，その先は周期解でなく不規則な変化をするようになる．いわゆるカオス解の出現である．

「カオス」については，名前を聞いたことがある人も多いと思う．たまに柑橘類の「カボス」と間違っている人もいるが，そっちではなくて，初期値の僅かな差や微妙な雑音，数値的な丸め誤差で結果が大きくかわる（バタフライ効果）というアレである．乱数を入れなくても，初期値に微妙に依存して，いつまでもランダムに動き続けるのである．

カオスというと文字通り混沌とした状態を想像する人もいると思うが，そういうものばかりではない．メインの変動はほぼ2年周期であって，よく見るとその上に不規則な動きが乗っかっているという場合もある（\beta_1=0.27）．見かけ上，さっきの8周期とほとんど変わらないが，こんどはいくら長くやっていても周期がない．

f:id:ibaibabaibai_h:20160814221310p:plain

一方で，これはかなりデタラメ感がある（\beta_1=0.272）．

f:id:ibaibabaibai_h:20160814223025p:plain

もっと長くプロットすると

f:id:ibaibabaibai_h:20160814230900p:plain

ところが，これには問題があって，Matlabで計算したり，方程式を解くときのやり方（刻み幅）を変えると，結果が変わってしまったりする．それで，最初はプログラムが間違っていると思った．実は解き方だけでなく，\beta_1の値をホンの僅か変えても違う結果になるので，これは普通でいう間違いというよりはカオス系特有の不安定な振る舞いなのだろう．

基本的には，（１）2年周期がほぼ保たれて僅かに乱れるカオス，（２）2年周期が完全に壊れたカオス，あるいは壊れた状態と2年周期がほぼ保たれた状態を間欠的に移動する，のふたつがあって，ちょっとした計算法やパラメータの違いで入れ替わるように見える．ただし（２）がホンモノであることは十分に確認できていない．　

麻疹の時系列はカオスか？

ここで，はっと思いついて実際のデータを見ると，いままで雑音というか外部的なランダムネスのように思えた部分が，（ほぼ2年周期の解が少し不規則になった）カオスなのではないかと思えてくる．

「見た目の通りの大きさの雑音が入っている」と思っていた部分は，「微小な雑音がシステムの性質によって拡大して生じた」ものなのかもしれないのだ．

80年代にさまざまな分野の研究者を襲った「カオス病」の発症である．しかし，実際のデータがカオス解なのか，単に安定な周期解（リミット・サイクル）に雑音が乗ったものを見ているのか，見分けるのは簡単ではない．

実際の麻疹と水疱瘡の時系列データをいろいろな技術を駆使して調べた結果が下の論文にあるが，これによると，水疱瘡の時系列は周期解＋外部雑音だが，麻疹はおそらくカオスである可能性が高いとのこと．

カオスだと実際的に何が違うのかは，本当のところよくわからないのだが，麻疹マニアとしては，何となくうれしい気がするのである．

http://www.ganino.com/games/Science/science%20magazine%201989-1990/root/data/Science%201989-1990/pdf/1990_v249_n4968/p4968_0499.pdf

(おまけ１）線形と非線形

大まかにいうと，関数fが線形というのは，入力がxのときの出力がf(x)，yのときの出力がf(y)のとき，x+yに対する出力f(x+y)が単純な和f(x)+f(y)になることを意味する．もしx, yがただの数字ならこれは正比例とか直線関係ということになる（もともと「線形」の「線」は直線の意味である）．実際には入力x, yや出力はもっと複雑にベクトルだったり，「足し算が定義される何か」だったりする．関数でなく方程式についての「線形」は「２つの解の和がまた解になる」ことである．

一般に，システムに周期的な変化を外から加えたとき，何がおきうるかは「線形」のシステムと「非線形」のシステムで違う．

線形のシステムにある振動数で外からの力を入力すると，そのままの振動数で振動する．「振動数」の代わりにその逆数の「周期」といっても同じである．ただし，外力に反応する強さは，そのシステムの特性による．「システムをちょっと引っぱたいてあとは自然に振動させたときの振動数」（固有振動数）に外力の振動数が近いほど，同じ大きさの外力に対する反応は大きい．これを共鳴とか共振といって，たとえば地震や風による建造物の破壊とか，テレビやラジオのチャンネルの選択とか，いろんなことに関係がある．

正確にいうと，最初のうちは固有振動数の振動が外力による振動に重なってみえることがあるが，「摩擦」のようなものが少しでもあると，固有振動に相当する部分は，外からエネルギーが供給されないのでしだいに消えてしまい，外力の周期の部分だけが残る．また「外力の整数倍の周期があらたにできる」ような現象は起きない．

ここで「共鳴」は知っていても「いつも外力と同じ周期で振動する」というのは改めて意識していない人がいるかも知れない．このわかりやすい例としては「高速道路などで使われる黄色いナトリウム灯で照らされた光景」がある．

前に書いたが，光の色は光の振動数に関係しており，物体の色は物体の中身がそれを照らす光に「共鳴」することで反射・吸収が起きて生じる．ナトリウム灯はほぼ決まった振動数なので，照らされた物体の中身もその振動数で揺れ，その振動数の光を出す．だから，どんな物体も黄色くみえる．さらによく見ると，普段黄色く見えるものは明るく，違う色のものは暗く見えるが，これは共鳴の効果である．

これに対して「非線形」の場合には，入れた外力の2倍，3倍・・の周期の振動が起きることがある．上で述べた麻疹の当たり年，外れ年はまさにその例で，もとの1年周期の上にその2倍の2年周期の成分が（位相の山と谷が合うように）かぶっているのである．

（だいぶ前のバージョンで入っていた違う色の光が出てくる話は「２倍の振動数」（半分の周期）で逆でした．これも非線形性のあらわれではあるのですが，ミスリードの可能性があるので削除しました）

（おまけ２）微分方程式の説明

上の計算に用いたモデルは下記のような微分方程式モデルである．これはいわゆるSEIRモデルの仲間（永久的な免疫を仮定しているので4つの変数のうちのRは除いてある）だが，麻疹の場合にはHamerとかSoperという人たちが早い時期に研究している（SEIRモデルの初出とされる論文が1927, Soperが1929, Hamerの研究はさらに前）．

Compartmental models in epidemiology - Wikipedia, the free encyclopedia　
「感染症と文明」のSEIRモデル : wrong, rogue and log

以下で，変数は S（未感染の人の数；Suseptible），E(暴露して潜伏期にある人の数；Exposed），I（発症して人にうつせる人の数；Infectious）．

\displaystyle\frac{dS}{dt}=\mu-\beta(t)S(t)I(t)-\mu S(t)

\displaystyle\frac{dE}{dt}=\beta(t)S(t)I(t)-(\mu+\alpha)E(t)

\displaystyle\frac{dI}{dt}=\alpha E(t)-(\mu+\gamma)I(t)

ここで，季節変動を係数\beta(t)の部分に入れる．新学期は複数あるし，いきなり階段状に変わりそうなのに，コサインかよ～というツッコミは自分もしたい．

\beta(t)=\beta_0 (1+\beta_1\cos(2\pi t))

式の意味は，式ごとよりも項ごとに見て行くのがわかりやすいと思う．

  • 感染：免疫なし→潜伏状態
    • 第1式と第2式 \pm\beta(t)S(t)I(t)
    • 上で\betaの季節変動が　\beta(t)=\beta_0 (1+\beta_1\cos(2\pi t)) で与えられる．tの単位は年．
  • 発症：潜伏状態→病気の状態
    • 第2式 -\alpha E(t)
  • 治癒：病気の状態→免疫ありの状態
    • 第3式 -\gamma I(t)
  • 出生：免疫なしの人が新しく出てくる
    • 第1式 \mu
  • 死亡：出生の分だけ除去される（病気でも死亡率は増えないと仮定）
    • 第1式-\mu S(t)，第2式 -\mu E(t)，第3式 -\mu I(t)　

（おまけ３）

改稿にあたって種々ご助言頂いた @ghost_orange様，@KuboBook様，およびS氏に感謝します．

上をR言語のパッケージdeSolveで解くコードを掲載しておくが，以下のリンクの上のほう（およびそのリンク）の内容をそっくり使わせて頂いて，式を入れ替えただけである．

Rパッケージ deSolveを使う
http://d.hatena.ne.jp/teramonagi/20150531/1433042797

場合によっては，数百年以上も最初の部分を捨てないと論文のアトラクタに収束しないようだが，そのあたりの現実の麻疹との対応がイマイチよくわからない．

library(deSolve)

#方程式とパラメータを設定
parameters=c(b0=1800,b1=0.2,m=0.02,a=35.84,g=100)
SEI=function(t, state, parameters) {
	with(as.list(c(state, parameters)), {
                b=b0*(1+b1*cos(2*pi*t))
		dS =m*(1-S)-b*S*I
		dE =b*S*I-(m+a)*E
		dI =a*E-(m+g)*I
		list(c(dS,dE,dI))
	})
}

#初期値と刻みを設定
infected=0.0001
initial=c(S=1-infected, E=0, I=infected)
times=seq(0, 1200, 0.005)

#解く
out <- ode(y = initial, times = times, func = SEI, parms = parameters)

#プロットする
olen=length(out[,4])
plen=olen*0.02
col=3

windows()
plot(log(out[(olen-plen):olen,2]),log(out[(olen-plen):olen,4]),col=col,pch=".",xlab="S",ylab="I",cex=2)

windows()
plot(out[(olen-plen):olen,1],(out[(olen-plen):olen,4]),col=col,pch=".",xlab="year",ylab="I",cex=2)

（メモ）空間パターン

化学反応や非線形素子のネットワークなどで，渦巻きや進行波などのパターンが生じることが知られている，友人が教えてくれた下記の論文では，提案の解析手法を使うと麻疹などの流行に空間パターン（進行波？）が見いだせるとのこと（詳細未見，友人からの情報）．

http://www.nature.com/nature/journal/v414/n6865/full/414716a.html