Actividad - COVID

Modelando la epidemia COVID

Este problema tiene como objetivo guiarte en el desarrollo y análisis de un modelo dinámico empleando el lenguaje de programación R. El objetivo es que te familiarices con las herramientas de análisis en R y que inicies el desarrollo de tus propios modelos dinámicos. Para este efecto tomaremos como referencia el caso de la epidemia COVID-19, basándonos en el modelo básico SI, descrito por Sterman (2013) en el apartado 9.2.1.

Descripción del caso:

Para iniciar tomaremos como referencia el diagrama “stock-flow” de Sterman (2013) mostrado en la siguiente figura.

Como primer paso haremos un listado del tipo de variables en el modelo:

Variables de Estado (Stock variables)

Population Susceptible to COVID
Population Infected with COVID

Variables de flujo (flow variables)

Infection Rate

Variables auxiliares endógenas (endogenous auxiliary variables)

Susceptible Contatcs
Probability of Contact with Infected Person
Contacts Between Infected and Uninfected People

Parámetros de simulación (variables en la frontera del sistema o exogenous auxiliary variables)

Contact Frequency
Total Population
Infectivity

Preguntas del caso:

4.1. Incluye tu modelo en un sólo chunk de código en el que se utilice la función plot para ver el comportamiento de las variables de estado.

library(deSolve)
parameters<-c(Infectivity = 0.1, # [1] dimmensionless
              Contact.Frequency = 2, # people/day
              Total.Population = 350 ) #people)

InitialConditions <- c(population.susceptible.to.COVID = 349 ,
                       population.infected.with.COVID = 1)

times <- seq(0 , #initial time, days
             120 , #end time, days
             0.25 ) #time step, days #porcentaje a la que va a correr del dia, en horas

intg.method<-c("rk4")

covid.epidemic <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    #Endogenous auxiliary variables
    Probability.of.Contact.with.Infected<-population.infected.with.COVID/Total.Population #dimensionless
    Susceptible.Contacts<-population.susceptible.to.COVID*Contact.Frequency #[people/time]
    Contacts.bt.Infected.and.Uninfected.People<-Susceptible.Contacts*Probability.of.Contact.with.Infected #[people/time]
      
    #Flow variables
    Infection.Rate<-Infectivity*Contacts.bt.Infected.and.Uninfected.People #[people/time]
      
    #State (stock) variables
    dpopulation.susceptible.to.COVID<-(-1)*Infection.Rate #Stock units: People/time
    dpopulation.infected.with.COVID<-Infection.Rate #Stock units: People/time
    
    list(c(dpopulation.susceptible.to.COVID,
           dpopulation.infected.with.COVID))
  })
}

out <- ode(y = InitialConditions,
           times = times,
           func = covid.epidemic,
           parms = parameters,
           method =intg.method )

plot(out,
     col=c("blue"))

4.2. ¿Qué sucede cuando inicializas la variable de estado “Population Infected with COVID” en cero? Explica brevemente que origina el comportamiento que observas, emplea la estructura del modelo para cimentar tu argumentación.

parameters<-c(Infectivity = 0.1, # [1] dimmensionless
              Contact.Frequency = 2, # people/day
              Total.Population = 350 ) #people)

InitialConditions <- c(population.susceptible.to.COVID = 349 ,
                       population.infected.with.COVID = 0)

times <- seq(0 , #initial time, days
             120 , #end time, days
             0.25 ) #time step, days #porcentaje a la que va a correr del dia, en horas

intg.method<-c("rk4")

covid.epidemic <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    #Endogenous auxiliary variables
    Probability.of.Contact.with.Infected<-population.infected.with.COVID/Total.Population #dimensionless
    Susceptible.Contacts<-population.susceptible.to.COVID*Contact.Frequency #[people/time]
    Contacts.bt.Infected.and.Uninfected.People<-Susceptible.Contacts*Probability.of.Contact.with.Infected #[people/time]
      
    #Flow variables
    Infection.Rate<-Infectivity*Contacts.bt.Infected.and.Uninfected.People #[people/time]
      
    #State (stock) variables
    dpopulation.susceptible.to.COVID<-(-1)*Infection.Rate #Stock units: People/time
    dpopulation.infected.with.COVID<-Infection.Rate #Stock units: People/time
    
    list(c(dpopulation.susceptible.to.COVID,
           dpopulation.infected.with.COVID))
  })
}

out <- ode(y = InitialConditions,
           times = times,
           func = covid.epidemic,
           parms = parameters,
           method =intg.method )

plot(out,
     col=c("blue"))

Respuesta: Si la variable de estado “Population Infected with COVID” se inicializa en cero, la propagación del virus no se activará dentro del modelo. Esto se debe a que la tasa de infección depende directamente del número de infectados. Al no haber individuos contagiados desde el inicio, la probabilidad de transmisión será cero, lo que impide la propagación del virus en la población. Como resultado, la población susceptible permanecerá constante en 350, mientras que el número de infectados se mantendrá en cero a lo largo del tiempo.

4.3. ¿Cómo cambia la dinámica de comportamiento del modelo si inicializas esta variable de estado a un valor positivo diferente de cero?.

parameters<-c(Infectivity = 0.1, # [1] dimmensionless
              Contact.Frequency = 2, # people/day
              Total.Population = 350 ) #people)

InitialConditions <- c(population.susceptible.to.COVID = 349 ,
                       population.infected.with.COVID = 12)

times <- seq(0 , #initial time, days
             120 , #end time, days
             0.25 ) #time step, days #porcentaje a la que va a correr del dia, en horas

intg.method<-c("rk4")

covid.epidemic <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    #Endogenous auxiliary variables
    Probability.of.Contact.with.Infected<-population.infected.with.COVID/Total.Population #dimensionless
    Susceptible.Contacts<-population.susceptible.to.COVID*Contact.Frequency #[people/time]
    Contacts.bt.Infected.and.Uninfected.People<-Susceptible.Contacts*Probability.of.Contact.with.Infected #[people/time]
      
    #Flow variables
    Infection.Rate<-Infectivity*Contacts.bt.Infected.and.Uninfected.People #[people/time]
      
    #State (stock) variables
    dpopulation.susceptible.to.COVID<-(-1)*Infection.Rate #Stock units: People/time
    dpopulation.infected.with.COVID<-Infection.Rate #Stock units: People/time
    
    list(c(dpopulation.susceptible.to.COVID,
           dpopulation.infected.with.COVID))
  })
}

out <- ode(y = InitialConditions,
           times = times,
           func = covid.epidemic,
           parms = parameters,
           method =intg.method )

plot(out,
     col=c("blue"))

Respuesta: Si la variable de estado “Population Infected with COVID” se inicializa con un valor positivo distinto de cero, el modelo reflejará una propagación activa del virus. A medida que los infectados aumentan, la tasa de infección crecerá, lo que reducirá progresivamente la población susceptible. Con el tiempo, más personas se contagiarán, provocando una dinámica en la que los infectados aumentan mientras la cantidad de individuos aún susceptibles disminuye.

4.4. ¿Cómo cambia la dinámica del sistema si aumenta el valor del parámetro “Contact Frequency”? ¿El valor de este parámetro modifica el valor final de la variable de estado “Population Infected with COVID”? Explica porque sí o porque no haciendo referencia a la estructura del modelo y a los resultados de la simulación.

Respuesta: Un aumento en Contact Frequency incrementa la velocidad de propagación del virus, lo que resulta en un crecimiento más rápido de la población infectada. Sin embargo, no altera el valor final de Population Infected with COVID, ya que la epidemia sigue un ciclo en el que eventualmente todos los individuos susceptibles se infectan o recuperan. Lo que cambia es el tiempo en que se alcanza este equilibrio, reduciendo la duración total del brote

parameters<-c(Infectivity = 0.1, # [1] dimmensionless
              Contact.Frequency = 5, # people/day
              Total.Population = 350 ) #people)

InitialConditions <- c(population.susceptible.to.COVID = 349 ,
                       population.infected.with.COVID = 1)

times <- seq(0 , #initial time, days
             120 , #end time, days
             0.25 ) #time step, days #porcentaje a la que va a correr del dia, en horas

intg.method<-c("rk4")

covid.epidemic <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    #Endogenous auxiliary variables
    Probability.of.Contact.with.Infected<-population.infected.with.COVID/Total.Population #dimensionless
    Susceptible.Contacts<-population.susceptible.to.COVID*Contact.Frequency #[people/time]
    Contacts.bt.Infected.and.Uninfected.People<-Susceptible.Contacts*Probability.of.Contact.with.Infected #[people/time]
      
    #Flow variables
    Infection.Rate<-Infectivity*Contacts.bt.Infected.and.Uninfected.People #[people/time]
      
    #State (stock) variables
    dpopulation.susceptible.to.COVID<-(-1)*Infection.Rate #Stock units: People/time
    dpopulation.infected.with.COVID<-Infection.Rate #Stock units: People/time
    
    list(c(dpopulation.susceptible.to.COVID,
           dpopulation.infected.with.COVID))
  })
}

out <- ode(y = InitialConditions,
           times = times,
           func = covid.epidemic,
           parms = parameters,
           method =intg.method )

plot(out,
     col=c("blue"))

4.5. ¿Cómo cambia el comportamiento del modelo si la variable de flujo “Infection Rate” cambia? Sigue los siguientes lineamientos para dar tu respuesta: Responde a esta pregunta describiendo brevemente los cambios que identificas al cambiar el valor de esta variable. Emplea un par de gráficos de comportamiento del modelo para dar soporte a tu respuesta.

Respuesta: Un cambio en la Infection Rate afecta la velocidad a la que las personas pasan de susceptibles a infectadas, modificando así el ritmo de la epidemia. Sin embargo, en un modelo como este, donde no hay reinfección ni intervenciones externas, el resultado final sigue siendo el mismo: la mayoría de la población eventualmente se infecta o se recupera. Lo que varía es la rapidez con la que se alcanza este estado.

4.6. El modelo que has desarrollado siguiendo el tutorial anterior es demasiado simple. Brevemente critica la formulación y estructura del modelo y lista las suposiciones del modelo que consideras son irrealistas.

Crítica del modelo: El modelo es una representación simplificada de la propagación de una epidemia y omite varios factores clave que influyen en la dinámica real de contagio. No considera la posibilidad de inmunidad, ni distingue entre recuperación y fallecimientos. Además, asume que la población es completamente homogénea en términos de interacción y susceptibilidad, lo que no refleja la variabilidad real en la sociedad. Tampoco incorpora medidas de intervención, como cuarentenas o vacunación, que pueden alterar significativamente el curso de una epidemia.

Suposiciones irreales:

Contacto homogéneo: Se asume que cada individuo tiene la misma probabilidad de interactuar con cualquier otro, ignorando patrones sociales y restricciones de movilidad.
Población fija: No hay nacimientos, muertes ni migraciones, lo que no refleja la dinámica demográfica real.
Sin recuperación ni inmunidad: El modelo inicial no contempla que los infectados puedan recuperarse ni volverse inmunes, lo que sobreestima el impacto a largo plazo.
Tasa de contacto constante: Se considera que la frecuencia de contacto entre personas no cambia con el tiempo, lo cual es irreal, ya que las medidas de salud pública o el comportamiento individual pueden modificar la propagación del virus.
No se distingue entre individuos expuestos e infectados, lo que simplifica en exceso el proceso de transmisión.

En el punto anterior identificaste algunas suposiciones irrealistas. Las siguientes preguntas tienen como objetivo que explores que sucede cuando se expanda el modelo para atender sus limitaciones.

Hasta el momento hemos asumido que la población se mantiene infectada con el virus COVID de manera indefinida. En epidemiologia esto se conoce como el modelo SI (i.e. Susceptible-Infectious). El modelo SI es apropiado para representar enfermedades crónicas para las que no existe una cura. Sin embargo, en el caso de muchas enfermedades infecciosas, incluyendo COVID, viruela o influenza, las personas infectadas pueden recuperarse o en los casos más lamentables morir.

El siguiente diagrama stock-flow expande la estructura del modelo base para describir el proceso de recuperación de la población infectada con COVID. Esta expansión del modelo en epidemiología es conocida como el modelo SIR (i.e. la “R” indica “Recovery”, explicada en el apartado 9.2.2 del libro de Sterman (2013)). Sigue las instrucciones siguientes para expandir el modelo del tutorial.

El diagrama stock-flow muestra que debes agregar tres variables nuevas: una nueva variable de estado “Population Recovered from COVID”, una nueva variable de flujo “Recovery Rate” (por simplicidad no distinguiremos entre los pacientes que se recuperan y aquellos que mueren) y un nuevo parámetro “Average Duration of Infection”.

El parámetro “Average Duration of Infection” indica el tiempo promedio (i.e. en días) que una persona permanece infectada con el virus COVID. Los epidemiólogos estiman que la fase de infección del COVID tiene una duración promedio de 7 a 21 días. Emplea tu criterio para elegir el valor de este parámetro.

Existen muchas formas de modelar la variable de flujo “Recovery Rate” pero la especificación empleada con mayor frecuencia es la siguiente: Recovery Rate=Population Infected with COVID/Average Duration of Infectivity

Para implementar exitosamente esta estructura en el modelo es necesario que especifiques que esta nueva variable de flujo afecta también a la variable de estado existente “Population Infected with COVID” de la siguiente manera (i.e. sintaxis en R):

dpopulation.infected.with.COVID<-Infection.Rate - Recovery.Rate

También es necesario agregar nueva variable de estado “Population Recovered from COVID”, esto lo puedes lograr empleando la siguiente especificación:

dpopulation.recovered.from.COVID<- Recovery.Rate

Recuerda que al agregar una nueva variable de estado es necesario que indiques en el vector de condiciones iniciales el valor inicial de esta variable y también indicar en la última línea de código de la función “covid.epidemic” que esta nueva variable de estado será impresa por la simulación:

list(c(dpopulation.susceptible.to.COVID, dpopulation.infected.with.COVID, dpopulation.recovered.from.COVID))

Emplea esta nueva versión del modelo para responder a las siguientes preguntas:

parameters<-c(Infectivity = 0.1, # [1] dimmensionless
              Contact.Frequency = 2, # people/day
              Total.Population = 350,
              average.duration.of.infectivity = 7) #people)

InitialConditions <- c(population.susceptible.to.COVID = 349 ,
                       population.infected.with.COVID = 1,
                       population.recovered.from.COVID = 0)



times <- seq(0 , #initial time, days
             120 , #end time, days
             0.25 ) #time step, days #porcentaje a la que va a correr del dia, en horas

intg.method<-c("rk4")

covid.epidemic <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    #Endogenous auxiliary variables
    Probability.of.Contact.with.Infected<-population.infected.with.COVID/Total.Population #dimensionless
    Susceptible.Contacts<-population.susceptible.to.COVID*Contact.Frequency #[people/time]
    Contacts.bt.Infected.and.Uninfected.People<-Susceptible.Contacts*Probability.of.Contact.with.Infected #[people/time]
      
    #Flow variables
    Infection.Rate<-Infectivity*Contacts.bt.Infected.and.Uninfected.People #[people/time]
    Recovery.Rate= population.infected.with.COVID/average.duration.of.infectivity
    
      
    #State (stock) variables
    dpopulation.susceptible.to.COVID<-(-1)*Infection.Rate #Stock units: People/time
    dpopulation.infected.with.COVID<-Infection.Rate - Recovery.Rate #Stock units: People/time
    dpopulation.recovered.from.COVID<- Recovery.Rate
    
    
    list(c(dpopulation.susceptible.to.COVID, dpopulation.infected.with.COVID, dpopulation.recovered.from.COVID))
  })
}

out <- ode(y = InitialConditions,
           times = times,
           func = covid.epidemic,
           parms = parameters,
           method =intg.method )

plot(out,
     col=c("red"))

4.7. ¿De qué manera cambia el comportamiento de la epidemia una vez que agregas estas nuevas variables al modelo? Al agregar la recuperación, la epidemia sigue un patrón típico del modelo SIR: la población susceptible no llega a cero porque algunos infectados se recuperan antes de contagiar a todos. Como resultado, la curva de infectados alcanza un pico y luego disminuye, mientras que la población recuperada crece con el tiempo. Esto refleja una epidemia autolimitada en la que no toda la población se infecta.

4.8. ¿Describe gráficamente y con un breve texto el efecto en el sistema de cambios (i.e. incremento y decremento) de las siguientes variables: “contact frequency” y “infectivity”? Enfatiza en las diferencias que percibes con respecto del comportamiento del modelo base.

library(deSolve)

# Función del modelo epidémico
covid.epidemic <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    Probability.of.Contact.with.Infected <- population.infected.with.COVID / Total.Population
    Susceptible.Contacts <- population.susceptible.to.COVID * Contact.Frequency
    Contacts.bt.Infected.and.Uninfected.People <- Susceptible.Contacts * Probability.of.Contact.with.Infected
    
    Infection.Rate <- Infectivity * Contacts.bt.Infected.and.Uninfected.People
    Recovery.Rate <- population.infected.with.COVID / average.duration.of.infectivity
    
    dpopulation.susceptible.to.COVID <- -Infection.Rate
    dpopulation.infected.with.COVID <- Infection.Rate - Recovery.Rate
    dpopulation.recovered.from.COVID <- Recovery.Rate
    
    list(c(dpopulation.susceptible.to.COVID, dpopulation.infected.with.COVID, dpopulation.recovered.from.COVID))
  })
}

# Condiciones iniciales
times <- seq(0, 120, 0.25)
InitialConditions <- c(population.susceptible.to.COVID = 349, 
                       population.infected.with.COVID = 1, 
                       population.recovered.from.COVID = 0)

# Escenarios de Contact Frequency e Infectivity
scenarios <- list(
  "Base" = c(Infectivity = 0.1, Contact.Frequency = 2),
  "High Contact" = c(Infectivity = 0.1, Contact.Frequency = 4),
  "Low Contact" = c(Infectivity = 0.1, Contact.Frequency = 1),
  "High Infectivity" = c(Infectivity = 0.2, Contact.Frequency = 2),
  "Low Infectivity" = c(Infectivity = 0.05, Contact.Frequency = 2)
)

# Simulación de cada escenario
library(ggplot2)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

results <- list()

for (scenario in names(scenarios)) {
  parameters <- c(scenarios[[scenario]], Total.Population = 350, average.duration.of.infectivity = 7)
  out <- ode(y = InitialConditions, times = times, func = covid.epidemic, parms = parameters, method = "rk4")
  df <- as.data.frame(out)
  df$Scenario <- scenario
  results[[scenario]] <- df
}

data <- bind_rows(results)

ggplot(data, aes(x = time, y = population.infected.with.COVID, color = Scenario)) +
  geom_line(size = 1) +
  labs(title = "Impacto de Contact Frequency e Infectivity en la Epidemia",
       x = "Tiempo (días)",
       y = "Población Infectada",
       color = "Escenario") +
  theme_minimal()

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Análisis de los resultados:

Efecto del Contact Frequency:

Un mayor contacto acelera la propagación del virus, alcanzando un pico de infección más alto y en menos tiempo.
Un menor contacto reduce la velocidad de propagación, aplanando la curva y extendiendo la epidemia por más tiempo.

Efecto de la Infectividad:

Una mayor infectividad provoca un aumento más rápido en los contagios y un pico más pronunciado.
Una menor infectividad retrasa y reduce la cantidad de personas infectadas en cada momento.

Comparado con el modelo base, ambos factores alteran la velocidad y magnitud de la epidemia, pero la infectividad tiene un impacto más fuerte en la cantidad máxima de infectados, mientras que el contact frequency afecta más la rapidez del brote.

Actividad - COVID

Realizado por: Diego Durán Martínez

Modelando la epidemia COVID