Sesión 09 - 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)

Variables de flujo (flow variables)

Variables auxiliares endógenas (endogenous auxiliary variables)

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

Esta clasificación es útil para organizar el espacio de trabajo en Rstudio.

Tutorial de modelado del caso en R:

Iniciaremos por definir este modelo dinámico como una función definida por el usuario siguiendo los siguientes pasos.

#Carga la librería deSolve empleando la función library() 
library("deSolve")

Declara el espacio de trabajo como una función definida por el usuario. En este caso sólo tienes que cambiar el nombre de la función.

La función tiene que llamarse como el problema que estamos trabajando (covid.epidemic). La función va a depender de tres cosas t (tiempo), state (variables de estado), parameters (variables endógenas).

covid.epidemic <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    #Endogenous auxiliary variables
      
    #Flow variables
      
    #State (stock) variables
      
    list(c())
  })
}

Ahora agregaremos las variables de estado y los flujos asociadas a ellas. Iniciaremos con la variable de estado “Population Susceptible to COVID” como se muestra en el chunk siguiente.

covid.epidemic <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    #Endogenous auxiliary variables
      
    #Flow variables
      
    #State (stock) variables
    dpopulation.susceptible.to.COVID<-
      
    list(c())
  })
}

#Esta función se llena de abajo para arriba, primero las variables de estado, luego las variables de flujo y al final las variables endógenas.

Nota que en R las variables no pueden ser nombradas con espacios. De manera que esta variable de estado es nombrada como: “population.susceptible.to.COVID” nota también que las variables de estado son precedidas por el símbolo diferencial “d” para indicar que esta es una variable de estado, la sintaxis completa es: “dpopulation.susceptible.to.COVID”. Después de este paso habrás creado exitosamente la primera variable de estado del modelo.

Una práctica muy recomendable es documentar tus modelos. En R puedes hacer esto empleando el símbolo “#” y escribiendo delante de éste una breve descripción de la variable que estas representando. Además de describir la variable que estas creando es muy recomendable escribir también las unidades de medición de tu variable. El siguiente chunk muestra un ejemplo:

covid.epidemic <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    #Endogenous auxiliary variables
      
    #Flow variables
      
    #State (stock) variables
    #The population susceptible to COVID is equal to the 
    #population susceptible prior to the onset of the disease
    #less all of those that have contracted it
    dpopulation.susceptible.to.COVID<-(-1)*Infection.Rate #Stock units: People/time
    
    #Esta linea indica que variables se imprimen como 
    #resultado de la integración del modelo, todas las
    #variables de estado deben estar listadas
    list(c())
  })
}

Como se muestra en el stock-flow diagram, la variable “Infection.Rate” está conectada a la variable de estado como un flujo de salida (i.e. outflow variable), esta estructura se modela como se muestra en el chunk anterior. Nota que esta variable de flujo afecta con signo negativo a la variable de estado. Este signo negativo específica a esta variable de flujo como un flujo de salida.

Podemos seguir un proceso similar para modelar y documentar la segunda variable de estado “population.infected.with.COVID” como se muestra en la siguiente figura. Nota que en este caso la variable de flujo “Infection.Rate” es modelada como un flujo de entrada (i.e. inflow variable) y por esta razón no se incluye un signo negativo y nota también que la variable de estado “population.infected.with.COVID” es especificada precedida del signo diferencial “d”.

covid.epidemic <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    #Endogenous auxiliary variables
      
    #Flow variables
      
    #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())
  })
}

Como se muestra en el stock-flow diagram, en este modelo solo existe una sola variable de flujo que está conectada a las variables de estado. Esta variable de flujo es determinada por dos variables auxiliares: “Contacts between Infected and Uninfected People” e “Infectivity”. Para este caso especificamos la variable de flujo “Infection.Rate” como la multiplicación simple de estas dos variables auxiliares tal como se muestra en la siguiente figura. Nota que al agregar esta nueva variable también hemos incluido la documentación que describe esta variable y sus unidades de medición.

covid.epidemic <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    #Endogenous auxiliary variables
      
    #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())
  })
}

Al definir la variable de flujo “Infection.Rate” hemos agregado dos nuevas variables que debemos especificar. La variable “Contacts between Infected and Uninfected People” es una variable auxiliar endógena que es determinada por la variable “Susceptible Contacts” y la variable “Probability of Contact with Infected Person” esta interacción es especificada de la siguiente manera:

covid.epidemic <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    #Endogenous auxiliary variables
    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())
  })
}
#Infectivity es hexógena.

Al definir esta nueva variable hemos agregado dos nuevas variables auxiliares endógenas que debemos definir. La primera variable “Susceptible Contacts” es determinada por la variable de estado “Population Susceptible to COVID” y por la variable “Contact Frequency” y es especificada como se muestra en la figura siguiente. Nota que en esta ocasión al emplear la variable de estado para definir otra variable endógena auxiliar no es necesario usar el símbolo diferencial antes de la variable de estado.

covid.epidemic <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    #Endogenous auxiliary variables
    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())
  })
}

La última variable endógena auxiliar por modelar es la variable “Probability of Contact with Infected Person”. De manera similar al caso anterior, esta variable es determinada por la variable de estado “population infected with COVID” dividida por la variable exógena auxiliar “Total Population”

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())
  })
}

Esto concluye la especificación de todos los elementos endógenos del modelo. El siguiente paso es especificar los valores de los parámetros (i.e. variables auxiliares exógenas), las condiciones iniciales de las variables de estado, el horizonte temporal de análisis y el método de integración para la simulación.

En nuestro modelo hemos especificado tres parámetros: “Contact Frequency”, “Total Population” e “Infectivity”.

En R podemos especificar los parámetros del modelo empleando un vector como se muestra en la siguiente figura. Nota que el nombre de los parámetros es idéntico al nombre empleando en la especificación descrita en los pasos anteriores. También nota que cada parámetro está asociado a un valor numérico único para el cual correremos el modelo de simulación. Finalmente nota que para cada parámetro se especifican sus unidades de medición.

parameters<-c(Infectivity = 0.1, # [1] dimmensionless
              Contact.Frequency = 2, # people/day
              Total.Population = 350 ) #people)
#todas estas son variables exógenas

De igual manera podemos emplear un vector en R para definir las condiciones iniciales de cada variable de estado, esto se muestra en la siguiente figura. Nuevamente el nombre de las variables de estado es idéntico al empleado en la especificación del modelo y cada variable de estado es inicializada a un valor único.

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

Ahora es necesario especificar el vector de tiempo que será usado para simular el modelo. Esta especificación se lleva a cabo como se muestra en la siguiente figura. Nota que para especificar este vector de tiempo empleamos la función seq(). Esta función crea una secuencia numérica y requiere tres parámetros: valor inicial, valor final e intervalo de crecimiento. En nuestro modelo dinámico estos tres parámetros representan el tiempo inicial de la simulación, el tiempo final de la simulación y la resolución temporal de la simulación. Nota que para cada parámetro hemos indicado la unidad de tiempo correspondiente.

times <- seq(0 , #initial time, days
             120 , #end time, days
             0.25 ) #time step, days

Aún no discutimos con detalle las propiedades de los diferentes métodos de integración que podemos emplear para simular nuestro modelo. Por lo pronto, para este ejercicio, elegiremos el método Runge-Kutta de Orden 4. El chunk siguiente muestra la forma de especificar este método de integración.

intg.method<-c("rk4")

El último paso en el proceso de especificación del modelo es elegir las variables que serán “impresas” por la simulación. Este es un paso muy importante ya que nos permite elegir las variables que deseamos analizar. La figura siguiente muestra la forma de especificar las variables a imprimir por el modelo. Nota que esto es especificado en la última línea de código de la función contiene nuestro modelo dinámico. Nota también que en este caso hemos elegido imprimir sólo las variables de estado y que ambas están precedidas por el signo diferencial “d” y son concatenadas empleando la función c().

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,#el mismo orden que tengas en las variables de estado, tiene que ser el mismo orden en el que deben de aparecer en las condiciones iniciales. 
           dpopulation.infected.with.COVID))
  })
}

Esto concluye la especificación del modelo de simulación. El siguiente paso es ejecutar este código y llevar a cabo la simulación. Corre los chunks que contienen la librería deSolve, los vectores que guardamos como parameters, InitialConditions, times e intg.method, y también corre el chunk de la función covid.epidemic.

Una vez realizado esto, en la consola de R ya están cargados tanto el modelo como todos los parámetros necesarios para llevar a cabo la simulación, pero aún no hemos generado datos de la simulación. Para hacer esto, crearemos una base de datos “out” que contiene los resultados de la simulación empleando la función “ode”. Nota que la función “ode” emplea como parámetros de entrada condiciones iniciales de las variables de estado del modelo, el vector de tiempo, la función covid.epidemic que describe nuestro modelo, las variables exógenas (i.e. parámetros) y el método de integración. Todos estos elementos los hemos definido ya en los pasos anteriores.

out <- ode(y = InitialConditions,
           times = times,
           func = covid.epidemic,
           parms = parameters,
           method =intg.method )
##Cambio de variables de estado a través del tiempo por 120 días cada 6 horas.

Si has hecho esto correctamente verás un nuevo objeto llamado “out” listado en el panel superior derecho.

El paso final es analizar gráficamente los resultados para esto emplea la función “plot”

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

Preguntas del caso

Pregunta 4.1

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")


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

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

times <- seq(0 , #initial time, days
             120 , #end time, days
             0.25 ) #time step, days
##Cambio de variables de estado a través del tiempo por 120 días cada 6 horas.

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,#el mismo orden que tengas en las variables de estado, tiene que ser el mismo orden en el que deben de aparecer en las condiciones iniciales. 
           dpopulation.infected.with.COVID))
  })
}

intg.method<-c("rk4")

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

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

EN

En el gráfico de población suceptible, vemos que al inicio toda la población (350). Antes del día 20 empezamos a ver una disminución drástica. Esto se debe a que más personas se están infectando de COVID. También vemos que el número de personas suceptible se mantiene constante hasta en niveles muy bajos, casi 0. Esto puede decirnos que para ese momento, la mayoría de la población esta infcetada.

En el gráfico de la población infectada por COVID vemos que al inicio había 0 personas infectadas, sin embargo, estas aumentan rápidamente. Alcanza su punto máximo a los 40 días y después de esto, el número de infectados se estabiliza en un nivel de casi 350 personas.

Pregunta 4.2

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.

library("deSolve")


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

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

times <- seq(0 , #initial time, days
             120 , #end time, days
             0.25 ) #time step, days
##Cambio de variables de estado a través del tiempo por 120 días cada 6 horas.

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,#el mismo orden que tengas en las variables de estado, tiene que ser el mismo orden en el que deben de aparecer en las condiciones iniciales. 
           dpopulation.infected.with.COVID))
  })
}

intg.method<-c("rk4")

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

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

Esto significa que no existe una fuente de entrada, es decir, una fuente de contagio en la población (no existen personas infectadas). Al comenzar con cero personas infectadas, la variable population.infected.with.COVID en el vector InitialConditions es igual a 0. La población susceptible al COVID (population.susceptible.to.COVID) permanecerá igual (349) porque no existe otra variable que la afecte y se quedará como una línea horizontal.

Pregunta 4.3

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?

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 = 2)

times <- seq(0 , #initial time, days
             120 , #end time, days
             0.25 ) #time step, days
##Cambio de variables de estado a través del tiempo por 120 días cada 6 horas.

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,#el mismo orden que tengas en las variables de estado, tiene que ser el mismo orden en el que deben de aparecer en las condiciones iniciales. 
           dpopulation.infected.with.COVID))
  })
}

intg.method<-c("rk4")

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

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

Si se inicializa variable de estado population.infected.with.COVID, sí se inicializan estas variables a un valor positivo diferente de cero, esta acumulará y subirá la cantidad de personas infectadas. ¿Qué quiere decir esto? Que al haber personas infectadas desde el inicio, se generan contactos entre individuos susceptibles al COVID y aquellos que ya estan contagiados, lo que acelera la tasa de transmisión, y aumentará el número de personas infectadas mucho más rápido.Como consecuencia, la población suceptible disminuirá porque ya pasan a la población de infectados. Comparando con el caso anterior con una inicialización de 0, donde no existen personas infectadas, la tasa de infección aumentará más rápido al principio si ya se inicia con personas infectadas.

Pregunta 4.4

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.

library("deSolve")

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

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

times <- seq(0 , #initial time, days
             120 , #end time, days
             0.25 ) #time step, days
##Cambio de variables de estado a través del tiempo por 120 días cada 6 horas.

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,#el mismo orden que tengas en las variables de estado, tiene que ser el mismo orden en el que deben de aparecer en las condiciones iniciales. 
           dpopulation.infected.with.COVID))
  })
}

intg.method<-c("rk4")

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

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

Para este ejemplo usamos la variable de “contact frequency” de 8 personas por día. Si se auemnta la variable de “contact frequency” la velocidad de la propagación de contagios será mayor, por lo que al final habrá mayor número de personas infectadas. El valor de este parámetro sí modifica el valor final de la variable de estado “Population Infected with COVID”. ¿Por que? Porque la velocidad lo es todo. Es decir, si hay más “contact frequency”, hay más contacto suceptible, que a su vez causa más contacto entre personas infectadas y no infectadas, causando una mayor tasa de infectados de manera más rápida. Eventualmente, la mayoría de población suceptible terminará contagiada.

Pregunta 4.5

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.

library("deSolve")
library("ggplot2")
library("tidyr")
library("dplyr")

# Modelo base
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
    dpopulation.susceptible.to.COVID <- (-1) * Infection.Rate
    dpopulation.infected.with.COVID <- Infection.Rate
    list(c(dpopulation.susceptible.to.COVID, dpopulation.infected.with.COVID))
  })
}

InitialConditions <- c(population.susceptible.to.COVID = 349, 
                       population.infected.with.COVID = 2)
times <- seq(0, 
             120,
             0.25)

intg.method <- "rk4"

Contact.Frequency <- 2
Total.Population <- 350

#Alta
parameters_high <- c(Infectivity = 0.2, Contact.Frequency = Contact.Frequency, Total.Population = Total.Population)
out_high <- ode(y = InitialConditions, times = times, func = covid.epidemic, parms = parameters_high, method = intg.method)
df_high <- as.data.frame(out_high) %>% pivot_longer(cols = starts_with("population"), names_to = "State", values_to = "Population") %>% mutate(Infectivity = "High (0.2)")

#Baja
parameters_low <- c(Infectivity = 0.05, Contact.Frequency = Contact.Frequency, Total.Population = Total.Population)
out_low <- ode(y = InitialConditions, times = times, func = covid.epidemic, parms = parameters_low, method = intg.method)
df_low <- as.data.frame(out_low) %>% pivot_longer(cols = starts_with("population"), names_to = "State", values_to = "Population") %>% mutate(Infectivity = "Low (0.05)")

df_combined <- rbind(df_low, df_high)

# Graficar
ggplot(df_combined, aes(x = time, y = Population, color = State)) +
  geom_line() +
  facet_wrap(~ Infectivity, ncol = 1) +
  scale_color_manual(values = c("population.infected.with.COVID" = "red",
                               "population.susceptible.to.COVID" = "blue")) +
  labs(title = "Impact of the Change in Infection Rate",
       x = "time",
       y = "population",
       color = "Estado") +
  theme_minimal() +
  theme(plot.title = element_text(face = "bold"))

El comportamiento del modelo cambia significativamente si la variable de flujo “Infection Rate” se modifica. La “Infection Rate” directamente controla la velocidad a la que la población susceptible se convierten en infectada. Si se aumenta la variable de “infection rate”, más rápido disminuirá el número de personas suceptibles y un aumento más ráido en el número de personas infectadas, es decir, alcanzará su pico en un tiempo más rápido. Si se disminuye esta variable de “infection rate” ocurrirá el efceto contrario a lo anteriormente mencionado. El punto de equilibrio también se verá afectado. si se diminuye la variable, este equilibrio tardará más y se extendera, y el caso contrario si se aumenta.

Pregunta 4.6

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.

#El modelo supone que las personas infectadas permanecen infectadas indefinidamente. En un escenario real, las personas pueden recuperarse (no se considera la tasa de recuperación) o morir (no se considera la tasa de mortalidad), incluirlas mejoraría la precisión del modelo. Además, no considera directamente las medidas de protección como el distanciamiento, cubrebocas, lavado de manos, etc, lo cuál representa un papel fundamental porque eso podría desacelerar el crecimiento de la tasa de contagios.

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))

Pregunta 4.7

4.7. ¿De qué manera cambia el comportamiento de la epidemia una vez que agregas estas nuevas variables al modelo?

library("deSolve")

#para las de estado
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

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 #Nueva ecuación
    
    list(c(dpopulation.susceptible.to.COVID,#el mismo orden que tengas en las variables de estado, tiene que ser el mismo orden en el que deben de aparecer en las condiciones iniciales. 
           dpopulation.infected.with.COVID, 
           dpopulation.recovered.from.COVID))
  })
}

parameters<-c(Infectivity = 0.1, # [1] dimmensionless
              Contact.Frequency = 2, # people/day
              Total.Population = 350, 
              Average.Duration.of.Infectivity=14) #people)
#todas estas son variables exógenas

intg.method<-c("rk4")

out <- ode(y = InitialConditions,
           times = times,
           func = covid.epidemic,
           parms = parameters,
           method =intg.method )
##Cambio de variables de estado a través del tiempo por 120 días cada 6 horas.

#Graficamos
plot(out,
     col=c("blue"))
par(mar = c(5, 4, 4, 8), xpd = TRUE)

matplot(out[, "time"], out[, -1], type = "l", lty = 1, col = c("blue", "red", "green"),
        xlab = "Time (days)", ylab = "Population", main = "COVID Epidemic")

legend("topright", inset = c(-0.7, 0), legend = c("Susceptible", "Infected", "Recovered"), 
       col = c("blue", "red", "green"), lty = 1, cex = 0.9, bty = "n")

Tenemos cuatro gráficas que nos describen el comportamiento de todas las variable que agregamos. Para la población suceptible a COVID vemos que este grupo diminuye de manera pronunciada a los 40-50 días y después su decenso va más lento. ¿Por qué vemos esto? Porque entre más personas se infecten y se recuperan, quedarían inmunes por lo que ya no podrían volver a contagiarse. Después tenemos a la población infectada, la cuál tiene un aumento notable a partir de los 20 días, alcanzando su pico máximo de personas contagiadas a los 50 días y a partir de ahí va en descenso.¿Por qué pasa esto? Porque ahora tenemos una nueva variable llama “Tasa de población recuperada”, que claramente disminuye las personas contagiadas de COVID eventualmente porque se curan o mueren. Por último, tenemos a la población que se recuperó del COVID, lo cual está ligado a lo que vimos en las otras gráficas, donde la población infectada disminuyó porque las personas empezaron a curarse o a morir. Ahora vemos que la población recuperada aumenta notablemente a partir de los 40 días para después seguir creciendo de forma sostenida y un poco más estable. Eventualmente la pandemia se estabiliza porque la población va recuperándose cada vez más, y la cantidad de personas contagiadas es mínima, al igual que la población suceptible.

Pregunta 4.8

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.

# Aumento de Contact Frequency y Infectivity
comportamiento_1 <- c(Infectivity = 0.5, # Aumentamos infectividad
                  Contact.Frequency = 3, # Aumentamos la frecuencia de contacto
                  Total.Population = 350,
                  Average.Duration.of.Infectivity = 14)

# Disminución de Contact Frequency y Infectivity
comportamiento_2 <- c(Infectivity = 0.05, # Reducimos infectividad
                  Contact.Frequency = 1, # Reducimos la frecuencia de contacto
                  Total.Population = 350,
                  Average.Duration.of.Infectivity = 14)


# Realizar las simulaciones con los parámetros modificados

# Simulación con aumento de Contact Frequency e Infectivity
out_1 <- ode(y = InitialConditions, 
             times = times, 
             func = covid.epidemic, 
             parms = comportamiento_1,
             method = "rk4")

# Simulación con disminución de Contact Frequency e Infectivity
out_2 <- ode(y = InitialConditions, 
             times = times, 
             func = covid.epidemic, 
             parms = comportamiento_2,
             method = "rk4")

# Graficar los resultados para comparar

# Gráfico con Contact Frequency y Infectivity aumentados
plot(out_1,
     col=c("blue"))
par(mar = c(5, 4, 4, 8), xpd = TRUE)
matplot(out_1[, "time"], out_1[, -1], type = "l", lty = 1, col = c("blue", "red", "green"), 
        xlab = "Time (days)", ylab = "Population", main = "Increased Contact Frequency and Infectivity", cex.main = 0.8)
legend("topright", inset = c(-0.7, 0), legend = c("Susceptible", "Infected", "Recovered"), col = c("blue", "red", "green"), lty = 1, cex = 0.9, bty = "n")
# Gráfico con Contact Frequency y Infectivity disminuidos
plot(out_2,
     col=c("blue"))

En el caso donde las variables de “contact frequency” e “infectivity” tienen un incremento, vemos que la población suceptible tiene un decremento notable a partir de los 10 días y posteriormente se mantiene. En cuanto a la población infectada,vemos un pico notable a los 10 días para después tener un decremento.Por último, en la población recuperada de COVID vemos un incremento rápido desde el inicio, para después seguir manteniéndose en un nivel alto. ¿Qué significa todo esto? Que al aumentar las variables de “contact frequency” e “infectivity”, se va a propagar el COVID más rápido, logrando más personas infectadas pero también más personas recuperadas o fallecidas (porque hay más población contagiada que puede cuidarse para recuperarse o fallecer por complicaciones). En el caso donde las variables de “contact frequency” e “infectivity” tienen un decremento, vemos que la población suceptible también lo tiene pero de manera lenta (formando una curva). Este mismo comportamiento se observa en la población infectada. Por último, la población recuperada presenta un incremento gradual. ¿Qué nos dice todo esto? Que al reducir las variables de “contact frequency” e “infectivity”, la pandemia y los contagios se desaceleran en comparación con el caso donde se incrementan, en el cual existe una explosión del virus.