📂동역학

로지스틱 성장 모델: 집단 성장의 한계

모델

$$ \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}}$ 의 확률로 성분을 가지게 된다.

격자 공간에서의 확산이라는 간단한 액션만을 사용했는데, 격자 공간 자체가 한정되어 있으므로 자연스럽게 성장에 제한이 걸리는 것을 확인할 수 있다.

다음은 이론적인 성장 추이와 비교한 것이다.

보다시피 시뮬레이션과 이론이 정확히 일치하지는 않는다. 우선 곡선을 그린다기보단 선형으로 성장하며, 초기 성장이 너무 가파르다. 그러나 여기서 눈여겨보아야할 것은 어쨌든 환경 용량을 다 채워가면서 꺾이는 부분이다. 이를 현실에 빗대자면 어떤 배양접시에서 세균이 번식하는 것으로 볼 수 있다. 아무리 계속 성장하고 싶더라도 배양접시에서 넘쳐흐르도록 무한히 성장할 수는 없는 것이다.

이론과 시뮬레이션이 정확히 일치하지 않는 이유는 이 시뮬레이션이 충분히 정교하지 않고 모델 역시 너무 단순하기 때문이다. 이러나 저러나 시뮬레이션에서는 그냥 운이 없었던 것일 수도 있고 격자라는 것이 주는 영향도 절대 무시할 수 없다. 단 이 경우는 단지 운 때문이 아니다.

실제로 이론적인 추이를 비교해보면 이 포스트에서 소개된 시각화는 오히려 곰페르츠 성장 모델에 가까워보인다. 실제로 격자 공간에서 종양이 성장했다고 생각해보면 체적에 대한 표면적의 비가 작아진 것으로 보이고, 세포 자체는 늘어났지만 실제로 종양의 크기를 키울 수 있는 세포가 점점 줄어가는 현상이 나타났다.

코드

다음은 움짤을 만들기 위한 줄리아 코드다.

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)