1 Introducción

Este documento presenta un modelo SEIR (Susceptible-Expuesto-Infeccioso-Recuperado) que incorpora:

  • Demografía: Natalidad y mortalidad (μ ≠ 0)
  • Pérdida de inmunidad: Individuos recuperados regresan al estado susceptible (R → S)
  • Oscilaciones endémicas:

1.1 Ecuaciones del Modelo

Las ecuaciones diferenciales del modelo son:

\[\frac{dS}{dt} = \mu N - \frac{\beta S I}{N} - \mu S + \omega R\]

\[\frac{dE}{dt} = \frac{\beta S I}{N} - \sigma E - \mu E\]

\[\frac{dI}{dt} = \sigma E - \gamma I - \mu I\]

\[\frac{dR}{dt} = \gamma I - \mu R - \omega R\]

Donde:

  • \(\beta\) = tasa de transmisión
  • \(\sigma\) = tasa de incubación (1/período latente)
  • \(\gamma\) = tasa de recuperación (1/período infeccioso)
  • \(\mu\) = tasa de natalidad/mortalidad
  • \(\omega\) = tasa de pérdida de inmunidad
  • \(N\) = población total

2 Configuración del Modelo

2.1 Cargar Paquetes

# Instalar paquetes si es necesario
# install.packages(c("deSolve", "ggplot2", "dplyr", "gt", "pracma"))

library(deSolve)
library(ggplot2)
library(dplyr)
library(gt)
library(pracma)

2.2 Definir el Modelo

# Función del modelo SEIR con demografía y pérdida de inmunidad
seir_demographics <- function(time, state, parameters) {
  with(as.list(c(state, parameters)), {
    N <- S + E + I + R
    
    # Ecuaciones diferenciales
    dS <- mu * N - beta * S * I / N - mu * S + omega * R
    dE <- beta * S * I / N - sigma * E - mu * E
    dI <- sigma * E - gamma * I - mu * I
    dR <- gamma * I - mu * R - omega * R
    
    return(list(c(dS, dE, dI, dR)))
  })
}

2.3 Parámetros del Modelo

# Parámetros sarampión
parameters <- c(
  beta = 2.3,        # Tasa de transmisión (contactos/año)
  sigma = 1/7.5,     # Tasa de incubación (1/período latente en días)
  gamma = 1/6.5,     # Tasa de recuperación (1/período infeccioso en días)
  mu = 1/(75*365),         # Tasa de natalidad/mortalidad (1/duración vida)
  omega = 1/180       # Tasa de pérdida de inmunidad (1/duración inmunidad)
)

# Calcular R0
R0 <- parameters["beta"] / (parameters["gamma"] + parameters["mu"])

cat("Número básico de reproducción (R₀):", round(R0, 2), "\n")
## Número básico de reproducción (R₀): 14.95

2.4 Condiciones Iniciales

# Población de 1 millón
N <- 1e6

# Condiciones iniciales (cerca del equilibrio endémico)
initial_state <- c(
  S = 0.198 * N,     # 14% susceptibles
  E = 0.001 * N,    # 0.1% expuestos
  I = 0.001 * N,    # 0.1% infecciosos
  R = 0.8 * N     # 85.8% recuperados
)

# Verificar que suman N
cat("Población total inicial:", sum(initial_state), "\n")
## Población total inicial: 1e+06

3 Simulación del Modelo

# Tiempo de simulación: 
times <- seq(0, 365, by = 0.5)

# Resolver el sistema de ecuaciones diferenciales
output <- ode(
  y = initial_state, 
  times = times, 
  func = seir_demographics, 
  parms = parameters
)

# Convertir a data frame
output_df <- as.data.frame(output)

# Calcular proporciones
output_df <- output_df %>%
  mutate(
    S_prop = S / N,
    E_prop = E / N,
    I_prop = I / N,
    R_prop = R / N,
    N_total = S + E + I + R
  )

# Mostrar primeras filas
head(output_df)
##   time        S        E         I        R    S_prop      E_prop       I_prop
## 1  0.0 198000.0 1000.000 1000.0000 800000.0 0.1980000 0.001000000 0.0010000000
## 2  0.5 200005.8 1156.158  995.2348 797842.9 0.2000058 0.001156158 0.0009952348
## 3  1.0 202003.1 1304.566 1000.5752 795691.8 0.2020031 0.001304566 0.0010005752
## 4  1.5 203989.8 1447.860 1014.8628 793547.4 0.2039898 0.001447860 0.0010148628
## 5  2.0 205964.0 1588.330 1037.1843 791410.5 0.2059640 0.001588330 0.0010371843
## 6  2.5 207923.6 1727.986 1066.8296 789281.5 0.2079236 0.001727986 0.0010668296
##      R_prop N_total
## 1 0.8000000   1e+06
## 2 0.7978429   1e+06
## 3 0.7956918   1e+06
## 4 0.7935474   1e+06
## 5 0.7914105   1e+06
## 6 0.7892815   1e+06

4 Visualización de Resultados

4.1 Dinámica de Todos los Compartimentos

ggplot(output_df, aes(x = time)) +
  geom_line(aes(y = S_prop, color = "Susceptibles"), linewidth = 1.2) +
  geom_line(aes(y = E_prop, color = "Expuestos"), linewidth = 1.2) +
  geom_line(aes(y = I_prop, color = "Infecciosos"), linewidth = 1.2) +
  geom_line(aes(y = R_prop, color = "Recuperados"), linewidth = 1.2) +
  scale_color_manual(values = c(
    "Susceptibles" = "#3498db",
    "Expuestos" = "#f39c12",
    "Infecciosos" = "#e74c3c",
    "Recuperados" = "#2ecc71"
  )) +
  labs(
    title = "Modelo SEIR con Demografía y Pérdida de Inmunidad",
    subtitle = paste0("R₀ = ", round(R0, 2)),
    x = "Tiempo (días)",
    y = "Proporción de la población",
    color = "Compartimento"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    legend.position = "bottom",
    plot.title = element_text(face = "bold", size = 16),
    plot.subtitle = element_text(size = 12),
    panel.grid.minor = element_blank(),
    legend.title = element_text(face = "bold")
  )

4.2 Enfoque en Infecciosos

ggplot(output_df, aes(x = time, y = I_prop)) +
  geom_line(color = "#e74c3c", linewidth = 1.5) +
  labs(
    title = "Oscilaciones Endémicas en el Compartimento Infeccioso",
    x = "Tiempo (años)",
    y = "Proporción de infecciosos (I)"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    plot.title = element_text(face = "bold", size = 16),
    panel.grid.minor = element_blank()
  )

5 Análisis de Periodicidad

5.1 Detección de Picos

# Encontrar picos locales en el compartimento infeccioso
picos <- findpeaks(
  output_df$I_prop, 
  minpeakheight = max(output_df$I_prop) * 0.2
)

if (!is.null(picos)) {
  indices_picos <- picos[, 2]
  tiempos_picos <- output_df$time[indices_picos]
  
  # Calcular períodos entre picos
  if (length(tiempos_picos) > 1) {
    periodos <- diff(tiempos_picos)
    
    cat("=== ANÁLISIS DE PERIODICIDAD ===\n")
    cat("Número de picos detectados:", length(tiempos_picos), "\n")
    cat("Tiempos de los picos (días):\n")
    print(round(tiempos_picos, 2))
    cat("\nPeríodos entre picos (días):\n")
    print(round(periodos, 2))
    cat("\nPeríodo promedio:", round(mean(periodos), 2), "días\n")
    cat("Desviación estándar:", round(sd(periodos), 2), "días\n")
  }
}
## === ANÁLISIS DE PERIODICIDAD ===
## Número de picos detectados: 4 
## Tiempos de los picos (días):
## [1]  46 141 236 331
## 
## Períodos entre picos (días):
## [1] 95 95 95
## 
## Período promedio: 95 días
## Desviación estándar: 0 días

5.2 Visualización de Picos

ggplot(output_df, aes(x = time, y = I_prop)) +
  geom_line(color = "#e74c3c", linewidth = 1.2) +
  geom_point(
    data = output_df[indices_picos, ], 
    aes(x = time, y = I_prop), 
    color = "#c0392b", 
    size = 4, 
    shape = 21, 
    fill = "yellow",
    stroke = 2
  ) +
  labs(
    title = "Oscilaciones Endémicas con Picos Marcados",
    subtitle = paste0("Período promedio (días): ", round(mean(periodos), 2)),
    x = "Tiempo (días)",
    y = "Proporción de infecciosos (I)"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    plot.title = element_text(face = "bold", size = 16),
    plot.subtitle = element_text(size = 12),
    panel.grid.minor = element_blank()
  )