Modelo epidémico SIR

La idea

El modelo SIR de Kermack y McKendrick divide una población de tamaño N en tres compartimentos: Susceptibles, Infectados y Recuperados. Los contagios ocurren al ritmo β S I / N, proporcional a la probabilidad de que un susceptible se encuentre con un infectado, y las recuperaciones al ritmo γ I, lo que equivale a suponer que el periodo infeccioso dura en promedio 1/γ. El sistema es S′ = −β S I/N, I′ = β S I/N − γ I, R′ = γ I, y la suma S + I + R permanece constante: nadie entra ni sale de la población.

El parámetro decisivo es el número reproductivo basal R₀ = β/γ: el número medio de contagios que provoca un infectado en una población totalmente susceptible. La epidemia sólo despega si I′ > 0 al principio, es decir si R₀ · S/N > 1. La cantidad Rₑ = R₀ · S/N es el número reproductivo efectivo, y desciende a medida que se agotan los susceptibles. Cuando Rₑ cruza el valor 1 se alcanza el pico de infectados, que ocurre precisamente cuando S = γN/β = N/R₀; a partir de ahí la curva baja aunque no se haga nada más.

De aquí sale el umbral de inmunidad de grupo: si una fracción p = 1 − 1/R₀ de la población es inmune desde el inicio, Rₑ nunca supera 1 y el brote no arranca. La epidemia, además, nunca infecta a todo el mundo: el tamaño final R(∞) satisface la ecuación implícita R(∞)/N = 1 − e^(−R₀ R(∞)/N), que se resuelve iterando. El sistema no tiene solución cerrada, así que se integra numéricamente con Runge-Kutta; ampliándolo con compartimentos de exposición, vacunación o pérdida de inmunidad se obtienen los modelos SEIR y SIRS que se usan en salud pública.

Fórmulas clave

• Sistema SIRS′ = −β S I / N, I′ = β S I / N − γ I, R′ = γ I
• Número reproductivo basalR₀ = β / γ
• Número reproductivo efectivoRₑ = R₀ · S/N
• Condición de broteI crece ⇔ S/N > 1/R₀
• Inmunidad de grupop = 1 − 1/R₀
• Pico de infectadosS en el pico = N/R₀; I máx / N = 1 − (1 + ln R₀)/R₀
• Tamaño finalR(∞)/N = 1 − e^(−R₀ R(∞)/N)

Dónde se aplica

Durante la pandemia de COVID-19 las agencias de salud publicaron a diario estimaciones de Rₑ calculadas con variantes de este modelo para decidir cuándo endurecer o relajar medidas; «achatar la curva» es literalmente reducir β para bajar el pico I máx aunque la epidemia dure más. El umbral p = 1 − 1/R₀ fija las coberturas de vacunación objetivo: con el sarampión, cuyo R₀ ronda 15, hace falta inmunizar a más del 93 % de la población. Las mismas ecuaciones se aplican a la difusión de rumores y de malware en redes, donde los «recuperados» son los nodos parcheados o los usuarios que dejan de compartir.

Con programación (Python)

import math

N, beta, gamma = 1000.0, 0.4, 0.1
R0 = beta / gamma

def f(t, s):
    S, I, R = s
    nuevos = beta * S * I / N
    return [-nuevos, nuevos - gamma * I, gamma * I]

def rk4(f, t, s, h):
    a = f(t, s)
    b = f(t + h / 2, [s[i] + h / 2 * a[i] for i in range(3)])
    c = f(t + h / 2, [s[i] + h / 2 * b[i] for i in range(3)])
    d = f(t + h, [s[i] + h * c[i] for i in range(3)])
    return [s[i] + h * (a[i] + 2 * b[i] + 2 * c[i] + d[i]) / 6 for i in range(3)]

print("R0 =", R0, " inmunidad de grupo =", round(100 * (1 - 1 / R0), 2), "%")
s, t, h, pico = [N - 1, 1.0, 0.0], 0.0, 0.25, (0.0, 0.0)
for n in range(1, 801):
    s = rk4(f, t, s, h)
    t += h
    if s[1] > pico[1]:
        pico = (t, s[1])
print("Pico: día", round(pico[0], 2), "con", round(pico[1], 2), "infectados")
print("I máx teórico:", round(N * (1 - (1 + math.log(R0)) / R0), 2))
print("Al final: S =", round(s[0], 2), " R =", round(s[2], 2))
z = 0.5
for _ in range(200):
    z = 1 - math.exp(-R0 * z)
print("Tamaño final teórico:", round(z * N, 2))

La versión en JavaScript, editable y ejecutable en el navegador, está en la página interactiva del tema.

Temas relacionados

Crecimiento logístico · Depredador y presa (Lotka-Volterra) · Crecimiento y decaimiento exponencial