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

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

Malthusian growth model

모델

01\_malthusian\_growth\_integration.gif

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

malthusian\_growth\_integration1.gif

시뮬레이션인만큼 처음에는 살짝 주춤할수도 있지만 번식률이 사망률보다 크기 때문에 결국은 폭발적인 성장이 일어난다.

$r<0$ 인 경우

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

malthusian\_growth\_integration2.gif

이론적인 절멸 수순을 거의 정확하게 밟아가는 것을 볼 수 있다.

$r=0$ 인 경우

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

malthusian\_growth\_integration3.gif

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

코드

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

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
댓글