📂동역학

맬서스 성장 모델: 이상적인 집단 성장

모델

$$ \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 # 사망률

미분방정식으로 표현했을 땐 초기값이 고정점이 되어 변동이 없어야하지만 시뮬레이션에선 그 순간순간의 운 때문에 절멸 직전까지도 갔다가 다시 회복하기도 한다. 이 움직임은 일종의 브라우니안 모션으로 보이기도 하는데, 실제로 인구수를 늘리거나 번식률과 사망률을 내리면 조금 더 안정적인 평형상태를 유지할 것이다.

• 실제로는 기하 브라운 운동이다.

코드

다음은 이 포스트에 쓰인 줄리아 코드다.

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