knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)
library(deSolve)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.2     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.4     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(rootSolve)

1 Introducción a los Modelos Epidemiológicos Mecanicistas

En esta parte del curso, trabajaremos con modelos construidos a partir del conocimiento de las relaciones entre los elementos que componen un sistema biológico-epidemiológico. El modelo básico que vamos a estudiar se conoce como SIR (Susceptible-Infectado-Recuperado) y se construye mediante la formulación de ecuaciones diferenciales ordinarias (ODE). Estas ecuaciones se resuelven utilizando procedimientos de integración numérica, ya que se trata de sistemas de ecuaciones diferenciales no-lineales que no poseen soluciones analíticas sencillas.

El modelo SIR fue formulado por Kermack and McKendrick (1927) para entender el funcionamiento de un sistema epidemiológico para una infección que se extiende y luego disminuye.

1.1 Ecuaciones del Modelo SIR

El planteamiento del modelo SIR se puede expresar mediante el siguiente sistema de ecuaciones diferenciales:

\[ \frac{dS}{dt} = \mu N - \mu S - \beta \frac{IS}{N} \quad (1) \] \[ \frac{dI}{dt} = \beta \frac{IS}{N} - (\mu + \gamma)I \quad (2) \] \[ \frac{dR}{dt} = \gamma I - \mu R \quad (3) \]

En las ecuaciones (1) y (2), podemos reemplazar el valor de \(N\) con 1 cuando \(S\), \(I\), y \(R\) se expresan como fracción de la población total \(N\).

1.2 Descripción y Supuestos del Sistema de Ecuaciones

Supuestos del modelo SIR con dinámica demográfica:

  • Población total:
    • \(N(t)\) es el tamaño poblacional. En equilibrio demográfico, la tasa per cápita de natalidad es \(\mu\) y la tasa per cápita de mortalidad es \(\mu\), por lo que, en ausencia de enfermedad letal, \(dN/dt = μN − μN = 0\) y \(N = S + I + R\) permanece constante.
    • Se asume mortalidad natural igual en todas las clases (\(S, I, R\)) y mortalidad inducida por la infección igual a 0.
  • Composición por clases:
    • No hay inmunidad perinatal: todos los nacidos entran a la clase susceptible \(S\) (flujo \(\mu N → S\)).
    • No hay latencia explícita: los individuos que se infectan pasan directamente a la clase infecciosa \(I\) (no hay clase “\(E\)”).
    • Los recuperados pasan a \(R\) y son inmunes de por vida (no hay pérdida de inmunidad ni retorno a \(S\)).
  • Transmisión:
    • El término de incidencia es de tipo mezcla homogénea (mass action estándar): \(\beta (S I)/N\).
      • \(S/N\) es la fracción de susceptibles; se asume mezcla al azar entre susceptibles e infecciosos.
      • \(\beta\) es el parámetro de transmisión efectivo (tasa de contactos potencialmente infecciosos por unidad de tiempo multiplicada por la probabilidad de transmisión por contacto).
    • La infectividad de un individuo no cambia a lo largo de su periodo infeccioso \(\beta\) constante en el tiempo y entre individuos).
  • Duración de la infección y salida de \(I\):
    • Los individuos permanecen infecciosos un tiempo promedio \(1/\gamma\), tras el cual se recuperan y pasan a \(R\).
    • Se suele asumir \(\mu << \gamma\), de modo que la recuperación ocurre típicamente antes que la muerte natural; aun así, la salida total de \(I\) es \((\gamma + \mu)I\) si se incluye mortalidad natural.
  • Procesos estocásticos y estructura:
    • El modelo es determinista y de compartimentos, sin estructura por edad, por espacio ni por red de contactos.
    • Las probabilidades de recuperación y muerte natural son constantes en el tiempo (tasas constantes), lo que implica tiempos de permanencia con distribución exponencial en cada compartimento.
  • Cantidad clave: número reproductivo básico

    • \(R_0 = \beta / (\gamma + \mu)\)

    • Interpretación: número esperado de infecciones secundarias causadas por un caso típico en una población completamente susceptible, con salidas de \(I\) por recuperación y muerte natural.

    • Umbral epidémico: si \(R_0 > 1\), la infección puede invadir; si \(R_0 < 1\), decae.

    • Umbral de inmunidad de manada: proporción crítica de inmunizados \(p_c = 1 - 1/R_0\) para que \(R_e < 1\) y la infección no pueda mantenerse.

Fracción susceptible en equilibrio endémico (si \(R_0 > 1\))

  • \(S^* / N = 1 / R_0\)
  • Prevalencia endémica \(I^* > 0\) determinada por balance entre entrada por transmisión y salida por \(\gamma + \mu\).

Aclaraciones y límites del modelo

  • Mezcla homogénea: en poblaciones reales suele haber heterogeneidad (edad, contactos, espacio); esto puede alterar \(\beta\) y \(R_0\) efectivos.
  • Incidencia estándar vs. densidad-dependiente: aquí se usa \(\beta SI/N\) (frecuencial). En contextos con contactos que crecen con la densidad, puede usarse \(\beta SI\) (densidad-dependiente).
  • Sin latencia: para enfermedades con periodo de incubación no infeccioso, el modelo SEIR es más apropiado (introduce una clase expuesta \(E\) con tasa de progresión \(\sigma\)).
  • Inmunidad de por vida: si la inmunidad se pierde, usar SIRS (con tasa de pérdida \(\omega: R → S\)).
  • Mortalidad por infección: si existe, añadir término \(\alpha I\) en \(dN/dt\) y en \(dI/dt\) (salida adicional de \(I\)); entonces \(R_0 = \beta / (\gamma + \mu + \alpha\)).
  • Parámetros constantes: en brotes reales, \(\beta\) puede variar por estacionalidad, cambios de comportamiento o intervenciones; puede modelarse \(\beta(t)\).

Interpretación práctica de \(\beta\) y \(\gamma\)

  • \(\gamma ≈ 1/(duración\ media\ como\ infeccioso\)). Por ejemplo, si la fase infecciosa dura 4 días, \(\gamma ≈ 1/4\ día^{-1}\).
  • \(\beta\) puede estimarse con \(R_0\) y \(\gamma\): \(\beta \approx R_0 (\gamma + \mu)\). En poblaciones humanas, \(\mu\) suele ser muy pequeño en escalas de brote, y a veces se omite.

2 Integración Numérica del Modelo SIR (Epidemia Cerrada)

Para resolver el sistema de ecuaciones diferenciales e integrar los cambios en función del tiempo, usaremos el paquete deSolve de R. Para simplificar, comenzaremos con un modelo de epidemia cerrada, haciendo que \(\mu = 0\).

En primer lugar, hay que definir una función en R en la cual están incluidos: tiempo (t), vector de las variables de estado (y = S, I, R), y los parámetros del modelo (parms = \(\beta, \mu, \gamma, N\)). A esta función la llamaremos sirmod.

sirmod <- function(t, y, parms){
   # Pull state variables from y vector
   S <- y[1]
   I <- y[2]
   R <- y[3]

   # Pull parameter values from parms vector
   beta <- parms["beta"]
   mu <- parms["mu"]
   gamma <- parms["gamma"]
   N <- parms["N"]

   # Define equations
   dS <- mu * (N  - S)  - beta * S * I / N
   dI <- beta * S * I / N - (mu + gamma) * I
   dR <- gamma * I - mu * R

   res <- c(dS, dI, dR)
   # Return list of gradients
   list(res)
}

Vamos a usar de ejemplo la infección con el virus SARS-CoV-2, causante de la enfermedad COVID-19. Este virus se transmite directamente entre humanos, y la infección confiere inmunidad por un tiempo prolongado. La duración de la fase infecciosa es de aproximadamente 5 días, y el valor de \(R_0\) estimado para la cepa original (Wuhan) es de aproximadamente 2.5.

Definimos los parámetros, los valores iniciales de las variables de estado y la escala de tiempo y resolución. La escala de tiempo es de 60 días con resolución de 0.5 días (para suavizar las curvas de la gráfica final). Modelaremos usando la fracción de individuos en cada clase (estandarizamos con \(N = 1\)), sin nacimientos ni muertes (\(\mu = 0\)), un período infeccioso de 5 días (\(\gamma = 1/5\)), y una tasa de transmisión de 0.5 (\(\beta = 0.5\)), obtenida de \(R_0 = 2.5\). Asumimos una fracción de infectados iniciales de 0.1% de la población inicial, y el resto susceptible.

times <- seq(0, 60, by = 0.5)
parms <- c(mu = 0.0, N = 1, beta = 0.5, gamma = 0.2)
start <- c(S = 0.999, I = 0.001, R = 0)

A continuación se aplica la función ode para encontrar soluciones a las ecuaciones diferenciales ordinarias, y los resultados se convierten a un data.frame para usarlos con más facilidad.

out <- ode(y = start, times = times, func = sirmod, parms = parms)
out <- as.data.frame(out)
head(round(out, 3))
##   time     S     I     R
## 1  0.0 0.999 0.001 0.000
## 2  0.5 0.999 0.001 0.000
## 3  1.0 0.998 0.001 0.000
## 4  1.5 0.998 0.002 0.000
## 5  2.0 0.998 0.002 0.001
## 6  2.5 0.997 0.002 0.001

2.1 Visualización de la Dinámica del Modelo SIR

Usando plot base de R

par(mar = c(5, 4, 4, 2) + 0.1) # Reset margins to default for base plot
plot(x = out$time, y = out$S, ylab = "Fracción de N",
     xlab = "Días", type = "l", lwd = 3, col = "blue")
lines(x = out$time, y = out$I, lwd = 3, col = "red")
lines(x = out$time, y = out$R, lwd = 3, col = "green")
legend("right", legend = c("S", "I", "R"), col = c("blue", "red", "green"), lty = 1, lwd = 3)

FIGURA 1: Dinámica del Modelo SIR (Epidemia Cerrada) para COVID-19 durante 60 días. Se muestran las variaciones (fracción de la población) en susceptibles (S), infectados (I) y recuperados (R).

Usando tidyverse y ggplot2

Para facilitar el análisis y la graficación de los resultados, rearreglaremos los datos usando el paquete tidyverse.

new.out <- as.data.frame(out) %>%
  gather(key, value, -time)
head(new.out)
##   time key     value
## 1  0.0   S 0.9990000
## 2  0.5   S 0.9987301
## 3  1.0   S 0.9984166
## 4  1.5   S 0.9980526
## 5  2.0   S 0.9976310
## 6  2.5   S 0.9971414

En este nuevo objeto, new.out, la variable key contiene los grupos S, I, y R, y value el respectivo valor del grupo en ese tiempo (time). Ahora construiremos una gráfica usando el paquete ggplot2:

ggplot(new.out, aes(x = time, y = value, color = key)) +
  geom_line(linewidth = 1) +
  labs(x = "Tiempo (Días)", y = "Fracción de la población") +
  scale_color_manual(values = c("S" = "blue", "I" = "red", "R" = "green"),
                     labels = c("S" = "Susceptibles", "I" = "Infectados", "R" = "Recuperados")) +
  theme_minimal() +
  theme(legend.title = element_blank())

FIGURA 2: Dinámica del Modelo SIR para COVID-19 (Epidemia Cerrada) durante 60 días. Se muestran las variaciones (fracción de la población) en susceptibles (S), infectados (I) y recuperados (R).

2.2 Valores Clave de la Epidemia

Un dato importante es el valor máximo de infectados, cuándo ocurren y el valor umbral de \(S\). \(S\) se define como el valor de susceptibles en el pico de la infección, es decir, cuando \(I\) alcanza su valor máximo. Luego de que los susceptibles caen por debajo de este valor, la infección comienza a decaer.

# Valor máximo de I y tiempo usando un pipe
max_I_data <- new.out %>%
  filter(key == "I") %>%
  filter(value == max(value)) %>%
  mutate(maxI = round(value, 2)) %>%
  select(time, maxI)

# Valor de S_umbral
S_at_max_I <- out$S[which.max(out$I)]

# Tabla de resultados con gt
library(gt)
result_table <- data.frame(max_I_data, S_umbral = round(S_at_max_I, 3))
param_table <- gt(result_table) %>%
  tab_header(
    title = "Tabla 1. Parámetros Epidemiológicos"
  ) %>%
  fmt_number(
    columns = time,
    decimals = 1
  ) %>%
  fmt_number(
    columns = S_umbral,
    decimals = 3
  )
param_table
Tabla 1. Parámetros Epidemiológicos
time maxI S_umbral
24.5 0.23 0.393

EJERCICIO 1

Examinar los resultados al cambiar los parámetros y valores iniciales y relacionarlos a escenarios epidemiológicos hipotéticos.

  • ¿Cuál es el efecto de disminuir la cantidad inicial de susceptibles, \(S\), y con cuáles escenarios epidemiológicos podemos relacionarlos?

  • Igual con \(\beta\) y \(\gamma\). ¿Qué significado epidemiológico tiene cambiarlos?

3 \(R_0\), Razón Reproductiva Básica

Para patógenos que se transmiten directamente, \(R_0\) es, por definición, el número esperado de casos secundarios que surgen de un caso específico en una población completamente susceptible. Es uno de los valores más importantes en los estudios epidemiológicos y manejo de enfermedades infecciosas.

Algunas de las inferencias que se obtienen del cálculo de \(R_0\):

3.1 Valor Teórico de \(R_0\)

A partir de la ecuación (2), y asumiendo \(\mu = 0\):

\[ \frac{dI}{dt} = \beta \frac{IS}{N} - \gamma I \]

Si consideramos \(N=1\) (fracción de la población):

\[ \frac{dI}{dt} = \beta IS - \gamma I \]

Obtenemos:

\[ \frac{dI}{dt} = I(\beta S - \gamma) \]

Una solución para el caso en que no ocurra epidemia (\(\frac{dI}{dt}=0\)) es que:

\[ \beta S - \gamma = 0 \]

Y de aquí podemos obtener la fracción de susceptibles necesaria para que ocurra la epidemia (umbral), que llamaremos \(S_u\):

\[ S_u = \frac{\gamma}{\beta} \]

A esta razón \(\frac{\gamma}{\beta}\) se la conoce como tasa de remoción relativa, y tiene que ser suficientemente baja para permitir que la enfermedad progrese. Sin embargo, el valor que más se utiliza en epidemiología es el inverso, y es al que se conoce como:

\[ R_0 = \frac{\beta}{\gamma} \]

Sabiendo que el \(R\) efectivo es \(R_E = R_0 \cdot s\), podemos encontrar el valor umbral teórico de \(R\) para que proceda la infección, reemplazando \(R_0\) por \(\frac{\beta}{\gamma}\) y \(s\) por \(S_u\):

\[ R_E = R_0 \cdot S_u = \frac{\beta}{\gamma} \cdot \frac{\gamma}{\beta} = 1 \]

3.2 Análisis de \(R_0\)

En la siguiente gráfica podemos observar las variaciones de \(R_E\) y el punto en el cual la infección decae al alcanzarse el umbral \(S_u\).

# Calculate R0
R0_val <- parms["beta"] / (parms["gamma"] + parms["mu"])
cat("R0 calculado:", round(R0_val, 2), "\n")
## R0 calculado: 2.5
# Adjust margins to accommodate a second right axis
par(mar = c(5, 5, 2, 5))

# Plot state variables
plot(x = out$time, y = out$S, ylab = "Fracción de la población",
     xlab = "Tiempo", type = "l", col = "blue", xlim = c(0, 60))
lines(x = out$time, y = out$I, col = "red")
lines(x = out$time, y = out$R, col = "green")

# Add vertical line at turnover point (max I)
xx <- out$time[which.max(out$I)]
lines(c(xx, xx), c(1/R0_val, max(out$I)), lty = 3)
ss <- out$S[which.max(out$I)]
lines(c(0, xx), c(ss, ss), lty = 3)

# Prepare to superimpose 2nd plot
par(new = TRUE)

# Plot effective reproductive ratio (w/o axes)
plot(x = out$time, y = R0_val * out$S, type = "l", lty = 2,
     lwd = 2, col = "black", axes = FALSE, xlab = NA,
     ylab = NA, ylim = c(0, 4), xlim = c(0, 60))
lines(c(xx, 26), c(1, 1), lty = 3) # Line for RE = 1

# Add right-hand axis for RE
axis(side = 4)
mtext(side = 4, line = 3, expression(R[E]), cex = 1.2) # Adjusted line for mtext

# Add legend
legend("right", legend = c("S", "I", "R", expression(R[E])),
       lty = c(1, 1, 1, 2), col = c("blue", "red", "green", "black"),
       lwd = c(1, 1, 1, 2))

FIGURA 3: Análisis de \(R_0\) y \(R_E\). La línea negra discontinua representa el valor umbral de \(R_E = 1\). La línea vertical punteada indica el punto en el cual la infección alcanza su pico máximo (\(I_{max}\)) y la línea horizontal punteada indica el valor umbral de susceptibles (\(S_u\)).

4 Tamaño Final de la Epidemia

El “tamaño final de la epidemia” (a menudo denotado como \(f\) o \(R^*\)) se refiere a la proporción total de la población que ha sido infectada (y por lo tanto se ha recuperado o ha fallecido, si el modelo lo incluye) al final de una epidemia. Es decir, es la fracción de individuos que han pasado por la clase de infectados a lo largo de todo el curso de la epidemia. Una aproximación es \(f \approx 1 - e^{-R_0}\). Se puede obtener un valor más exacto de manera computacional.

Usando el paquete rootSolve para encontrar el estado de equilibrio (cuando la epidemia ha terminado, \(I \approx 0\)):

parms_eq <- c(mu = 0.0, N = 1, beta = 3.5, gamma = 1/2) # Using beta = 3.5 for a larger R0
equil <- runsteady(y = c(S = 1 - 1E-5, I = 1E-5, R = 0),
                   times = c(0, 1E5), func = sirmod, parms = parms_eq)
round(equil$y, 3)
##     S     I     R 
## 0.001 0.000 0.999

El valor de R en el equilibrio (equil$y["R"]) indica la fracción de la población que se recuperó (o fue infectada) al final de la epidemia.

4.1 Gráfica del Tamaño Final de la Epidemia vs. \(R_0\)

# Candidate values for R0 and beta
R0_values <- seq(0.1, 10, length = 100)
betas <- R0_values * (1/2) # Assuming gamma = 1/2

# Vector of NAs to be filled with numbers
f <- rep(NA, length(R0_values))

# Loop over i from 1 to length(R0_values)
for(i in seq_along(R0_values)){
     equil <- runsteady(y = c(S = 1 - 1E-5, I = 1E-5, R = 0),
                        times = c(0, 1E5), func = sirmod,
                        parms = c(mu = 0, N = 1, beta = betas[i], gamma = 1/2))
     f[i] <- equil$y["R"]
}

plot(R0_values, f, type = "l", xlab = expression(R[0]),
     ylab = "Tamaño Final de la Epidemia (f)",
     main = expression("Tamaño Final de la Epidemia vs. " * R[0]))
curve(1 - exp(-x), from = 1, to = 10, add = TRUE, col = "red", lty = 2)
legend("bottomright", legend = c("Modelo SIR", expression(1 - e^{-R[0]})),
       col = c("black", "red"), lty = c(1, 2))

FIGURA 4: Tamaño Final de la Epidemia vs. \(R_0\). La línea negra representa el tamaño final de la epidemia calculado mediante el modelo SIR, mientras que la línea roja discontinua representa la aproximación \(1 - e^{-R_0}\).

5 Epidemia Abierta (con reclutamiento de susceptibles)

Ahora consideraremos el caso en el que hay reclutamiento de nuevos susceptibles (\(\mu > 0\)). En este escenario, la dinámica de la epidemia puede mostrar varias ocurrencias o un estado endémico.

Definimos un período de tiempo más largo, una tasa de mortalidad/natalidad (\(\mu\)) y valores iniciales que pueden llevar a una dinámica diferente.

beta <- 2.5
gamma <- 1/5
times_open <- seq(0, 730, by = .1) # 2 years (365 days * 2)
parms_open <- c(mu = 0.002, N = 1, beta = beta, gamma = gamma) # mu = 0.001 per day
start_open <- c(S = 0.199, I = 0.001, R = 0.8) # Example initial conditions

out_open <- as.data.frame(ode(y = start_open, times = times_open,
                              func = sirmod, parms = parms_open))

par(mfrow = c(1, 1)) # Make room for single plot
plot(times_open, out_open$I, ylab = "Fracción de Infectados", xlab = "Tiempo (días)",
     type = "l", col = "red")
# I at 730 days print on graph at 700
text(550, 2*out_open$I[which(out_open$time == 730)],
     labels = round(out_open$I[which(out_open$time == 730)], 3), pos = 4)

FIGURA 5: Dinámica del Modelo SIR (Epidemia Abierta) durante 2 años. Se muestra la fracción de infectados (I) a lo largo del tiempo, con reclutamiento de susceptibles.

EJERCICIO 2. Diversidad de Enfermedades Infecciosas

EJERCICIO 3. Modelo SIR para Influenza en Humacao

El problema: Si se declara una epidemia de influenza en Humacao, es muy probable no contar con las camas suficientes para los enfermos. Hay que comprobarlo.

Solución: Vacunar a la población, ¿cuánto?

Datos a buscar:

Usando el modelo SIR, simular la epidemia de influenza en Humacao con los datos obtenidos. Luego, incorporar la vacunación (modificando la fracción inicial de susceptibles o la tasa de transmisión) para determinar el porcentaje de vacunación necesario para evitar el colapso hospitalario.

6 Estimación Empírica de \(R_0\)

En situaciones reales, \(R_0\) puede ser estimado empíricamente a partir de datos epidemiológicos. Una forma común es utilizando la tasa de crecimiento inicial de los casos durante la fase exponencial de una epidemia. Si \(r\) es la tasa de crecimiento exponencial de los casos y \(D\) es el período infeccioso promedio, entonces: \[ R_0 \approx 1 + rD \]
Otra forma es utilizando la proporción de susceptibles al inicio de la epidemia y el tamaño final de la epidemia: \[ R_0 = -\frac{\ln(1 - f)}{s_0} \] donde \(f\) es el tamaño final de la epidemia y \(s_0\) es la proporción de susceptibles al inicio.

6.1 Ejemplo de Estimación Empírica de \(R_0\)

Supongamos que durante la fase inicial de una epidemia, los casos crecieron a una tasa de \(r = 0.2\) por día, y el período infeccioso promedio es de \(D = 5\) días. Entonces: \[ R_0 \approx 1 + (0.2)(5) = 2 \]
Si al final de la epidemia, el 70% de la población ha sido infectada (\(f = 0.6\)) y al inicio el 100% era susceptible (\(s_0 = 1\)), entonces: \[ R_0 = -\frac{\ln(1 - 0.7)}{1} = -\ln(0.3) \approx 1.20 \]

7 Referencias