포식자-피식자 위상공간은 피식자 개체수를 가로축, 포식자 개체수를 세로축에 놓고 시간 순서대로 점을 찍어 이은 그래프다. 같은 모델에서 난수 시드만 100번 갈아 끼우면, 개체수는 그중 11번만 제자리로 돌아온다.
늑대와 양, 5일의 지연
화면 안에서 양이 풀을 뜯고 늑대가 양을 쫓는다. 옆에 그려지는 개체수 곡선 두 개는 60일 내내 크게 출렁이는데, 늑대 곡선이 양 곡선을 며칠 늦게 따라온다. 여기까지는 두 종 시뮬레이션을 소개할 때 흔히 보는 화면이다. 그런데 모델도 파라미터도 그대로 두고 난수 시드만 갈아 끼워 100번을 돌려 보면, 이 화면은 11번밖에 나오지 않는다.
숫자로 보면 이렇다. 피식자는 첫날 67마리에서 시작해 48일째 202마리까지 불어났다가 56일째 11마리까지 줄어든다. 포식자는 다르게 움직인다. 첫날 18마리에서 45일째 5마리까지 내려갔다가 54일째 54마리로 올라선다. 두 곡선을 하루씩 밀어 가며 상관계수를 재면 5일 밀었을 때 0.8864로 가장 높고, 밀지 않으면 −0.1548이다. 오늘의 포식자 수를 오늘의 피식자 수로는 설명할 수 없다.
그림 1는 이 시뮬레이션의 출하 설정, 시드 42를 60일 돌린 결과다. 앞에서 말한 100번은 이 설정에서 난수 시드만 갈아 끼운 것이고, 두 종이 이렇게 함께 살아남아 출렁인 시드가 11개다. 나머지 89개는 60일 안에 무너진다. 그중 82개에서는 포식자가 먼저 사라지고, 남은 7개는 두 종이 함께 무너진다. 무너지는 쪽은 마지막 절에서 센다. 살아남은 쪽부터 뜯어보면, 곡선을 이렇게 밀어 놓는 것이 무엇인지는 개체 하나의 행동 규칙에 다 들어 있다.
국소 규칙만 아는 개체들
개체 하나가 아는 것은 자기 주변뿐이다. 피식자는 감지 범위 안의 풀을 향해 움직이다가도, 포식자가 감지 범위에 들어오면 먹이를 두고 반대 방향으로 도망친다. 도망 반경(20)이 먹이 감지 반경(18)보다 넓어서, 먹이가 훨씬 가까이 있어도 도망이 먼저다. 포식자는 감지 범위 안의 피식자를 쫓고, 따라잡으면 사냥한다.
두 종 모두 매 틱 대사로 에너지를 잃고, 에너지가 바닥나면 굶어 죽는다. 하루는 200틱이고, 하루가 끝나는 시점에 에너지가 번식 임계를 넘은 개체만 새끼를 남긴다. 시작은 피식자 70마리·포식자 14마리다.
어느 규칙에도 "주기를 만들라"는 말이 없다. 개체는 자기 감지 반경 밖을 보지 못하고, 전체 개체수라는 숫자를 아는 개체도 없다. 그런데도 전역 개체수가 규칙적으로 출렁인다.
5일이라는 지연은 규칙이 순서대로 이어지는 데 걸리는 시간이다. 포식자가 피식자를 한 마리 잡는다고 포식자가 곧바로 늘지는 않는다. 사냥으로 얻은 에너지가 임계까지 쌓여야 하고, 그다음 하루 경계에서야 새끼가 태어난다. 피식자가 많아진 시점과 포식자가 많아지는 시점 사이에 이 축적 시간이 통째로 끼어든다.
그림 2가 위상 지연이 생기는 자리다. 설계자가 넣은 주기가 아니라, 먹은 것이 새끼가 되기까지 걸리는 시간만큼 두 곡선이 어긋난다. 지연 0에서 상관계수가 음수로 나오는 것도 같은 이유다. 시간축 그림에서는 이 어긋남이 "며칠 밀렸다"로만 보인다.
포식자-피식자 위상공간으로 접기 — 닫힌 고리의 등장
지금까지 본 그림은 가로축이 시간이다. 축을 바꾸면 같은 데이터가 다른 모양이 된다. 새로 재는 숫자는 없고, 이미 가진 60일치 관측을 다른 좌표에 옮겨 놓을 뿐이다.
위상공간(phase space)은 시간축을 지우고 상태 변수끼리 축을 잡는 좌표계다. 여기서는 가로축이 피식자 개체수, 세로축이 포식자 개체수다. 축이 둘이니 이 좌표계는 평면이 되고, 그래서 위상평면이라고도 부른다. 하루치 관측 한 쌍이 점 하나가 되고, 60일이면 점 60개를 시간 순서대로 이어 선 하나가 된다. 시간은 축에서 사라지고 선을 따라가는 순서로만 남는다.
이렇게 그리면 선이 한 점을 감고 돈다. 피식자가 늘고 포식자가 늘고, 피식자가 줄고 포식자가 줄어드는 순서가 평면에서는 회전으로 나타난다. 출발점 근처로 되돌아와 원처럼 이어지는 이 모양을 닫힌 고리(closed orbit) 라고 부른다.
그림 3에서 궤적은 중앙값으로 잡은 중심(피식자 28.0·포식자 13.0) 둘레를 2.691바퀴 감는다. 이렇게 감은 바퀴 수를 권선수(winding number)라고 부른다. 60일을 2.691로 나누면 한 바퀴에 22.3일이다. 시간축에서 "며칠 주기로 출렁인다"고 말하던 것이 평면에서는 "몇 바퀴 돈다"가 된다. 다만 선이 정확히 겹치지는 않는다. 그 어긋남은 뒤에서 잰다.
아래는 그림 1의 지연을 다시 계산하는 코드다. 두 시계열을 각각 전체 구간 평균으로 중심화한 뒤, 지연만큼 밀어 겹친 조각끼리 내적하고 그 조각의 노름으로 나눈다. 이 글의 상관 수치는 전부 이 규약으로 잰 값이다. 겹친 구간을 다시 중심화하는 피어슨 상관계수로 계산하면 최적 지연은 그대로 5일이지만, 값 자체는 조금 달라진다.
import numpy as np
from horongs_sims.predprey import PredPreyModel, PredPreyParams
def daily_counts(seed, days):
"""하루 = day_ticks(200) 틱 + 일 경계 처리. 개체수는 일 경계 직후 기록."""
p = PredPreyParams()
m = PredPreyModel(p, seed)
prey, pred = [], []
for _ in range(days):
for _ in range(p.day_ticks):
m.step()
m.day_boundary()
prey.append(m.n_prey)
pred.append(m.n_pred)
return np.array(prey, float), np.array(pred, float)
x, y = daily_counts(seed=42, days=60)
xc, yc = x - x.mean(), y - y.mean()
for lag in range(8):
a, b = (xc, yc) if lag == 0 else (xc[:-lag], yc[lag:])
r = a @ b / (np.linalg.norm(a) * np.linalg.norm(b))
print(f"lag={lag} r={r:+.4f}")lag=0 r=-0.1548
lag=1 r=+0.1139
lag=2 r=+0.4185
lag=3 r=+0.6884
lag=4 r=+0.8486
lag=5 r=+0.8864
lag=6 r=+0.7919
lag=7 r=+0.5994지연 0에서 −0.1548이던 값이 한 칸씩 밀 때마다 올라가 5일에서 0.8864로 꼭대기를 찍고 다시 내려간다. 5일이라는 숫자는 이 목록에서 가장 큰 값이 나온 자리를 고른 것이다. 여기까지는 관측이다.
개체를 지우면 남는 것 — 평균장 로트카-볼테라 방정식
고리를 눈으로 보는 것과 식으로 적는 것은 다른 일이다. 개체 하나하나의 위치와 에너지를 그대로 들고서는 두 줄짜리 식이 나오지 않는다. 식으로 적으려면 개체를 지워야 한다.
평균장(mean-field) 근사는 개체 하나하나를 지우고 개체수를 연속적인 양으로 다루는 방법이다. 누가 누구를 어디서 만나는지는 묻지 않는다. 단위 시간당 만남의 횟수를 두 개체수의 곱으로 놓을 뿐이다. 여기에 피식자가 스스로 늘어나는 항과 포식자가 스스로 죽는 항을 붙이면 로트카-볼테라(Lotka-Volterra, 이하 LV) 방정식이 나온다.
식 1에 들어간 기호는 여섯 개다.
- $x$ — 피식자 개체수
- $y$ — 포식자 개체수
- $\alpha$ — 피식자 성장률
- $\beta$ — 포식률
- $\delta$ — 사냥이 포식자 증식으로 바뀌는 비율
- $\gamma$ — 포식자 자연사망률
개체·공간·에너지가 전부 사라지고 비율 네 개만 남았다. 이 네 개에 어떤 값이 들어가든, 식의 모양에서 바로 읽어낼 수 있는 것이 하나 있다.
두 식이 동시에 0이 되는 지점을 고정점(fixed point)이라고 부른다. 첫 식을 0으로 놓고 풀면 $y^{\ast} = \alpha/\beta$가 나온다. 둘째 식에서는 $x^{\ast} = \gamma/\delta$다. 개체수가 이 두 값에 정확히 놓이면 두 종 모두 늘지도 줄지도 않는다. 평균장 궤적은 이 점을 중심으로 돌게 되는데, 어떤 고리를 그리며 도는지 그 고리가 매번 같은 고리인지는 아직 말하지 않았다.
보존량과 중립 고리 — 교과서 고리는 왜 영원히 도는가
교과서에 실린 LV 위상평면 그림에서 궤적은 감쇠하지 않는다. 몇 바퀴를 돌든 같은 고리를 그대로 다시 돈다. 이유는 식 1에 보존량이 하나 숨어 있기 때문이다.
보존량은 궤적을 따라가도 값이 변하지 않는 양이다. 마찰 없는 진자의 역학적 에너지가 그런 양이다. LV에서는 다음 조합이 그 역할을 한다.
식 2를 시간으로 미분하고 식 1를 대입하면 결과가 정확히 0이다. 그래서 평균장 궤적은 $V$가 같은 값을 갖는 등고선 위를 벗어나지 못한다. 이 등고선이 고정점을 둘러싼 닫힌 곡선이다. 그래서 궤적은 같은 고리를 계속 돈다.
이런 고리를 중립 고리(neutral cycle) 라고 부른다. 중립은 안쪽으로 수렴하지도 바깥으로 발산하지도 않는다는 뜻이다. 어느 고리를 돌지는 초기 조건 하나가 정하고, 그 뒤로는 바뀌지 않는다. 다음 절에서 실측 데이터에 맞춰 볼 네 파라미터를 미리 넣고 평균장 방정식을 60일 수치 적분한 뒤 매 스텝의 $V$를 재면, 표준편차가 0.0이고 변동계수가 0.0%다.
다만 이것은 이상화된 수학적 성질이다. 개체가 정수로 세어지지 않고, 만남이 확률이 아니며, 공간도 없는 세계에서만 성립한다. 앞에서 돌린 시뮬레이션은 셋 다 아니다.
우리 데이터로 적합해보기 — 요약이 깨지는 지점
지금까지 네 파라미터는 기호였다. 이제 실측 데이터에 맞춰 값을 정한다. 새로 필요한 수학은 하나뿐인데, 두 미분방정식을 직선 적합 문제로 바꾸는 변형이다. 식 1의 두 식을 각각 왼쪽 변수로 나누면 왼쪽이 로그의 미분이 된다.
식 3의 오른쪽은 둘 다 일차식이다. 왼쪽은 60일 실측 개체수에서 중앙차분으로 구할 수 있다. 그러면 남는 것은 점들에 직선 하나를 맞추는 최소제곱 문제이고, 풀면 α = 0.3028 /일, β = 0.01842, δ = 0.00337, γ = 0.21369 /일이 나온다.
여기서 적합도를 함께 적어야 한다. 피식자 식의 결정계수 R²는 0.5523, 포식자 식은 0.4693이다. 이 직선들은 일별 변화의 절반 안팎만 설명한다. 적합값 네 개만 인용하고 이 두 수를 빼면, 이 글이 마지막에 문제 삼는 바로 그 인용이 된다.
앞에서 구한 고정점 공식에 이 네 값을 넣으면 (63.45, 16.44)다. 단, 반올림하지 않은 적합값을 넣어야 이 수가 나온다. 위에 인쇄한 반올림 파라미터끼리 나누면 소수 둘째 자리가 어긋난다. 이제 앞에서 세워 둔 세 가지(고정점·주기·보존량)를 실측 궤적에 그대로 대 볼 수 있다.
먼저 시간평균이다. LV에는 한 주기 동안의 시간평균이 정확히 고정점과 같다는 정리가 있다. 전제는 "한 주기"다. 앞에서 센 60일 창은 2.691바퀴라 정수 주기가 아니고, 하필 피식자가 11마리까지 내려간 깊은 골에서 끊긴다. 이 창으로 자르면 평균장 ODE조차 자기 고정점을 1.19% 벗어난다 — 실측 궤적이 아니라 정리가 성립하는 그 방정식의 해가 그렇다.
정수 바퀴에 가장 가까운 창으로 다시 잰다. 1일부터 52일까지가 2.033바퀴다. 이 창의 실측 평균은 피식자 61.731·포식자 15.923으로, 예측값보다 각각 2.71%·3.12% 낮다. 요약은 궤적의 중심을 3% 안쪽으로 맞힌다. 60일 창 전체로 재면 피식자가 11.3% 낮고 포식자가 5.9% 높은데, 그 차이의 상당 부분은 모형이 틀려서가 아니라 정수 주기가 아닌 곳에서 창을 끊었기 때문이다.
주기도 비슷하게 맞는다. 고정점 근처에서 선형화한 주기 $2\pi/\sqrt{\alpha\gamma}$는 24.701일이고, 적분한 궤도를 직접 재면 24.727일이다. 앞에서 권선수로 환산한 실측 주기는 22.3일이다. 차이는 10% 안쪽이다.
갈리는 곳은 진폭이다. 같은 파라미터로 첫날 상태(피식자 67·포식자 18)에서 출발해 60일을 적분하면 평균장 피식자는 55.97과 71.56 사이만 오간다. 중심 대비 ±12% 폭이다. 같은 기간 실측 피식자는 11마리에서 202마리 사이를 오가는데, 최대를 최소로 나누면 18.4배다. 포식자도 평균장은 14.80~18.19인데 실측은 5~54로 10.8배다. 요약은 중심과 박자를 맞히고 크기를 놓친다.
고리가 닫히는지도 같은 식으로 잰다. 60일 실측 개체수를 하루씩 식 2에 넣어 $V$를 뽑고, 그 표준편차와 변동계수를 본다. 상수로 남지 않는다.
import numpy as np
from horongs_sims.predprey import PredPreyModel, PredPreyParams
p = PredPreyParams()
m = PredPreyModel(p, 42)
prey, pred = [], []
for _ in range(60):
for _ in range(p.day_ticks):
m.step()
m.day_boundary()
prey.append(m.n_prey)
pred.append(m.n_pred)
x = np.array(prey, float)
y = np.array(pred, float)
alpha, beta, delta, gamma = 0.3028, 0.01842, 0.00337, 0.21369
V = delta * x - gamma * np.log(x) + beta * y - alpha * np.log(y)
print(f"V mean = {V.mean():.4f}")
print(f"V std = {V.std():.4f}")
print(f"V CV = {100 * V.std() / abs(V.mean()):.2f}%")
print(f"V range= [{V.min():.3f}, {V.max():.3f}]")V mean = -1.0689
V std = 0.0863
V CV = 8.08%
V range= [-1.216, -0.752]평균장에서 0.0%였던 변동계수가 8.08%다. 시드 42를 절멸 직전까지 밀면 95일 구간 변동계수는 11.57%로 커지고, 하루당 $V$ 변화의 표준편차는 0.0751이다. 이 긴 구간을 넷으로 나눠 평균을 내면 앞에서 뒤로 갈수록 계속 올라간다. 고리가 조금씩 바깥으로 벌어진다는 뜻이다. 그리고 96일째에 포식자가 0이 된다.
원인으로 가장 먼저 의심할 것은 개체수다. 개체가 유한하면 태어나고 죽는 사건 하나하나가 잡음이 되고, 그 잡음의 상대 크기는 개체수의 제곱근에 반비례해 작아진다. 그렇다면 계를 키울수록 표류가 줄어야 한다. 밀도와 상호작용 반경을 고정한 채 면적과 개체수를 함께 키워 세 규모를 돌렸다. 평균 총개체수는 95.2 · 247.4 · 524.1이고, 하루당 $V$ 표류 폭은 0.0603 · 0.0697 · 0.0585다.
유한 개체 잡음이 원인이라면 이 세 점의 로그-로그 기울기가 −0.5로 나와야 한다. 실제 기울기는 사실상 0이다. 총개체수를 5.5배 키우는 동안 표류가 줄었다는 신호가 없다. 고리가 닫히지 않는 주된 이유는 개체가 적어서가 아니다. 평균장 LV라는 요약이 이 모형을 제대로 기술하지 못하기 때문이다. 앞에서 본 R² 0.5523·0.4693이 이미 같은 신호였다.
시드 스윕 — 공존은 예외다
여기까지는 시드 하나다. 파라미터를 그대로 두고 난수 시드만 1번부터 100번까지 바꿔 60일씩 돌렸다. 판정 순서는 이렇다. 먼저 절멸했는지 보고, 다음으로 인위적 상한에 붙었는지 본다. 권선수가 1바퀴 미만이면 "진동 없음"으로 걸러내고, 남은 것만 공존 사이클로 셌다.
그림 4의 공존율은 11.0%다. 에이전트 모델이 무너지는 모양은 크게 셋으로 나뉜다 — 폭주, 개체수 붕괴, 인위적 상한 포화다. 이 스윕에서 상한 포화는 0건이고, 나머지 89건은 전부 붕괴 쪽이다.
실패한 89건의 절멸 시점은 중앙값 30일이다. 가장 이른 것이 10일, 가장 늦은 것이 59일이다. 10일 안에 무너진 경우는 2.2%뿐이다. 실패의 대부분은 초기 배치를 잘못 뽑아 곧바로 무너진 것이 아니라, 사이클을 몇 번 돌다가 붕괴한 것이다.
붕괴 뒤에는 폭주가 온다. 포식자가 사라진 뒤 피식자의 최종 개체수는 중앙값 247마리다. 이 모형의 소프트캡이 400마리이니, 포식이 없어지면 피식자는 상한 쪽으로 밀려 올라간다.
공존으로 판정된 11개도 60일이라는 관측창 덕을 본다. 그 시드들을 200일까지 늘려 돌리면 전부 절멸한다. 출하 시드 42는 96일째에 포식자가 0이 되고, 그때까지 4.57바퀴를 돌았다. 마지막 6일 포식자 개체수는 46 → 27 → 13 → 3 → 1 → 0이다. 60일에서 끊었기 때문에 "지속 진동"으로 보였을 뿐이다.
파라미터를 흔들면 무엇이 사이클을 떠받치는지 드러난다. 기본 설정에서 공존한 시드 셋(42·4·11)을 골라 은신처만 바꿔 60일씩 다시 돌렸다. 은신처를 끄면 셋 다 포식자가 절멸하고, 절반인 7로 줄여도 결과는 같다. 1.5배인 21로 늘리면 셋 중 하나만 공존한다. 어느 쪽으로 옮겨도 기본값보다 나빠진다.
초기 포식자 수는 이 방식으로 재면 안 된다. 그 세 시드는 애초에 기본값 14에서 공존한다는 이유로 뽑은 시드라, 14에서 3/3이 나오는 것은 설계상 보장돼 있다. 값마다 같은 시드 집합을 쓰되, 시드 1번부터 30번까지 전수로 다시 셌다. 공존율은 초기 포식자 11에서 20.0%, 14에서 10.0%, 17과 20에서 각각 16.7%다. 30개 시드로 10.0%와 16.7%를 가를 수는 없으니, 말할 수 있는 것은 기본값이 특별히 유리한 자리는 아니라는 것까지다.
이 모형에는 은신처(prey_refuge=14)가 들어 있다. 피식자가 14마리 이하로 내려가면 포식이 멈추는 규칙이다. 교과서 LV에는 없는 인위적 안정자이며, 피식자 절멸을 막아 반등을 만든다. 이 규칙을 끄면 기본 설정에서 공존하던 시드마저 전부 절멸한다 — 사이클은 이 규칙 하나에 얹혀 있다.
여기까지가 모형 하나를 끝까지 뜯어본 결과다. 같은 질문은 남이 돌린 모형에도 그대로 던질 수 있다. 논문이나 영상에서 위상평면 고리 한 장을 받았을 때, 그림만 봐서는 알 수 없는 것이 다섯 가지다.
- 시드가 몇 개인가. 하나라면 그 고리는 11%에 해당하는 그림일 수 있다.
- 실패한 런은 몇이고 어떻게 실패했나. 초기 배치 운으로 곧바로 무너진 것과 사이클을 몇 바퀴 돌다 붕괴한 것은 다른 이야기다.
- 관측창이 얼마나 긴가. 시드 42는 60일에서 공존이고 96일에서 절멸이다.
- 평균장 파라미터를 인용했다면 적합도는 얼마인가. R²가 없는 α·β·δ·γ는 측정값이 아니라 그림 설명이다.
- 인위적 안정자가 들어갔는가. 은신처 같은 규칙 하나가 사이클 전체를 떠받치고 있을 수 있다.
다섯 항목이 다 채워진 결과라면 고리 한 장을 그대로 받아도 된다. 비어 있다면 그 그림은 아직 주장이 아니라 표본 하나이고, 그 고리가 규칙이 만든 것인지 시드가 만든 것인지는 채워진 다음에야 갈린다.
더 읽기
- U. Täuber, "Population oscillations in spatial stochastic Lotka-Volterra models: A field-theoretic perturbational analysis" (2012) — https://arxiv.org/abs/1206.2303
- C. Mandl, S. Chaturvedi, M. van Gerven, "Tuning Agent-Based Predator-Prey Models Toward Lotka-Volterra Dynamics" (2026) — https://arxiv.org/abs/2606.13639 (Semantic Scholar 사본: https://www.semanticscholar.org/paper/067695bfc26e2325d27b091239f7ca7422d1bf09)
- 이 글이 잰 시뮬레이션 코드 — 호롱스 자체 자산
sims/horongs_sims/predprey.py(PredPreyModel, 기본 파라미터·시드 42)