길레스피 알고리즘

1. 개요

길레스피 알고리즘 Gillespie algorithm (SSA)
푸는 것	화학 마스터 방정식의 정확한 표본 경로
상태	각 화학종의 분자 개수 (정수)
한 걸음	대기시간 $\tau\sim\mathrm{Exp}(a_0)$, 반응 $j$ 를 확률 $a_j/a_0$ 로 선택
비용	반응 한 번당 한 단계 — 시간당 $a_0$ 단계
원전	길레스피 (1976, 1977); 확률론 쪽 선례: 도브 (1945), 켄들 (1950)
개량	다음-반응법, 합성-기각, 타우 도약, 하이브리드

반응속도식은 분자 수가 $10^{23}$ 개일 때의 평균 이야기다. 세포 안에 mRNA 가 다섯 개뿐이면 「농도」라는 말부터 성립하지 않는다.

길레스피 알고리즘(Gillespie algorithm), 정식 이름 확률적 시뮬레이션 알고리즘(stochastic simulation algorithm, SSA)은 잘 섞인 반응 용기 안의 화학 반응계를 분자 하나하나의 반응 사건 단위로 시뮬레이션하는 몬테카를로 방법이다. 상태는 각 화학종의 분자 수 $\mathbf{x} = (x_1, \dots, x_S) \in \mathbb{Z}^S_{\ge 0}$ 이고, 반응이 한 번 일어날 때마다 정수만큼 바뀐다. 다니엘 길레스피가 1976년 계산물리 논문과 1977년 물리화학 논문으로 내놓았다. 같은 계산 구조는 확률 과정 이론에서 도브(1945)·켄들(1950)이 이미 썼지만¹, 화학 반응계에 맞춘 정식화와 「이것이 화학 마스터 방정식의 정확한 표본을 준다」는 논증이 길레스피의 몫이다.

결정론적 반응속도론의 상미분방정식은 분자 수가 많을 때 이 확률 과정의 평균이 따르는 근사다. 분자 수가 적거나(유전자 발현, 이온 채널, 촉매 표면), 요동 자체가 관심사이거나(쌍안정 스위치, 멸종), 비선형이라 평균과 평균의 함수가 다를 때 SSA 가 필요하다. 이벤트 구동 시뮬레이션 문서가 SSA 를 이벤트 구동의 한 사례로 한 단락 소개했고, 포아송 분포 문서가 근사판인 타우 도약을 다뤘다. 이 문서는 알고리즘 자체와 그 정당성, 검증 방법, 그리고 비용 문제를 다룬다.

2. 화학 마스터 방정식과 성향 함수

반응 채널 $j = 1, \dots, M$ 마다 상태 변화 벡터 $\nu_j$ 와 성향 함수(propensity) $a_j(\mathbf{x})$ 가 있다. $a_j(\mathbf{x})\,dt$ 는 「현재 상태가 $\mathbf{x}$ 일 때 다음 $dt$ 동안 반응 $j$ 가 한 번 일어날 확률」이다. 잘 섞인 기체·희박 용액 가정에서 이것이 물리적으로 정당하다는 것이 길레스피(1992)의 유도다.

성향 함수는 질량작용 법칙의 조합적 판이다. 확률적 속도 상수를 $c_j$ 라 할 때

반응	성향 $a_j(\mathbf{x})$
$\varnothing \to \cdots$	$c_j$
$X \to \cdots$	$c_j\,x$
$X + Y \to \cdots$	$c_j\,xy$
$2X \to \cdots$	$c_j\,x(x-1)/2$

이다. 마지막 줄의 $x(x-1)/2$ 는 「서로 다른 두 분자 쌍의 개수」라 $x^2$ 가 아니다. 결정론적 속도 상수 $k$ 와의 관계는 부피 $V$ 와 아보가드로 수 $N_A$ 를 통해 단분자 $c = k$, 이분자 이종 $c = k/(N_AV)$, 이분자 동종 $c = 2k/(N_AV)$ 다(동종의 2 는 속도식을 $-\tfrac{d[X]}{dt}$ 로 쓰느냐 반응 속도로 쓰느냐의 관례에 따라 달라지므로 출처의 관례를 반드시 확인할 것)².

상태 $\mathbf{x}$ 에 있을 확률 $P(\mathbf{x}, t)$ 는 화학 마스터 방정식(CME)

$$\frac{\partial P(\mathbf{x}, t)}{\partial t} = \sum_{j=1}^M \Big[a_j(\mathbf{x} - \nu_j)P(\mathbf{x} - \nu_j, t) - a_j(\mathbf{x})P(\mathbf{x}, t)\Big]$$

를 따른다. 상태 공간이 종 수에 대해 지수적으로 커서(종 10개, 각 0~100개면 $10^{20}$ 상태) 직접 풀기는 대개 불가능하다. SSA 는 이 방정식을 풀지 않고 이 방정식이 기술하는 연속시간 마르코프 연쇄의 표본 경로를 정확히 생성한다. 통계량은 경로를 여러 번 뽑아 평균한다.

3. 직접법

현재 상태에서 총 성향을 $a_0 = \sum_j a_j(\mathbf{x})$ 라 하자. 다음 반응까지의 대기시간 $\tau$ 와 그 반응의 번호 $j$ 의 결합 밀도는

$$p(\tau, j) = a_j\,e^{-a_0\tau}$$

이고, 이것은 $\tau \sim \mathrm{Exp}(a_0)$ 와 $P(j) = a_j/a_0$ 의 독립 곱으로 인수분해된다. 그래서 한 걸음은:

1. 모든 $a_j(\mathbf{x})$ 와 $a_0$ 를 계산한다.
2. 균등 난수 $r_1, r_2$ 를 뽑는다.
3. $\tau = -\ln(r_1)/a_0$.
4. $\sum_{i<j} a_i < r_2a_0 \le \sum_{i\le j} a_i$ 인 $j$ 를 찾는다.
5. $t \leftarrow t + \tau$, $\mathbf{x} \leftarrow \mathbf{x} + \nu_j$. 1로.

근사가 하나도 없다. 시간 간격 $dt$ 를 고르는 것도 없다. 그래서 「정확한」 확률적 시뮬레이션이라 부른다 — 오차는 오직 표본 수에서 오는 몬테카를로 오차다.

길레스피는 같은 논문에서 첫-반응법(first reaction method)도 제시했다. 채널마다 $\tau_j \sim \mathrm{Exp}(a_j)$ 를 따로 뽑아 가장 작은 것을 고른다. 지수분포 최솟값의 성질로 직접법과 같은 분포를 주지만 난수를 $M$ 개 쓰므로 그 자체로는 더 비싸다. 이것이 다음 절의 개량의 출발점이다.

4. 비용과 개량

직접법 한 걸음의 비용은 성향 갱신 $O(M)$ 과 선택 $O(M)$ 이다. 그리고 걸음 수는 반응 사건 수와 같아서, 시뮬레이션 시간당 평균 $a_0$ 걸음이 든다. 비용 문제는 두 갈래다.

채널이 많을 때($M$ 큼).

• 다음-반응법(깁슨·브룩 2000). 첫-반응법의 예약 시각들을 우선순위 큐(인덱스 힙)에 넣고, 반응 의존 그래프로 실제로 성향이 바뀐 채널만 갱신한다. 바뀐 채널의 예약 시각은 새 난수 없이 $\tau_j^{\text{new}} = (a_j^{\text{old}}/a_j^{\text{new}})(\tau_j^{\text{old}} - t) + t$ 로 재사용할 수 있다(지수분포의 무기억성). 걸음당 $O(\log M)$.
• 합성-기각(슬레포이·톰프슨·플림프턴 2008). 성향을 크기에 따라 2의 거듭제곱 구간 묶음으로 나누고, 묶음을 고른 뒤 묶음 안에서 기각으로 고른다. 묶음 수가 성향의 동적 범위의 로그로 묶이므로 걸음당 기대 비용이 사실상 $O(1)$ 이다.
• 선택 단계만 바꾸려면 누적합을 이진 트리(펜윅 트리 등)로 들고 다녀 $O(\log M)$ 으로 할 수도 있다.

사건이 너무 잦을 때($a_0$ 큼). 분자 수가 많거나 빠른 가역 반응이 있으면 대부분의 걸음이 의미 없는 왕복에 쓰인다. 이것이 확률적 시뮬레이션에서의 강성이다.

• 타우 도약(길레스피 2001). 간격 $\tau$ 동안 성향이 거의 일정하다고 보고 채널별 발생 횟수를 포아송으로 한 번에 뽑는다. 음수 분자 문제와 $\tau$ 선택은 포아송 분포 문서와 이항 분포의 이항 도약 절을 볼 것. $\tau$ 선택의 표준은 카오·길레스피·페촐드(2006)의 「성향의 상대 변화가 $\epsilon$ 이하」 규칙이다.
• 화학 랑주뱅 방정식. 포아송을 정규로 근사하면 확률미분방정식이 되고, 잡음까지 없애면 결정론적 반응속도식이 된다. SSA → 타우 도약 → CLE → ODE 가 분자 수가 늘어날 때의 근사 사다리다.
• 하이브리드·느린-빠른 분리. 빠른 반응은 준평형으로 두고 느린 반응만 SSA 로 돌리거나, 분자가 많은 종은 ODE·적은 종은 SSA 로 섞는다. 메커니즘 축소의 준정상 상태 가정과 같은 발상이다.

5. 검증 — 닫힌꼴이 있는 문제부터

SSA 코드는 「결과가 그럴듯해 보이는」 것으로는 검증이 안 된다. 성향 공식의 계수 하나만 틀려도 결과는 여전히 매끄럽게 요동친다. 해가 알려진 문제로 분포까지 맞춰야 한다.

출생-사망 과정. $\varnothing \xrightarrow{k} X$, $X \xrightarrow{\gamma} \varnothing$. CME 의 정상해는 정확히 $\mathrm{Poisson}(k/\gamma)$ 다 — 평균과 분산이 모두 $k/\gamma$, 판노 인자(분산/평균) 1. $k = 10$, $\gamma = 1$ 로 시간 $2\times10^4$ 동안 직접법을 돌려(초기 100 은 버림) 시간 가중 평균을 내면 평균 10.009, 분산 10.078, 판노 인자 1.007 이었다. 이 동안 반응 사건이 약 40만 번 일어났는데, 정상 상태 $a_0 = k + \gamma\bar x \approx 20$ 이니 단위 시간당 20 걸음으로 맞는다. 성향을 $\gamma x$ 대신 $\gamma x^2$ 처럼 틀리게 넣으면 판노 인자가 1 에서 바로 벗어난다.

확률적 멸종 — ODE 가 못 보는 것. 감염병 구획 모형 SIR 을 반응계로 쓰면 $S + I \xrightarrow{\beta/N} 2I$, $I \xrightarrow{\gamma} R$ 이다. $N = 1000$, 기초감염재생산수 $R_0 = \beta/\gamma = 2$, 초기 감염자 1명으로 SSA 를 4000번 돌리면 49.9% 의 경로에서 감염자 수가 100 명도 못 가 소멸한다. 분기 과정 이론의 예측 「대유행 확률 $= 1 - 1/R_0 = 0.5$」 와 일치한다. 대유행으로 간 경로들의 최종 감염 비율 평균은 0.795 로, 결정론적 최종 크기 방정식 $z = 1 - e^{-R_0z}$ 의 해 0.7968 과 맞는다. 결정론적 ODE 는 초기 감염자가 1 명이든 0.001 명이든 언제나 대유행을 예측한다 — 멸종 확률 50% 라는 정보 자체가 평균 방정식에는 없다.

이 두 시험은 각각 「요동의 크기」와 「평균이 못 보는 사건」을 검사한다. 쌍안정 계(슐뢰글 모형, 토글 스위치)의 상태 전이 시간 분포는 세 번째 표준 시험이다.

6. 어디에 쓰이나

• 유전자 발현 잡음. mRNA·단백질이 몇 개~수십 개인 계에서 세포 간 변이(내재적 잡음)의 표준 도구. 「전사 폭발」(bursting)이 만드는 과분산 단백질 분포를 모형과 비교할 때 SSA 표본이 기준선 역할을 한다.
• 이온 채널. 막의 채널 개수가 유한하면 개폐가 확률적이고, 그 요동이 활동전위 타이밍을 흔든다(이온 채널 참고).
• 역학·생태. 감염병, 개체군 동역학의 멸종·고정 확률.
• 운동 몬테카를로. 표면 반응·결정 성장·확산의 운동 몬테카를로(kinetic Monte Carlo, BKL 알고리즘)는 구조가 같다. 사건 목록에서 성향 비례로 하나를 고르고 지수 대기시간만큼 시계를 진행한다. 재료 쪽과 화학 쪽이 같은 알고리즘을 서로 다른 이름으로 부르는 경우다³.

7. 관련 문서

• 반응속도론 · 메커니즘 축소 · 화학 마스터 방정식 · 운동 몬테카를로
• 포아송 분포 · 이항 분포 · 확률미분방정식 · 몬테카를로 방법 · 난수 생성기
• 이벤트 구동 시뮬레이션 · 우선순위 큐 · 구획 모형 · 이온 채널 · 강성 방정식

8. Footnotes

1. 켄들의 1950년 논문 제목이 「단순 출생-사망 과정의 인공적 실현」이다. 난수로 표본 경로를 만들어 보인, 사실상 전산기 보급 이전의 SSA 다. 지금 노트북으로 1초에 수십만 사건을 뽑는 걸 보면 켄들이 킹받을 만하다. ↩

2. 동종 이분자 반응의 2 배 인자는 SSA 코드 버그 랭킹 상위권 단골이다. 평균이 ODE 와 정확히 두 배 어긋나는데, 비선형이라 「대충 비슷한데?」 하고 넘어가기 딱 좋다. 출생-사망처럼 해가 알려진 이량체화 문제로 한 번만 돌려 보면 잡힌다. ↩

3. 같은 알고리즘이 화학에서는 길레스피, 재료에서는 BKL(보츠·칼로스·레보위츠 1975), 확률론에서는 그냥 「연속시간 마르코프 연쇄 시뮬레이션」이라 불린다. 세 분야가 한동안 서로의 논문을 거의 인용하지 않고 각자 발전시켰다는 게 학계 사일로의 국룰 예시. 지금도 문헌 검색을 할 땐 세 이름을 다 넣어야 한다. ↩