📂力学系

マルサス成長モデル：理想的な集団成長

モデル

$$ \dot{N} = rN $$

変数

• $N(t)$: $t$時点における集団の個体数を表す。

パラメータ

• $r \in \mathbb{R}$ : 固有成長率^{intrinsic Rate of Increase}として、$0$より大きければ成長し、$0$より小さければ衰退する。繁殖率^{birth rate}$b$と死亡率^{death rate}$d$の差$r:=b-d$として定義されることもある。

説明

個体群動態^{population dynamics}は動力学が数理生物学へとつながる最初の通路であり、集団の個体数や種の共存といった主題に対する数学的アプローチである。マルサス成長モデルはそうした成長モデルの中でも最も単純なモデルであって、理想的に(何の妨げもなく)繁殖できる種の個体数は単純な常微分方程式で表すことができる。与えられた方程式は分離可能な1階微分方程式であるから、その解は初期人口$N_{0}$に対して

$$ N(t) = N_{0} e^{rt} $$

のように求められ、式に指数関数が現れるだけあって指数的に増加するという表現がぴったり合うため、指数成長モデル^{exponential Growth model}とも呼ばれる。マルサス成長モデルという名は、人口論の著者トマス・ロバート・マルサスの名にちなんだものである。

一見、個体数モデルがどこに役立つのか疑問に思うかもしれないが、これより進んだモデルは細菌やウイルス、あるいは特定の塩基配列の数と関連して実際にバイオ分野で使われており、個体数という意味を捨てて「動く量的変数」そのものとしてアプローチすることで数多くの応用モデルの根幹となることもある。たとえばインテルの共同設立者であるゴードン・ムーアは「半導体集積回路の性能が2年ごとに倍に増加する」と述べたが、これは文字どおり指数的成長を意味し、現在はムーアの法則 と呼ばれている。ムーアがマルサス成長モデルを適用して述べたという意味ではなく、成長という現象が概してこのように説明されるがゆえに重要だということである。バイオ分野でなくとも、この成長というものは経済・経営でサービスを利用する顧客や無リスク資産の元利合計に適用されることもあり、逆に原子核の放射線崩壊を原子核数の逆成長として見ることもできる。

導出

集団の成長をモデリングするには、二分法で繁殖する細菌を想像してみるとよい。個体数が増加する速度は、それぞれの成長速度とともに、全体の個体数そのものにも比例せざるを得ないだろう。たとえば細菌が$10$匹いる培養皿$A$と$20$匹いる培養皿$B$を考えてみると、同じ時間が経って個体数が二倍に増えたとき$A$には$20$匹、$B$には$40$匹いることになる。言い換えれば

$$ (N = 10 \text{일 때의 증가량}) = 10 \\ (N = 20 \text{일 때의 증가량}) = 20 $$

であるから、文字$N$を使って表してみると

$$ (N \text{의 증가량}) = N $$

が成立すると仮定できるのである。ここで成長速度を考慮できるように式を修正すると、ある実数$r \in \mathbb{R}$に対して

$$ (N \text{의 변화량}) = rN $$

のように表せるようになる。$r>0$なら人口の変化量が正であるということなので個体数が増加し、$r<0$なら変化量が負であるので時間が経つにつれて個体数が減少するだろう。さてこの変化量というものを言葉で書かずに数式で表現するには微分が必要である。$t$時点で$N$が変化する度合いは$dN / dt$であるから、左辺を修正すると次を得る。

$$ {{dN} \over {dt}} = r N $$

■

限界

微分方程式という言葉や、数式が怖く感じられるかもしれないが、いざ一つひとつ紐解いてみるとすべて常識的な前提から出発する導出であることに共感できるだろう。しかし常識的な前提だという説明が色あせるほど、このモデルの最大の問題点は現実をまったく反映できないということである。易しく単純であるため教科書の最初の例題としては適切かもしれないが、これだけでは実質的な応用はまったくできない。

実際の集団の成長は、いわゆるマルサスの罠^{malthusian Trap}とも呼ばれる限界と向き合わざるを得ない。培養皿に入った細菌の例を引き続き考えてみると、培養皿という小さな系の中で、彼らが永遠に成長と繁殖を続けられるほどの栄養が供給されるわけでもなく、その空間自体が有限でもある。結局ある瞬間を境に成長勢いは折れ、理想的な成長は止まってしまうだろう。

このモデルは最も単純であるがゆえに実際の個体数をフィッティングする形で使われる場合はほとんどなく、それ自体が1次項としての意味だけを持つか、非現実性を克服するためのさまざまなモデルが使われる。たとえば死を反映したり、複数の種族が互いに競争する形の方法を考えてみることができるだろう。

視覚的理解

マルサス成長モデルをエージェントベースのシミュレーションで実装し、視覚的に理解してみよう。

アクション

(すべてのエージェントが、毎ターン)

• 繁殖: bの確率で自分の位置に新しいエージェントを作る。
• 死亡: dの確率でシミュレーションから除外される。

b # 번식률
d # 사망률
replicated = (rand(N) .< b) # 번식 판정
new\_coordinate = coordinate[replicated,:]
coordinate = coordinate[rand(N) .> d,:] # 사망 판정
coordinate = cat(coordinate, new\_coordinate, dims = 1);

エージェントの繁殖と死亡という単純なアクションのみを使い、システムで固有成長率$r$に該当するパラメータが変わるにつれてシステムがどのように変わるかをチェックすれば十分である。

注意すべきなのは、シミュレーションに使われる繁殖確率bとシステムの繁殖率$b$、死亡確率dとシステムの死亡率$d$は、たいてい同じではないということである。微分方程式で表現される自律システムとは正確には一致しないため、シミュレーションに対するフィッティングはそこそこ適切なパラメータを直観的に見つける形で行われた。

$r>0$の場合

N0 = 50 # 초기 인구수
b = 0.05 # 번식률
d = 0.02 # 사망률

シミュレーションであるだけに最初は少しもたつくこともあるが、繁殖率が死亡率より大きいため結局は爆発的な成長が起こる。

$r<0$の場合

N0 = 50 # 초기 인구수
b = 0.04 # 번식률
d = 0.05 # 사망률

理論的な絶滅の手順をほぼ正確にたどっていくのを見ることができる。

$r=0$の場合

N0 = 50 # 초기 인구수
b = 0.05 # 번식률
d = 0.05 # 사망률

微分方程式で表現したときには初期値が固定点となり変動がないはずだが、シミュレーションではその瞬間瞬間の運のせいで絶滅直前まで行ってから再び回復することもある。この動きは一種のブラウン運動にも見えるが、実際に人口を増やしたり繁殖率と死亡率を下げたりすれば、もう少し安定した平衡状態を維持するだろう。

• 実際には幾何ブラウン運動である。

コード

次はこの投稿に使われたJuliaコードである。

cd(@__DIR__) # 파일 저장 경로

@time using Plots
@time using Random
@time using Distributions
@time using LinearAlgebra
@time using DifferentialEquations

#---

function malthusian_growth!(du,u,p,t)
  N = u[1]
  r = p
  du[1] = dN = r*N
end

u0 = [50.0]
p = 2.65

tspan = (0.,1.8)
prob = ODEProblem(malthusian_growth!,u0,tspan,p)
sol = solve(prob; saveat=(0.0:0.1:1.8))

compare = plot(sol,vars=(0,1),
  linestyle = :dash,  color = :black,
  size = (400,300), label = "Theoretical", legend=:topleft)

#---

N0 = 50 # 초기 인구수
b = 0.05 # 번식률
d = 0.02 # 사망률
max_iteration = 180 # 시뮬레이션 기간
gaussian2 = MvNormal([0.0; 0.0], 0.03I) # 2차원 정규분포

Random.seed!(0)
time_evolution = [] # 인구수를 기록하기 위한
let
  coordinate = rand(gaussian2, N0)'
  N = N0

  anim = @animate for t = (0:max_iteration)/100
    row2 = @layout [a{0.6h}; b]
    figure = plot(size = [300,500], layout = row2)

    plot!(figure[1], coordinate[:,1], coordinate[:,2], Seriestype = :scatter,
      markercolor = RGB(1.,94/255,0.), markeralpha = 0.4, markerstrokewidth	= 0.1,
      aspect_ratio = 1, title = "t = $t",
      xaxis=true,yaxis=true,axis=nothing, legend = false)
      xlims!(figure[1], -10.,10.)
      ylims!(figure[1], -10.,10.)

    replicated = (rand(N) .< b) # 번식 판정
    new_coordinate = coordinate[replicated,:]
    coordinate = coordinate[rand(N) .> d,:] # 사망 판정
    coordinate = cat(coordinate, new_coordinate, dims = 1);

    N = size(coordinate, 1)
    push!(time_evolution, N)
    coordinate = coordinate + rand(gaussian2, N)'

    if t < 0.9
      plot!(figure[2], sol,vars=(0,1),
        linestyle = :dash,  color = :black,
        label = "Theoretical", legend=:bottomright)
    else
      plot!(figure[2], sol,vars=(0,1),
        linestyle = :dash,  color = :black,
        label = "Theoretical", legend=:topleft)
    end
    plot!(figure[2], 0.0:0.01:t, Time_evolution,
      color = RGB(1.,94/255,0.), linewidth = 2, label = "Simulation",
      yscale = :log10, yticks = 10 .^(1:4))
    ylims!(figure[2], 0.,min(time_evolution[end]*2,10000.))
  end
  gif(anim, "malthusian_growth_integration1.gif", fps = 18)
end

#---

function malthusian_growth!(du,u,p,t)
  N = u[1]
  r = p
  du[1] = dN = r*N
end

u0 = [50.0]
p = -1

tspan = (0.,1.8)
prob = ODEProblem(malthusian_growth!,u0,tspan,p)
sol = solve(prob; saveat=(0.0:0.1:1.8))

compare = plot(sol,vars=(0,1),
  linestyle = :dash,  color = :black,
  size = (400,300), label = "Theoretical", legend=:topleft)


#---

N0 = 50 # 초기 인구수
b = 0.04 # 번식률
d = 0.05 # 사망률
max_iteration = 180 # 시뮬레이션 기간
gaussian2 = MvNormal([0.0; 0.0], 0.03I) # 2차원 정규분포

Random.seed!(0)
time_evolution = [] # 인구수를 기록하기 위한
let
  coordinate = rand(gaussian2, N0)'
  N = N0

  anim = @animate for t = (0:max_iteration)/100
    row2 = @layout [a{0.6h}; b]
    figure = plot(size = [300,500], layout = row2)

    plot!(figure[1], coordinate[:,1], coordinate[:,2], Seriestype = :scatter,
      markercolor = RGB(1.,94/255,0.), markeralpha = 0.4, markerstrokewidth	= 0.1,
      aspect_ratio = 1, title = "t = $t",
      xaxis=true,yaxis=true,axis=nothing, legend = false)
      xlims!(figure[1], -10.,10.)
      ylims!(figure[1], -10.,10.)

    replicated = (rand(N) .< b) # 번식 판정
    new_coordinate = coordinate[replicated,:]
    coordinate = coordinate[rand(N) .> d,:] # 사망 판정
    coordinate = cat(coordinate, new_coordinate, dims = 1);

    N = size(coordinate, 1)
    push!(time_evolution, N)
    coordinate = coordinate + rand(gaussian2, N)'


    plot!(figure[2], sol,vars=(0,1),
      linestyle = :dash,  color = :black,
      label = "Theoretical", legend=:topright)
    plot!(figure[2], 0.0:0.01:t, Time_evolution,
      color = RGB(1.,94/255,0.), linewidth = 2, label = "Simulation",
      yscale = :log10, yticks = 10 .^(1:4))
    ylims!(figure[2], 0., 50.)
  end
  gif(anim, "malthusian_growth_integration2.gif", fps = 18)
end


#---

function malthusian_growth!(du,u,p,t)
  N = u[1]
  r = p
  du[1] = dN = r*N
end

u0 = [50.0]
p = 0

tspan = (0.,3.)
prob = ODEProblem(malthusian_growth!,u0,tspan,p)
sol = solve(prob; saveat=(0.0:0.1:3.0))

compare = plot(sol,vars=(0,1),
  linestyle = :dash,  color = :black,
  size = (400,300), label = "Theoretical", legend=:topleft)


#---

N0 = 50 # 초기 인구수
b = 0.05 # 번식률
d = 0.05 # 사망률
max_iteration = 300 # 시뮬레이션 기간
gaussian2 = MvNormal([0.0; 0.0], 0.03I) # 2차원 정규분포

Random.seed!(0)
time_evolution = [] # 인구수를 기록하기 위한
let
  coordinate = rand(gaussian2, N0)'
  N = N0

  anim = @animate for t = (0:max_iteration)/100
    row2 = @layout [a{0.6h}; b]
    figure = plot(size = [300,500], layout = row2)

    plot!(figure[1], coordinate[:,1], coordinate[:,2], Seriestype = :scatter,
      markercolor = RGB(1.,94/255,0.), markeralpha = 0.4, markerstrokewidth	= 0.1,
      aspect_ratio = 1, title = "t = $t",
      xaxis=true,yaxis=true,axis=nothing, legend = false)
      xlims!(figure[1], -10.,10.)
      ylims!(figure[1], -10.,10.)

    replicated = (rand(N) .< b) # 번식 판정
    new_coordinate = coordinate[replicated,:]
    coordinate = coordinate[rand(N) .> d,:] # 사망 판정
    coordinate = cat(coordinate, new_coordinate, dims = 1);

    N = size(coordinate, 1)
    push!(time_evolution, N)
    coordinate = coordinate + rand(gaussian2, N)'


    plot!(figure[2], sol,vars=(0,1),
      linestyle = :dash,  color = :black,
      label = "Theoretical", legend=:topleft)
    plot!(figure[2], 0.0:0.01:t, Time_evolution,
      color = RGB(1.,94/255,0.), linewidth = 2, label = "Simulation",
      yscale = :log10, yticks = 10 .^(1:4))
    ylims!(figure[2], 0., 100.)
  end
  gif(anim, "malthusian_growth_integration3.gif", fps = 18)
end