1 . 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

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

1.1 . 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, manteniendo el template que hemos usado en otros modelos

covid.epidemic <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
  # La función 'with' junto con as.list() convierte los vectores 'state' 
  # y'parameters' en una lista e sto permite referirnos a cada variable
  # directamente por su nombre sin necesidad de índices
    
  
    #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())
  })
}

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

NOTA : El factor(−1) está allí para indicar que la población susceptible disminuye a medida que se producen

nuevas infecciones. En otras palabras, el flujo de infección (Infection.Rate)

resta personas del compartimento de susceptibles

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:

NOTA : Al ser una variable endógena, su valor depende de una variable interna del modelo

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

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: 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
      # Su valor depende de otra
      # la variable de estado no utiliza "d"
    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.

1.2 . Variables exógenas del modelo

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)

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 de la 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().

  • dpopulation.susceptible.to.COVID: Representa la tasa de disminución de la población susceptible (flujo de salida).

  • dpopulation.infected.with.COVID: Representa la tasa de aumento de la población infectada (flujo de entrada).

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

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 )

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

2 . Pregunta 1 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")

# Condiciones iniciales del modelo

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

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

#Función
covid.epidemic <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    
    #Endogenous auxiliary variables - sus valores dependen de otro
    
    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
    # (-1) indica que la población susceptible disminuye a medida que se producen nuevas infecciones
    dpopulation.susceptible.to.COVID<-(-1)*Infection.Rate #Stock units: People/time
    
    dpopulation.infected.with.COVID<-Infection.Rate #Stock units: People/time
    
    # 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.
    list(c(dpopulation.susceptible.to.COVID, 
           dpopulation.infected.with.COVID))
  })
}

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

# Modelo de integración
intg.method<-c("rk4")

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

#Grafica del comportamiento de las variables de estado
plot(out,
     col=c("blue"))

3 . Pregunta 2 del caso

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

# Condiciones iniciales del modelo

InitialConditions <- c(population.susceptible.to.COVID = 349 , 
                       population.infected.with.COVID = 0) #Iniciando la variable en 0
# Vectores de Tiempo
##Cambio de variables de estado a través del tiempo por 120 días cada 6 horas.
times <- seq(0 ,    #initial time, days
             120 ,  #end time, days
             0.25 ) #time step, days

#Función
covid.epidemic <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    
    #Endogenous auxiliary variables - sus valores dependen de otro
    
    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
    # (-1) indica que la población susceptible disminuye a medida que se producen nuevas infecciones
    dpopulation.susceptible.to.COVID<-(-1)*Infection.Rate #Stock units: People/time
    
    dpopulation.infected.with.COVID<-Infection.Rate #Stock units: People/time
    
    # 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.
    list(c(dpopulation.susceptible.to.COVID, 
           dpopulation.infected.with.COVID))
  })
}

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

# Modelo de integración
intg.method<-c("rk4")

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

#Grafica del comportamiento de las variables de estado
plot(out,
     col=c("red"))

Si se inicia con la variable de estado de population.infected.with.COVID en 0, el modelo no generará infecciones lo que significa que nunca hubo un ¨paciente cero* y no hay forma que se transmita dicho virus.

En otras palabras, la simulación no generará ninguna infección a lo largo del tiempo porque la tasa de infección depende directamente del número de personas infectadas. Por lo que, si no hay individuos infectados, tanto la probalididad de contacto con una persona infectada, como el contanto con personas infectadas y susceptibles sigue siendo cero. Eso significa que no existiría una población infectada inicial.

4 .Pregunta 3 del caso

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?.

La referencia a esta pregunta, se puede ver en la primera gráfica de la pregunta 1. Lo que pasa es que el modelo sí mostrará una pandemia, por lo que la dinámica empieza a cambiar dado que sí habrá por lo menos una persona infectada, lo que la población susceptible comenzará a contagiarse a través del contacto con la población infectada. Esto significa que la tasa de infección¨ tendrá un efecto positivo.

De hecho se puede observar como la dinámica de *la tasa de infección o **infection.rate* es la que determina la dinámica, ya que cuando existe un número positivo, la dinámica se activa.

5 .Pregunta 4 del caso

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

# Condiciones iniciales del modelo

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

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

#Función
covid.epidemic <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    
    #Endogenous auxiliary variables - sus valores dependen de otro
    
    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
    # (-1) indica que la población susceptible disminuye a medida que se producen nuevas infecciones
    dpopulation.susceptible.to.COVID<-(-1)*Infection.Rate #Stock units: People/time
    
    dpopulation.infected.with.COVID<-Infection.Rate #Stock units: People/time
    
    # 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.
    list(c(dpopulation.susceptible.to.COVID, 
           dpopulation.infected.with.COVID))
  })
}

Aumento en el parametro de Contact.Frecuency

# Variables Exogenas del modelo
parameters<-c(Infectivity = 0.1, # [1] dimmensionless
              Contact.Frequency = 10, # people/day
              #mayor al valor inicial
              Total.Population = 350 ) #people)

# Modelo de integración
intg.method<-c("rk4")

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

#Grafica del comportamiento de las variables de estado
plot(out,
     col=c("purple"))

En la gráfica se puede observar que el pico de contagios ocurre antes, con un cambio mayor al valor inicial. Para comprender esta dinámica es necesario mencionar que el Contact.Frecuency es el promedio de los contactos diarios que tiene la población susceptible, por lo que se usa para calcular cuántas veces una personas no infectada puede estar expuesta al virus.

La dinámica anterior, explica porque a mayor numero de Contact.Frecuency, hay mayores contactos. Lo que afecta a la velocidad a la que se propaga el virus, para ejemplificar lo siguiente podemos observar que pasa cuando el valor es menos a la cantidad inicial (2)

# Variables Exogenas del modelo
parameters<-c(Infectivity = 0.1, # [1] dimmensionless
              Contact.Frequency = 1, # people/day
              # menor al valor inicial
              Total.Population = 350 ) #people)

# Modelo de integración
intg.method<-c("rk4")

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

#Grafica del comportamiento de las variables de estado
plot(out,
     col=c("purple"))

En la siguiente gráfica se puede evidenciar como hay menos oportunidades de contagio

6 .Pregunta 5 del caso

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

# Condiciones iniciales del modelo

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

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

#Función
covid.epidemic <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    
    #Endogenous auxiliary variables - sus valores dependen de otro
    
    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
    # (-1) indica que la población susceptible disminuye a medida que se producen nuevas infecciones
    dpopulation.susceptible.to.COVID<-(-1)*Infection.Rate 
    
    #Stock units: People/time
    
    dpopulation.infected.with.COVID<-Infection.Rate #Stock units: People/time
    
    # 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.
    list(c(dpopulation.susceptible.to.COVID, 
           dpopulation.infected.with.COVID))
  })
}

# Variables Exogenas del modelo 
  # modelación de escenarios
parameters_fast<-c(Infectivity = 0.1, # [1] dimmensionless
              Contact.Frequency = 6, # people/day
              Total.Population = 350 ) #people)

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

# Modelo de integración
intg.method<-c("rk4")

# Simulaciones
out_fast <- ode(y = InitialConditions,
                times = times,
                func = covid.epidemic,
                parms = parameters_fast,
                method = intg.method)

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

# Gráficas comparativas
par(mfrow = c(2,1))  # dos gráficos verticales

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

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

Se puede observar que este cambio, determina la velocidad en la que se propaga una pandemia, la diferencia entre esta pregunta y la anterior, es que la variable de flujo de Infection.Rate contempla la variable del TIEMPO, por que la manera de afectarlo es através del parámetro de Contact.Frequency

Observa nuevamente la imagen:

La dinámica es la siguiente, si hay un cambio en Infection.Rate, eso significa que hay un cambio en Contacts.bt.Infected.and.Uninfected.People, por lo que hay un cambio en Susceptible.Contacts y esa variable está determinada por una variables de estado :

Susceptible.Contacts<-population.susceptible.to.COVID *Contact.Frequency #[people/time]

Observa que la variable de estado se está multiplicando con el parámetro de Contact.Frecuency, es por eso que se utiliza para cambiar el valor de Infection.Rate

7 .Pregunta 6 del caso

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.

En el modelo solo existen dos variables de estado, las personas susceptibles e infectadas. Sin embargo no hay una variable de recuperación o de salida para ninguna infección, tampoco se muestra la emergencia del caso; es decir, no considera la mortalidad o los casos graves. Tampoco se considera el tiempo de contagio o si la misma persona se puede volver a infectar.

Suposiciones Irrealistas - Que no sabes cuanto dura el contagio y se maneja como si fuera estático - Todos tienen la misma probabilidad de contagiarse, asuminedo que todas las personas son igual de susceptibles. - El modelo maneja a personas que no pueden volver a infectarse y eso no es del todo realista

8 . Expanción del modelo

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)

  • 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.

NOTA Para este parametro se toma en cuenta que la fase de infección toma de 7 a 21 días. Además diversos sitios como CDC USA indica que los sintomas se empiezan a sentir entre los 2 y 14 días, por lo que se escojerá un parametro de 10 días, contando los días de contagió y los días de los sintomas

# Parámetros del modelo
parameters <- c(
  Infectivity = 0.1,
  Contact.Frequency = 2,
  Total.Population = 350,
  #Nuevo parámetro 
  Average.Duration.of.Infection = 10
)

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

library(deSolve)

# Condiciones iniciales del modelo
InitialConditions <- c(
  population.susceptible.to.COVID = 349,
  population.infected.with.COVID = 1,
  # nueva condición inicial
  population.recovered.from.COVID = 1
)

# Tiempo de simulación
times <- seq(0, 120, 0.25)

# Función del modelo SIR
covid.epidemic <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    
    # Variables auxiliares
    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
    
    # Tasa de infección
    Infection.Rate <- Infectivity * Contacts.bt.Infected.and.Uninfected.People
    
    # Tasa de recuperación (nueva)
    Recovery.Rate <- population.infected.with.COVID / Average.Duration.of.Infection
    
    # Ecuaciones diferenciales (stocks)
    dpopulation.susceptible.to.COVID <- (-1) * Infection.Rate
    
    ## nuevas variables de estado
    dpopulation.infected.with.COVID <- Infection.Rate - Recovery.Rate
    dpopulation.recovered.from.COVID <- Recovery.Rate
    
    # Retorno del sistema
    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:

A continuación se muestra el código con los cambios en el modelo:

library(deSolve)

# Condiciones iniciales del modelo
InitialConditions <- c(
  population.susceptible.to.COVID = 349,
  population.infected.with.COVID = 1,
  population.recovered.from.COVID = 1
)

# Tiempo de simulación
times <- seq(0, 120, 0.25)

# Función del modelo SIR
covid.epidemic <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    
    # Variables auxiliares
    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
    
    # Tasa de infección
    Infection.Rate <- Infectivity * Contacts.bt.Infected.and.Uninfected.People
    
    # Tasa de recuperación (nueva)
    Recovery.Rate <- population.infected.with.COVID / Average.Duration.of.Infection
    
    # Ecuaciones diferenciales (stocks)
    dpopulation.susceptible.to.COVID <- (-1) * Infection.Rate
    dpopulation.infected.with.COVID <- Infection.Rate - Recovery.Rate
    dpopulation.recovered.from.COVID <- Recovery.Rate
    
    # Retorno del sistema
    list(c(dpopulation.susceptible.to.COVID,
           dpopulation.infected.with.COVID,
           dpopulation.recovered.from.COVID))
  })
}

# Parámetros del modelo
parameters <- c(
  Infectivity = 0.1,
  Contact.Frequency = 2,
  Total.Population = 350,
  Average.Duration.of.Infection = 10
)


# Modelo de integración
intg.method<-c("rk4")

# Simulación del modelo
out <- ode(
  y = InitialConditions,
  times = times,
  func = covid.epidemic,
  parms = parameters,
  method = intg.method
)

# Gráfica de los resultados
plot(out, 
     col = c("blue"))

9 Pregunta 7 del caso

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

Como se puede observar en las gráficas anteriores, antes de incorporar la variable de estado de population.recovered.from.COVID, así como sus variables y parámetros, el comportamiento que podiamos observar era un crecimiento exponencial para después permanecer constante. Sin embargo cuando se agrega estas nuevas variables vemos que tanto population.recovered.from.COVID como recovery.rate ofrecen una salida de la infección, lo cual limita y eventualmente detiene la propagación, es por eso que en las gráficas anteriores se pueden ver osilaciones de onda.

Aunque vemos que hay un pico, eso se debe a que la variable de estado population.infected.with.COVID sube al inicio, llegando a un pico de infecciones y de ahí disminuye cuando los infectados comienzan a recuperarse. Por eso ahora vemos que el el gráfico de population.susceptible.to.COVID va disminuyendo con el tiempo, un efecto contrario al que tiene population.recovered.from.COVID que va aumentando con respecto al tiempo.

10 . Pregunta 8 del Caso

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)

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

# Tiempo
times <- seq(0, 120, 0.25)

# Función del modelo SIR
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.Infection
    
    dpopulation.susceptible.to.COVID <- (-1) * 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))
  })
}

# 1. Modelo base
parameters_base <- c(Infectivity = 0.1, Contact.Frequency = 2, Total.Population = 350, Average.Duration.of.Infection = 14)
out_base <- ode(y = InitialConditions, times = times, func = covid.epidemic, parms = parameters_base, method = "rk4")

# 2. Mayor Contact Frequency
parameters_highCF <- c(Infectivity = 0.1, Contact.Frequency = 5, Total.Population = 350, Average.Duration.of.Infection = 14)
out_highCF <- ode(y = InitialConditions, times = times, func = covid.epidemic, parms = parameters_highCF, method = "rk4")

# 3. Mayor Infectivity
parameters_highInf <- c(Infectivity = 0.3, Contact.Frequency = 2, Total.Population = 350, Average.Duration.of.Infection = 14)
out_highInf <- ode(y = InitialConditions, times = times, func = covid.epidemic, parms = parameters_highInf, method = "rk4")

# Comparación gráfica

plot(out_base[, "time"], out_base[, "population.infected.with.COVID"], type = "l", col = "red", lwd = 2,
     main = "Modelo base – Contact.Frequency = 2, Infectivity = 0.1",
     xlab = "Días", ylab = "Población infectada")

plot(out_highCF[, "time"], out_highCF[, "population.infected.with.COVID"], type = "l", col = "orange", lwd = 2,
     main = "Mayor Contact.Frequency = 5 (Infectividad constante)",
     xlab = "Días", ylab = "Población infectada")

En esta gráfica se puede observar que a mayor contact.Frecuency, la virus se propaga más rápido, lo que hace que el pico se acelere y llegué antes, eso a su vez nos dice que hay una alta probabilidad de contagios al día. A diferencia de la gráfica que se presenta a continuación; el incremento de infectivity aunque tiene un comportamiento similar, esta gráfica nos indica que tan contagioso es cada encuentro, lo que muestra mayor riegos de transmición por interración.

plot(out_highInf[, "time"], out_highInf[, "population.infected.with.COVID"], type = "l", col = "purple", lwd = 2,
     main = "Mayor Infectivity = 0.3 (Frecuencia de contacto constante)",
     xlab = "Días", ylab = "Población infectada")

Se puede ver que ambos casos tiene comportamientos similares, y es que como se mencionaba en las preguntas anteriores, al cambiar contact.frecuency la dinámica se refleja en la aceleración de los contagios. Lo interesante de esto, es que el cambio tanto en contact.frecuency como en Infectivity afecta a infection rate es por eso que se ve ese comportamiento.