📂動力学

ロジスティック成長モデル：集団成長の限界

モデル

$$ \dot{N} = {{ r } \over { K }} N ( K - N) $$

変数

• $N(t)$：$t$時点での集団の個体数を表す。

パラメータ

• $r \in \mathbb{R}$：固有増加率^{intrinsic Rate of Increase}として、$0$より大きければ成長し、$0$より小さければ衰退する。出生率^{birth rate} $b$ と死亡率^{death rate} $d$ の差 $r:=b-d$ としても定義される。
• $K$：環境容量^{carrying Capacity}として、集団を収容できる環境の大きさを描写する。

説明

この微分方程式は、ロジスティック方程式またはベルフルスト方程式^{verhulst’s equation}とも呼ばれる。

ロジスティック成長モデルは、マルサスの成長モデルで個体数が現実にはありえないほど増加するだけのことを補うために考案された。個体数が多いほど人口成長にペナルティを与えることで、いわば方程式にマルサスの罠^{malthusian Trap}を添加するようなものだ。個体数 $N$ が環境容量 $K$ を超えると $(K-N) < 0$ になるため、個体数の変化量 $n '$ がマイナスになり、人口がむしろ減少することになる。もちろん、このアイデアが単に思いつきで導入されたわけではなく、どのようにペナルティを与えるかの導出過程を直接見てみよう。

導出

マルサスの成長モデル: $$ \dot{N} = r N $$

このような単純なモデルでは、個々の個体の平均的な成長率は、両辺を $N$ で割ることで得られる。指数成長モデルでは、集団の大きさに関係なく成長率が一定なので、次のように一定となる。

$$ {{ \dot{N} } \over { N }} = r $$

ここで右辺を修正しよう。個体数が多いほど成長率が減少するようにしたいが、どのように減少するかはすぐには明確な答えがないので、最もシンプルに直線的に減少すると仮定する。つまり、個々の成長率が定数関数ではなく直線を描くことになる。

図に示されているように、右辺を修正すると、$N$ 軸の切片が $K$ であり、$N’/N$ 軸の切片が $r$ の直線の方程式になるようにすると

$$ {{ \dot{N} } \over { N }} = {{ r } \over { K }} ( K - N) $$

両辺に $N$ を掛けると次の結果が得られる。

$$ \dot{N} = {{ r } \over { K }} N ( K - N) $$

■

固定点

ロジスティック成長モデルは、二つの固定点 $N=0, N=K$ を持ち、特に $N=K$ はリャプノフ安定である¹。

視覚的理解

ロジスティック成長モデルを格子ベースのシミュレーションで実装し、視覚的に理解してみよう。

アクション

（全てのセルで、毎ターン毎に）

• 計算：各セルを基準に上下左右にどれだけの成分があるか数える。$i$行 $j$列で計算された数を$n_{ij}$とする。
• 拡散：各セルは、$1-(1-\beta)^{n_{ij}}$の確率で成分を含むようになる。

格子空間での拡散という単純なアクションのみを使用しても、格子空間自体が有限であるため、自然に成長に制限がかかることがわかる。

次に、理論的な成長推移と比較したものである。

見ての通り、シミュレーションと理論が正確に一致はしない。曲線を描くというよりは、直線的に成長し、初期の成長があまりに急である。しかし、いずれにしても環境容量を満たしながら折れ曲がる部分に注意することが重要だ。これを実生活に当てはめれば、ある培養皿で細菌が増殖していくことに例えられる。どんなに成長を続けたいと思っても、培養皿から溢れるほど無限に成長することはできない。

理論とシミュレーションが正確に一致しない理由は、このシミュレーションが十分に洗練されておらず、モデルもあまりにも単純だからである。あるいは、単にシミュレーションで運が悪かっただけかもしれないが、格子が与える影響を絶対に無視できない。ただし、この場合は単に運のせいではない。

実際に理論的な推移を比較してみると、このポストで紹介された視覚化は、むしろゴンペルツ成長モデルに近いと思われる。実際に格子空間で腫瘍が成長すると考えると、体積に対する表面積の比が小さく見え、細胞自体は増えているが、実際に腫瘍の大きさを増やせる細胞が徐々に減っていく現象が起こっている。

コード

以下は、GIFを作成するためのJuliaのコードである。

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

@time using Plots
@time using Random
@time using DifferentialEquations

row_size = 20
column_size = 20
β = 0.1 # 번식률
γ = 0.1 # 사망률

#---
K = row_size*column_size

function logistic_growth!(du,u,p,t)
  N = u[1]
  r = p
  du[1] = dN = r/K*N*(K-N)
end

u0 = [1.0]
p = 0.8

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

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


#---
max_iteration = 180

Random.seed!(0);
time_evolution = Int64[]

stage_lattice = zeros(Int64, row_size, column_size)
stage_lattice[rand(1:row_size), rand(1:column_size)] = 1

let colormap_SI = cgrad([colorant"#EEEEEE", colorant"#111111"])
  anim1 = @animate for t = (0:max_iteration)/10
    heatmap(reverse(stage_lattice,dims=1), color=colormap_SI,
      xaxis=false,yaxis=false,axis=nothing, size = [400,400], legend = false)

    I_lattice = (stage_lattice .== 1)
    count_lattice =
      vcat(I_lattice[2:end, :], zeros(Int64, 1, column_size)) +
      vcat(zeros(Int64, 1, column_size), I_lattice[1:end-1, :]) +
      hcat(I_lattice[:, 2:end], zeros(Int64, row_size, 1)) +
      hcat(zeros(Int64, column_size, 1), I_lattice[:, 1:end-1])
    probability_lattice = 1 .- (1-β).^count_lattice
    hit_lattice = probability_lattice .> rand(row_size, column_size)
    stage_lattice[hit_lattice] .= 1
    if sum(stage_lattice) ≥ row_size*column_size
      colormap_SI = cgrad([colorant"#111111", colorant"#EEEEEE"])
    end

    push!(time_evolution, sum(stage_lattice))
  end; gif(anim1, "logistic_growth_lattice.gif", fps = 12)
end

anim2 = @animate for t = 1:max_iteration
  compare = plot(sol,vars=(0,1),
    linestyle = :dash,  color = :black,
    label = "Theoretical", legend=:bottomright)
  plot!(compare, 0.1:0.1:(t/10), time_evolution[1:t],
    color = colorant"#111111", linewidth = 2, label = "Simulation",
    size = [400, 300])
end; gif(anim2, "logistic_growth_time_evolution.gif", fps = 12)