SEIR 모형 문서 원본 보기

{{위키데이터 속성 추적}}
'''SEIR 모형'''은 [[역학 (의학)|역학]]에서 [[전염병]]의 확산을 설명할 때 [[수학]]적인 모델링을 통해 접근할 수 있는 방법 중 하나이다. 이 모형은 [[SIR 모형]]보다는 더 현실에 가깝다고 여겨지며, 이는 병에 막 감염된 [[개체 (생물)|개체]]가 당장은 다른 개체에게 병을 옮기지 않는다는 점을 고려하기 때문이다. 개체에 기반한 모형보다는 거시적으로 [[개체군]] 전체를 조명하는 방식에 가깝다.

== 방정식 ==
N 개체로 이루어진 개체군을 S, E, I, R의 네 그룹으로 나누어

: <math>S+E+I+R = 1.</math>

혹은

: <math>NS+NE+NI+NR = N.</math>

의 [[공식]]이 성립되도록 한다.

모든 개체는

: Susceptible (S, 감염대상군) &#x2192; Exposed (E, 접촉군) &#x2192; Infectious (I, 감염군) &#x2192; Recovered (R, 회복군)

의 단계를 거친다. 전염병의 확산은 다음과 같이 설명될 수 있다.

: <math>\begin{align}
\frac {\mathrm{d} S}{\mathrm{d}t} &= -\beta SI,\\
\frac {\mathrm{d} E}{\mathrm{d}t} &= \beta SI-aE,\\
\frac {\mathrm{d} I}{\mathrm{d}t} &= aE-\gamma I,\\
\frac {\mathrm{d} R}{\mathrm{d}t} &= \gamma I.
\end{align}</math>

이것은 [[상미분방정식]]의 비선형 시스템에 속한다. 문헌에서는 여러 형태의 공식을 발견할 수 있는데 예를 들면  <math>\frac {\beta S I}{N}</math>항을 <math>\beta S I</math> 대신 쓰는 식이다. 이외에도 <math>S+E+I+R</math>의 합을 반드시 <math>1 </math>이라고 정하는 대신 N을 전체 개체군이라 할 때 <math>N=S+E+I+R</math>라고 정의할 수도 있다.
{| class="wikitable"
!수량
!단위
!설명
|-
|<var>S</var>(<var>t</var>)
|1
|감염대상군 (Susceptible): 감염될 수 있으나 아직 감염되지 않은 개체의 수
|-
|<var>E</var>(<var>t</var>)
|1
|접촉군 (Exposed)'':'' 이미 감염되었으나 아직 병을 전염시키지는 않는 단계에 있는 개체의 수
|-
|<var>I</var>(<var>t</var>)
|1
|감염군 (Infectious): 병을 전염시킬 수 있는 개체의 수
|-
|<var>R</var>(<var>t</var>)
|1
|회복군 (Recovered 또는 Resistent): 병에서 회복한 개체의 수와 격리 중 사망한 개체의 수의 합
|-
|<var>t</var>
|d
|시간 (일)
|-
|<var>&#x3B2;</var>
|1/d
|감염률 (Transmission rate). 역수는 접촉 평균 시간과 같다.
|-
|<var>&#x3B3;</var>
|1/d
|회복률 (Recovery rate). 역수는 평균 전염 가능 시간과 같다.
|-
|<var>a</var>
|1/d
|역수는 평균 지연 시간 (latency)과 같다.
|}
평균 지연 시간은 접촉군 (Exposed)의 한 개체가 병을 전염시킬 수 있게 되어 감염군으로 분류되기까지의 평균 시간을 나타내며 [[잠복기]]와는 다른데 이는 증상의 발현과 병을 전염시킬 수 있는 능력이 갖춰지는 것을 동치할 수 없기 때문이다. 감염률은 접촉률 (contact rate)로 불리기도 하는데, 더 상세하게는 <math>\beta=pc</math> 로 <math>p</math>는 감염 확률을, c는 접촉률을 나타낸다.

== 기초감염재생산수와의 관계 ==
병의 확산이 멈추려면 <math>\frac {\mathrm{d} E}{\mathrm{d}t}=0</math> 와 <math>\frac {\mathrm{d} I}{\mathrm{d}t}=0</math> 의 조건이 성립해야 한다. 이 조건들을 방정식에 대입하면 <math>aE = \beta SI</math> 와 <math>aE=\gamma I</math>가 되므로 <math>\beta S = \gamma</math>라는 결과가 도출된다. 이는 [[기초감염재생산수]] <math>R_0=1</math>의 경우와 같고 <math>S_0=N</math>이라고 하면<ref name="driesche">[https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6002118/ Pauline van den Driessche, Reproduction numbers of infectious disease models], Infectious Disease Modelling, Band 2, August 2017, S. 288–303</ref><ref>Odo Diekmann, Hans Heesterbeek, Tom Britton, Mathematical tools for understanding infectious disease dynamics, Princeton UP, 2013, S. 35</ref>

: <math>R_0 = \frac {\beta N} {\gamma}</math>

가 된다. <math>N=1</math>일 때 [[정규화]]된 인명 수나 변형된 SEIR 공식을 사용하면 (SIR 모형에도 적용할 수 있으나 <math>R_0</math>를 구하는 공식에서 <math>\beta</math> 의 정의가 다르다)

: <math>R_0 = \beta/\gamma</math>.<ref name="idm-seir">Institute for Desease Modeling: [https://www.idmod.org/docs/hiv/model-seir.html SEIR Model] {{웹아카이브|url=https://web.archive.org/web/20200507213249/https://www.idmod.org/docs/hiv/model-seir.html}}</ref><ref>C. Hubbs (2020): [https://towardsdatascience.com/768292f04296 Social Distancing to Slow the Coronavirus]{{깨진 링크|url=https://towardsdatascience.com/768292f04296 }}</ref>

중간 단계로의 이행을 나타내는 매개변수 <math>a</math>는 접촉군에 속한 환자가 사망할 가능성을 무시할 경우 기초감염재생산수에 영향을 주지 않는다.<ref name="driesche" />

시간에 영향을 받는 순감염재생산수 (net reproduction number''')'''는 <math>q = 1 - S</math> 일 때, <math>R_q = R_0 \cdot S</math> 또는 <math>R_q = \left( 1 - q \right) \cdot R_0</math> 로 정의할 수 있다.

== 최대 감염자 비율 ==
기초감염재생산수는 개체군의 최대 감염자 비율에 큰 영향을 준다. 병에 감염된 개체 전체를 <math>J</math>라고 하고 <math>J(t) = E(t) + I(t)</math>이면 위의 [[미분방정식]]에 의해

: <math>
\frac{\mathrm{d}J}{\mathrm{d}S}
= \frac{\mathrm{d}E + \mathrm{d}I}{\mathrm{d}S}
= \frac{-\gamma I + \beta SI}{-\beta SI}= \frac{\gamma}{\beta S} - 1
</math>,

또는 기초감염재생산수 <math>R_0 = \frac{\beta}{\gamma}</math> 일 때

: <math>
\frac{\mathrm{d}J}{\mathrm{d}S} = \frac{1}{R_0 S} - 1.</math>

가 된다. 이 방정식을 <math>\mathrm{d}J = \frac{\mathrm{d}S}{R_0 S} - \mathrm{d}S</math> 로 바꾸고 [[변수분리법]]을 통한 [[적분]]에 의해 모든 ''t''에 적용되는

: <math>J(t) + S(t) - \frac{1}{R_0} \ln S(t) = \text{const}</math>

가 도출되는데 여기서 <math>\ln</math>은 [[자연로그]]를 뜻한다. 따라서

: <math>J(t) + S(t) - \frac{1}{R_0} \ln S(t) = J(0) + S(0) - \frac{1}{R_0} \ln S(0)</math>

혹은

: <math>J(t) = J(0) + S(0) - \frac{1}{R_0} \ln S(0) + \underbrace{\frac{1}{R_0} \ln S(t) - S(t)}_{=: \ f(S(t))}.</math>

가 성립한다. 도움 함수 <math>f(x) = \frac{1}{R_0} \ln x - x</math>는 <math>f'(x) = \frac{1}{R_0 x} - 1</math> 이므로 최대값이 <math>x = \frac{1}{R_0}</math>이다. 최대 감염자 비율 <math>J_{\max}</math>에는 <math>S(t) = \frac{1}{R_0}</math>일 때 도달하게 되며, 따라서 <math>J_{\max}</math>는 기초감염재생산수와 초기값 <math>J</math> 및 <math>S</math> 만의 영향을 받는다:

: <math>
J_{\text{max}}
= J(0) + S(0) - \frac{1}{R_0} \ln S(0) -\frac{1}{R_0} + \frac{1}{R_0} \ln \frac{1}{R_0}
= J(0) + S(0) + \frac{1}{R_0} \left( \ln \frac{1}{R_0 S(0)} - 1 \right).
</math>

아직 알려지지 않은 [[바이러스]] 등에 의해 새로이 등장하는 전염병의 경우 <math>S(0) = 1</math> 이고 <math>J(0) = 0</math>으로 설정되고 최대 감염자 비율과 기초감염재생산수의 관계는 다음과 같아진다.

: <math>
J_{\text{max}}
= 1 + \frac{1}{R_0} \left( \ln \frac{1}{R_0} - 1 \right).
</math>

이 공식은 SIR 모형의 최대 감염자 비율을 구하는 공식과 일치한다.

== 전염병 유행 말기 회복된 개체의 비율 ==

격리가 전혀 없이 전염병이 유행할 경우 개체군 내에서 감염된 전체 개체의 비율도 기초감염재생산수와 상관 관계가 있다. 미분방정식으로 나타낼 경우

: <math>\frac{\mathrm dS}{\mathrm dR} = \frac{S'(t)}{R'(t)} = -\frac{\beta S}{\gamma} = -R_0 S.</math>

이고 초기값이 <math>R(0)=0</math> 라면

: <math>S(t) = S(0){\mathrm e}^{-R_0 R(t)}.</math>

이다. <math>t\to\infty</math> 일 경우 <math>E=I=0</math> 이고 따라서 <math>S+R=1</math> 이 된다. 여기서

: <math>1-q = S(0){\mathrm e}^{-R_0 q}, \quad q:=\lim_{t\to\infty} R(t).</math>

를 도출할 수 있고 다음과 같은 공식에 이르게 된다.

: <math>-R_q {\mathrm e}^{-R_q} = (q-1)R_0{\mathrm e}^{(q-1)R_0} = -S(0)R_0{\mathrm e}^{-R_0}.</math>

<math>x:=-S(0) R_0 \mathrm {e}^{-R_0}</math> 이고 <math>y:=(q-1) R_0</math> 라고 했을 때 최종적인 공식은 <math>y \mathrm{e}^{y} = x</math>가 되고 [[람베르트 W 함수]] <math>W(x)</math>를 활용해 <math>y = W(x)</math> 로 바꾼 뒤 역치환하는 과정을 거쳐

: <math>q = 1+\frac{1}{R_0}W(-S(0)R_0{\mathrm e}^{-R_0})</math>

에 도달하게 된다.

W 함수의 필요한 부분을 계산하려면 <math>f(w):=w\mathrm e^w-x</math> 함수의 2차 [[테일러 정리|테일러 다항식]]에 &#x2212;1 을 넣고 영점을 구한 뒤 [[고정 소수점 반복]]을 <math>w\mapsto x\mathrm e^{-w}</math> 로 한차례 적용한다. 그 결과물은

: <math>w = x\exp(1-\sqrt{2\mathrm e x+2}).</math>

이다. W 함수의 가까운 [[근사치]]는

: <math>W(x)\approx \varphi^n(w),</math>

에 의해 얻을 수 있고 여기서 <math>\varphi^n</math>는 [[뉴턴 방법]]의 <var>n</var> 번째 반복인

: <math>\varphi(w) = w-\frac{f(w)}{f'(w)} = \frac{x\mathrm e^{-w}+w^2}{w+1}</math>

을 뜻한다. 현실의 전염병에 적용할 때는 <math>n=4</math> 이면 충분하다.

== 기하급수적 초기 단계 ==

전염병의 초기 단계에서는 확산이 기하급수적인 형태에 가깝게 일어난다. <math>I=I(0){\mathrm e}^{\lambda t}</math> 을 가정한다면 <math>I'=\lambda I</math> 역시 참이 된다. 그에 따라 <math>aE = \gamma I+\lambda I</math> 와 <math>R=\gamma I/\lambda</math> 이 성립하며 따라서

: <math>1-S = E+I+R = \frac{\gamma+\lambda}{a}I+I+\frac{\gamma}{\lambda}I = \Big(\frac{\gamma+\lambda}{a}+1+\frac{\gamma}{\lambda}\Big)I.</math>

라는 결과가 나온다. 이 방정식을 양쪽 다 [[미분]]한 후 <math>-S' = \beta SI</math> 로 치환하고 <math>I</math>로 나누면

: <math>\beta S = \Big(\frac{\gamma+\lambda}{a}+1+\frac{\gamma}{\lambda}\Big)\lambda</math>

혹은

: <math>a\beta S = (\lambda+a)(\lambda+\gamma).</math>

를 얻게 된다. 초기에는 <math>S\approx 1</math> 가 상당히 정확한 근삿값이므로 <math>S=1</math> 이 성립하고 매개변수 <math>a,\beta,\gamma</math> 와 [[증식상수]] <math>\lambda</math> 의 관계를 알아낼 수 있다.

다른 체계적 접근법으로는 <math>S</math>를 곧바로 '''1'''로 정의하고 미분방정식을 [[연립 일차 방정식]]으로 단순화하여

: <math>\begin{pmatrix}E'\\ I'\end{pmatrix} = \begin{pmatrix}-a & \beta\\ a & -\gamma\end{pmatrix}\begin{pmatrix}E\\ I\end{pmatrix}</math>

로 설명하는 것이 가능하다.<ref>Junling Ma: [https://www.sciencedirect.com/science/article/pii/S2468042719300491 ''Estimating epidemic exponential growth rate and basic reproduction number'']. In: Infectious Disease Modelling, Volume 5, 2020, S. 129–141, KeAi Publishing (17. Dec. 2019).</ref> 다른 접근 방법들도 고정계수와 [[고윳값과 고유 벡터|고유값&고유벡터]]를 활용한 선형 상미분방정식에 기반한다. 증식상수는 [[계수행렬]] <math>J:=\big(\begin{smallmatrix}-a & \beta\\ a & -\gamma\end{smallmatrix}\big)</math>의 고유값으로 이미 발견한 방정식은

: <math>0=\det\left(J-\lambda \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}\right)=(a+\lambda)(\gamma+\lambda)-a\beta.</math>

에 의해 구할 수 있다. 모든 기하급수적 증식과정에서와 마찬가지로 증식상수 <math>\lambda</math>는 더 이해가 쉽게

: <math>T_n = \frac{\ln(n)}{\lambda},</math>

를 통해 수 배의 증식이 이루어지기까지의 시간으로 나타낼 수 있으며, 이는 특히 [[배가시간]] <math>T_2</math>에도 적용 가능하다.

== 예시 ==
다음은 [[독일의 코로나19 범유행]]을 모델링할 때 쓰인 매개변수를 [[파이썬]]에서 그대로 적용한 예시이다.<ref name="dgepi-2020-03-19">{{웹 인용|url=http://docs.dpaq.de/16007-stellungnahme2020corona_dgepi_20200318__2_.pdf|title=Stellungnahme der Deutschen Gesellschaft für Epidemiologie (DGEpi) zur Verbreitung des neuen Coronavirus (SARS-CoV-2)|format=PDF|hrsg=[[Deutsche Gesellschaft für Epidemiologie]]|datum=2020-03-18|abruf=2020-06-18|확인날짜=2020-06-18|보존url=https://web.archive.org/web/20200320070905/http://docs.dpaq.de/16007-stellungnahme2020corona_dgepi_20200318__2_.pdf|보존날짜=2020-03-20|url-status=dead}} Die aktuelle Version siehe [https://www.dgepi.de/de/aktuelles/article/aktualisierte-stellungnahme-der-deutschen-gesellschaft-fuer-epidemiologie-dgepi-zur-verbreitung-des-neuen-coronavirus-sars-cov-2/109 hier] {{웹아카이브|url=https://web.archive.org/web/20200529023927/https://www.dgepi.de/de/aktuelles/article/aktualisierte-stellungnahme-der-deutschen-gesellschaft-fuer-epidemiologie-dgepi-zur-verbreitung-des-neuen-coronavirus-sars-cov-2/109}}.</ref> 전염병의 확산은 기초감염재생산수 <math>= 2.4</math>를 토대로 이루어지고, 이는 중요한 격리 조치들이 취해지지 않았을 경우에 해당한다. 기초값 문제를 수치적으로 푸는데는 [[오일러 방법]]이면 충분하다.

<syntaxhighlight lang="python">
from numpy import array as vector

# Explicit Euler method
def euler_method(f,t0,x0,t1,h):
    t = t0; x = x0
    a = [[t,x]]
    for k in range(0,1+int((t1-t0)/h)):
        t = t0 + k*h
        x = x + h*f(t,x)
        a.append([t,x])
    return a

def SEIR_model(beta,gamma,a):
    def f(t,x):
        S,E,I,R = x
        return vector([
            -beta*S*I,
            beta*S*I - a*E,
            a*E - gamma*I,
            gamma*I
        ])
    return f

def SEIR_simulation(beta,gamma,a,E0,I0,days,step=0.1):
    x0 = vector([1.0-E0-I0,E0,I0,0.0])
    return euler_method(SEIR_model(beta,gamma,a),0,x0,days,step)

def diagram(simulation):
    import matplotlib.pyplot as plot
    plot.style.use('fivethirtyeight')
    figure,axes = plot.subplots()
    figure.subplots_adjust(bottom = 0.15)
    axes.grid(linestyle = ':', linewidth = 2.0, color = "#808080")
    t,x = zip(*simulation())
    S,E,I,R = zip(*x)
    axes.plot(t,S, color = "#0000cc")
    axes.plot(t,E, color = "#ffb000", linestyle = '--')
    axes.plot(t,I, color = "#a00060")
    axes.plot(t,R, color = "#008000", linestyle = '--')
    plot.show()

def simulation1():
    N = 83200000 # Einwohnerzahl von Deutschland 2019/2020
    R0 = 2.4; gamma = 1/3.0
    return SEIR_simulation(
        beta = R0*gamma, gamma = gamma, a = 1/5.5,
        E0 = 40000.0/N, I0 = 10000.0/N, days = 140)

diagram(simulation1)
</syntaxhighlight>

{| class="wikitable"
|[[파일:SEIR-Simulation.svg|460x460픽셀|Dynamics of an infectious disease, according to the basic SEIR model]]
S, E, I, R의 네 군이 전체 개체군에서 차지하는 비율을 날짜별로 계산한 결과.
<span style="color: #0000cc"><var>S</var>는 청색</span>, <span style="color: #f0a000"><var>E</var>는 노란색 점선</span>, <span style="color: #a00060"><var>I</var>는 자홍색</span>, <span style="color: #006000"><var>R</var>는 녹색 점선으로 표시</span>.
|}
여기서 볼 수 있는 것은 전염병이 면역임계치 (critical immune threshold)

: <math>q_c = 1-\frac{1}{R_0} = 58\,\%</math>

에 도달한 다음에도 감염군에 의해 계속 전파된다는 점이다. 다만 전염병이 면역임계치에 도달한 후에는 약 한달간의 강력한 격리 조치를 통해 멈추도록 만들 수 있을 것이다.

== 인구통계역학을 포함할 경우 ==
[[파일:SEIR_state_diagram.svg|섬네일|상태도]]
사망률 <math>\mu</math>를 상수로 두고 이에 상응하는 출생률과 함께 모형을 확장시킨

: <math>\begin{align}
S' &= \mu-\beta SI-\mu S,\\
E' &= \beta SI-aE-\mu E,\\
I' &= aE-\gamma I-\mu I,\\
R' &= \gamma I-\mu R
\end{align}</math>

공식도 있다. 원래 SEIR 모형과는 다르게 전염병의 풍토병화를 전제하며

: <math>S'=E'=I'=R'=0</math>

라고 정의한 평형에 이르기까지 감염대상군 <math>S</math>에서는 [[진동]]이 올 수도 있다.

기초감염재생산수는 여기서 평형 상태를 놓고 볼 때 다음과 같이 정의할 수 있다.

: <math>R_0 = \frac{\beta a}{(\mu+a)(\mu+\gamma)}.</math>

== 시간의존적 전염률 (Time-dependent transmission rate) ==

전염병을 대하는 태도의 변화, 격리와 계절성은 전염률의 변화를 가져오고, 이는 전염률을 시간의존적 함수 <math>\beta(t)</math>로 모델링할 때 고려되는 사안이다. 이 때에도 다른 모형들은 달라지지 않는다.<ref>Gerardo Chowell, Cécile Viboud, Lone Simonsen, Seyed M. Moghadas: [https://royalsocietypublishing.org/doi/10.1098/rsif.2016.0659 ''Characterizing the reproduction number of epidemics with early subexponential growth dynamics''] {{웹아카이브|url=https://web.archive.org/web/20200619055502/https://royalsocietypublishing.org/doi/10.1098/rsif.2016.0659}}. In: J. R. Soc. Interface 13: 20160659 (17. August 2016). DOI: 10.1098/rsif.2016.0659.</ref> 계절성을 가장 단순하게 모델링하는 방법으로는 사인 진동 (sine oscillation) 형태로 날씨가 추운 달에는 전염률이 높고 따뜻한 달에는 전염률이 낮게 나오는 것을 가정할 수가 있다.<ref>M. Keeling, P. Rohani: ''Modeling Infectious Diseases in Humans and Animals''. Abschnitt 5.2. (S.&nbsp;159): ''Modeling forcing in childhood infectious diseases: Measles''.</ref>

명백한 시간의존도를 가진 미분방정식 체계는 자생적이지 않고 따라서 더 이상 직접적으로 역동적인 체계라고 할 수 없다. 다만 <math>t'=1</math> 의 추가를 통해 인위적으로 자생적인 체계를 갖출 수는 있다.

== 외부 링크 ==
* [http://gabgoh.github.io/COVID/index.html SEIR 모형 계산기] 전염병의 경과를 SEIR 모형과 변수의 자유로운 입력을 통해 시뮬레이션하기

== 참고 문헌 ==

* Matt J. Keeling, Pejman Rohani: ''Modeling Infectious Diseases in Humans and Animals''. Princeton University Press, 2008.
* Matthias an der Heiden, Udo Buchholz: [https://www.rki.de/DE/Content/InfAZ/N/Neuartiges_Coronavirus/Modellierung_Deutschland.html ''Modellierung von Beispielszenarien der SARS-CoV-2-Epidemie 2020 in Deutschland'']. Robert Koch-Institut, Abteilung für Infektionsepidemiologie (20. März 2020). [[doi:10.25646/6571.2|DOI:10.25646/6571.2]].

== 각주 ==
<references />

[[분류:역학 (의학)]]
[[분류:수리생물학]]