Scipyによる単純な個体群成長モデルの数値解法(ロトカ-ヴォルテラ方程式など)

ODE

目次

  1. 指数成長モデル
  2. ヴェアフルストのモデル（ロジスティック方程式）
  3. ロトカとヴォルテラのモデル

はじめに

常微分方程式によって成長モデルを記述することはパラメータの解釈が容易である、という点で教育的だ。基礎的な人口増加のモデルによれば、個体数Nは微小時間 d tにその個体数そのものに比例して増加すると考えられる。

\frac{dN(t)}{dt} = rN(t)

ここで rは増加率である。上式は最も簡単な部類の微分方程式であり、解はN_0を初期個体数として次のようになる。

N(t) = N_0 e^{rt}

これは指数成長モデル、あるいはマルサスのモデルと呼ばれ*1、人口が指数的に増加することを説明する論拠の一つとなっている。

指数成長モデル

指数成長モデルが現実の結果に適用できるかどうか確かめてみよう。 以下の例はアメリカの人口動態を示している（データは『微分方程式で数学モデルを作ろう』より*2）。

# 年度  米国の人口(単位: 百万人)
1790      3.9
1800      5.3
1810      7.2
1820      9.6
1830     12.9
1840     17.1
1850     23.2
1860     31.4
1870     38.6
1880     50.2
1890     62.9
1900     76.0
1910     92.0
1920    106.5
1930    123.2

手始めに、この数値に対して1790年と1800年のデータに対してフィットするモデルを作ろう。

決めるべきパラメータは N_0と \gammaである。 まず、

N(1790) = N_0 e^{\gamma(t - 1790)}

だから、t=1790において N(1790) = N_0であることは明らかである。よってデータより N_0 = 3.9としてよい。

続いて、

N(1800) = N_0 e^{\gamma(t - 1800)}

に N_0を代入して移項し計算すれば

\gamma = \frac{\ln{(5.3/3.9)}}{(1800 - 1790)} = 0.03067...

が求まる。Numpy+matplotlibによる計算とグラフ化コードは以下のようになる。

import numpy as np
from matplotlib import pyplot as plt

def expmodel(x):
    N0 = 3.9
    g = 0.031
    Nx = N0 * np.exp(g * (x - 1790))
    return  Nx

# データは
# http://cosmo.phys.hirosaki-u.ac.jp/wiki.cgi/maxima
# の、『人口動態と微分方程式 - Maxima による数式処理』 より
data = open("usa-euc.dat").readlines()[3:]

x, y = zip(*[map(float, d.split()) for d in data])
x = np.array(x)
y = np.array(y)
m = expmodel(x)

plt.plot(x, y, ".", label=r"$actual\ data$")
plt.plot(x, m, label=r"$N(t) = N_0 \exp{(\gamma(t_0 - t))}$")
plt.legend(loc="upper left")
plt.title(r"$Exponential\ model$")
plt.xlabel(r"$time\ [year]$")
plt.ylabel(r"$population\ [million]$")
plt.show()

f:id:xef:20130202110731g:plain

グラフを見ると、近い将来はある程度正確に予想できているが、1860年以降は実際の値との「ずれ」が甚だしくなっていることが見て取れる。このずれの評価をするために (1 - (model/actual))と比率を取ってみるのも良いかもしれないが、大した誤差評価ができるわけでもないのでここでは割愛する。

フェルフルストのモデル*3

指数モデルでは(\gamma > 0で）人口は単調、単純に増加する。しかしこれは長期的な人口動態を表現するためには向かないということがここまでの観察で分かった。人口増加は続くが、様々な要因によってそれは鈍化するのである。フェルフルストは、これまでの基礎モデルに人口過密の要因を取り入れた次の式を考案した。

式は N_{max}を人口の上限として、

\frac{dN(t)}{dt} = \gamma N(t) \biggl(1 - \frac{N(t)}{N_{max}}\biggr)

である。すなわち、人口が増えるごとに1 - N/N_{max}の割合で増率が逓減するように調整したのである。このような形式の微分方程式を、今日では広くロジスティック方程式と呼んでいる。ロジスティック方程式は一階非線形微分方程式であり、解析解が存在する。

この方程式を解くにはまず n(t) = N(t)/N_{max}として次のように書き直す。

\frac{dn(t)}{dt} = \gamma\ n(t)\ (1-n(t))

これを変数分離する。

\frac{1}{\gamma\ n(t)\ (1-n(t))}dn(t) = dt

左辺を部分分数分解。

\frac{1}{\gamma}\biggl(\frac{1}{n(t)} + \frac{1}{1-n(t)}\biggr)dn(t) = dt

この積分は簡単に計算できて、

\gamma^{-1}(\log{(n(t))} - \log{(n(t)-1)}) = t + N_0

n(t)について整理して最終的にN(t)についての式を得ればよい（飽きてきた）。

N(t) = N_0\ \frac{N_{max}}{N_0 + (N_{max} - N_0)e^{\gamma(t-t_0)}}

これがフェルフルスト・モデルの解である。

この式がアメリカの人口動態に最適フィットするような N_{max}を求めよう。

Scipyには汎用の最小二乗法フィッティングソルバscipy.optimize.leastsqが存在し、これを使うだけで簡単に最適なN_{max}を発見してくれる。

# -*- coding: utf- -*-
import numpy as np
from scipy.optimize import leastsq
# 実データを読み出す
data = open("usa-euc.dat").readlines()[3:]

x, y = zip(*[map(float, d.split()) for d in data])
x = np.array(x)
y = np.array(y)

def logistic_model(param, x, y):
    N0 = 3.9
    g = 0.031
    Nmax = param[0]
    y0 = 1790
    Nx = y - N0 * (Nmax/(N0+(Nmax-N0)*np.exp(-g*(x-y0))))
    return Nx

r, info = leastsq(logistic, [100], args=(x, y))
print r[0] # => 206.82438721

こうして得られた値を元に、モデルを再計算しグラフ化したものを以下に示す。

f:id:xef:20130202160208g:plain

指数モデルと比べてかなりよくフィッティングしていますね！

ロトカとヴォルテラのモデル

これまでは1種類の生物のみの個体数の変化を調べてきた。

ここではさらに話を進めて、2種類以上の生物が存在し、相互に影響しあう場合のそれぞれの個体数の変動を考えてみたい。

もっとも3種以上の生物が相互に影響しあう場合、個体数の変動はカオス的になる場合があり単純ではない。

さしあたり2種の生物について、すなわち被食者(prey)と捕食者(predetor)の関係を考えてみよう。

被食者xと捕食者yの数の増減関係を記述する微分方程式はロトカ-ヴォルテラ方程式と呼ばれ、次のように表現される。

\frac{d x}{d t} = x(\alpha - \beta y)
\frac{dy}{dt} = -y(\gamma - \sigma x)

ここでゴテゴテとついてくる係数 \alpha,\ \beta,\ \gamma,\ \sigmaについて説明する。まず\alphaはxの増加に関わる項であり、出生率とでも呼び習わされるべきものである。 \betaは xと yの積にかかりxの減少に関わる項であり、これはxの被食率であると考えられるだろう。 \gammaはyの死亡率、 \sigmaは捕食によるyの増加率である。

解析的な話には足を踏み入れずに、計算によって変化の動態を追う。数値解析にScipyのODEソルバを使っている。

# -*- coding: utf-8 -*-
import numpy as np
from numpy import array
from matplotlib import pyplot as plt
from scipy.integrate import odeint

def lotka_volterra(ys, ts, a, b, c, d):
    prey, pred = ys
    d_prey = prey * (a - b * pred)
    d_pred = pred * ( -c + d*b * prey)
    return array([d_prey, d_pred])

def main():
    prey = 10
    pred = 5

    a = 2.      # 被食者の出生率
    b = .1      # 被食者の被食率
    c = 1.5     # 捕食者の死亡率
    d = 0.75    # 捕食者の増加率

    ys = array([prey, pred])
    ts = np.linspace(0, 10., 5000)

    soln = odeint(lotka_volterra, ys, ts, args=(a, b, c, d))
    rs, fs = soln.T

    plt.plot(ts, rs, label=r"$dx/dt = x(\alpha - \beta y)$")
    plt.plot(ts, fs, label=r"$dy/dt = -y(\gamma - \sigma x)$")

    plt.title("$Lotka-Volterra\ Equation$")
    plt.xlabel("$time$")
    plt.ylabel("$Population$")
    plt.legend()
    plt.show()

if __name__ == "__main__":
    main()

f:id:xef:20130206144612g:plain

被食者の増加が指数的であるため変化が急峻であるが、ここにロジスティック的な抑制項を加えたらどうなるだろうか。

\frac{dx}{dt} = x(\alpha - \beta y - \epsilon x) 
# 抑制項 e を追加したバージョン
def lotka_volterra_logistic(ys, ts, a, b, c, d, e):
    prey, pred = ys
    d_prey = prey * (a - b * pred - e * prey)
    d_pred = pred * ( -c + d*b * prey)
    return array([d_prey, d_pred])

mainを変更。グラフの見栄えのために他のパラメータも少しいじってある。

+     e = 0.002   # 被食者増加の抑制項

-    soln = odeint(lotka_volterra, ys, ts, args=(a, b, c, d))
+    soln = odeint(lotka_volterra_logistic, ys, ts, args=(a, b, c, d, e))

ロジスティック項によりあたかも減衰振動のような動きを取り、時間が経つと一定の値に落ち着くことが分かる。

f:id:xef:20130206145738g:plain

後半やや尻切れだが力尽きたのでここまで。

Reference & Further Reading

マルサス・モデル

フェルフルスト・モデル

Lotka-Volterra型方程式

  • ScipyによるLotka-Volterra方程式の数値解法はCookbook/LoktaVolterraTutorialをほぼそのまま拝借している。といってもこれもODE solverを使っているだけなのでここから何かしら理論的な意味を見出すことは難しい。

  • 振動する生態系— Lotka-Volterra モデル(PDF)は短いものの、一からシミュレーションプログラムを作成するために順序だった説明を与えている。

書籍

*1:もっとも、指数的な増加を示す人口動態はマルサス以前から知られていたらしい。ソースはロジスティック式 #歴史 - Wikipediaですが

*2:ただしここに貼ってあるのは人口動態と微分方程式 - Maxima による数式処理からの手抜きコピペ

*3:ヴェアフルストという発音が見受けられたけど（と言うか最初に見たのがそれだったのでそういう発音だと思い込んでいた）ドイツ語って濁音発音しないよね。ピエール＝フランソワ・フェルフルスト - Wikipedia