Licencias

El código está bajo una licencia Aladdin Free Public License la cual permite el uso gratuito del mismo y su reproducción verbatim sin posibilidad de venta. El modelo y este documento se encuentran bajo una licencia de creative commons attribution sharealike no comercial.

Cualquier uso del modelo o el código debe dar atribución a sus creadores.

Modelo

Figura 1

Se desarrolló un modelo multigrupo SEIRS para describir la propagación del COVID-19 en México, asumiendo distintos escenarios de cuarentena . El modelo presentado toma en cuenta la heterogeneidad del estado de infección, es decir, los grupos de variantes clínicos que pueden ocurrir cuando se contrae la enfermedad. También se considera que la población se puede dividir en m clases, que pueden referirse tanto a grupos de edad, clases socioeconómicas o la mezcla de las anteriores. Cada grupo \(k\), donde \(1≤ k ≤m\), es a su vez particionado en compartimientos donde \(S_k\) representa a los individuos susceptibles; \(E_k\) a los individuos en periodo latente; \(I_k = \sum_{j=0}^{k}\) representa a los individuos que han sido infectados, donde \(I_{0,k}\) son individuos asintómaticos, \(I_{1, k}\) individuos con síntomas leves, \(I_{2,k}\) individuos en estado grave e \(I_{3,k}\) individuos en estado crítico; \(M_{k}\) representa el número de muertos de individuos en estado crítico y \(R_k\) a los individuos recuperados. Las cuarentenas son divididas en dos tipos: la primera representa una cuarentena voluntaria, donde un porcentaje de la población susceptible (\(Q_{S,k}\)), de la población asintomática (\(Q_{0,k}\)) y un porcentaje de la población en estado latente (\(Q_{E,k}\)) entra en aislamiento; la segunda cuarentena se compone de los pacientes infectados con síntomas leves \(Q_{2,k}\) y se considera que sólo un pequeño porcentaje de la población entrará a este grupo. La Figura 1 describe el modelo de compartimientos utilizado.

Las variables del modelo son:

Variables

\[\begin{align} S & = \text{ Susceptibles (pueden enfermar)} \\ E & = \text{Expuestos (incubando)} \\ I_0 & = \text{Infectados asintomáticos} \\ I_1 & = \text{Infectados con síntomas leves} \\ I_2 & = \text{Infectados graves que requieren cama de hospital} \\ I_3 & = \text{Infectados críticos que requieren terapia intensiva} \\ Q & = \text{Susceptibles en cuarentena voluntaria} \\ Q_1 & =\text{Asintomáticos en cuarentena voluntaria} \\ Q_2 & = \text{Infectados con síntomas leves en cuarentena} \\ R &= \text{Recuperados} \\ M & = \text{Muertos} \end{align}\]

Ecuaciones

Si consideramos que la población es dividida en m grupos (grupos de edad, nivel socioeconómico o cualquier otra variable de interés):

\[\begin{eqnarray} \dfrac{dS_k}{dt} & = & \Gamma_{S,k} - {\alpha_k(t)\,} \sum_{\ell = 1}^{m} \bigg( \sum_{j = 0}^{3} \gamma_{j,\ell} I_{j,\ell} \bigg) \cdot S_k \nonumber \\ &&{- \sum_{i = 1}^\infty q_{i,k} \cdot S_k(\tau_{i,k}^\text{init}) \cdot \Big(\delta(t - \tau_{i,k}^\text{init}) - \delta(t - \tau_{i,k}^\text{final}) \Big), } \label{eq:Sk} \\ \dfrac{dE_k}{dt} & = & \Gamma_{E,k} + {\alpha_k(t)\,} \sum_{\ell = 1}^{m} \bigg( \sum_{j = 0}^{3} \gamma_{j,\ell} I_{j,\ell} \bigg) \cdot S_k - \gamma_{E,k} \cdot E_k \nonumber \\ & & - {\sum_{i = 1}^\infty \Big( q_{i,k} \cdot E_{k} \cdot \delta(t - \tau_{i,k}^\text{init}) - Q_{E,k} \cdot \delta(t - \tau_{i,k}^\text{final}) \Big),} \label{eq:Ek} \\ \dfrac{dI_{0,k}}{dt} & = & \Gamma_{0,k} + p_{I_0,k} \cdot \gamma_{E,k} \cdot E_k - \beta_{0 \to R, k} \cdot I_{0,k} \nonumber \\ & & - {\sum_{i = 1}^\infty \Big( q_{i,k} \cdot I_{0,k} \cdot \delta(t - \tau_{i,k}^\text{init}) - Q_{0,k} \cdot \delta(t - \tau_{i,k}^\text{final}) \Big),} \\ \dfrac{dI_{1,k}}{dt} & =& \Gamma_{1,k} + (1 - p_{I_0,k}) \!\cdot\! \gamma_{E,k} \!\cdot\! E_k - \Big( p_{I_{2,k}} \!\cdot\! \beta_{1 \to 2, k} + (1 - p_{I_{2,k}}) \!\cdot\! \beta_{1 \to R, k} + p_{1,k} \Big) \!\cdot\! I_{1,k}, \\ \dfrac{dI_{2,k}}{dt} & = & p_{I_{2,k}} \cdot \beta_{1 \to 2, k}\cdot \big( Q_{1,k} + I_{1,k} \big) \nonumber \\ && - \Big( p_{I_{3,k}} (I_2) \cdot \beta_{2\to 3,k}(I_2) + \big( 1 - p_{I_{3,k}}(I_2) \big) \cdot \beta_{2 \to R,k} (I_2)\Big) \cdot I_{2,k}, \\ \dfrac{dI_{3,k}}{dt} & = & p_{I_{3,k}}(I_2) \cdot \beta_{2\to 3,k} (I_2) \cdot I_{2,k} \nonumber \\ &&- \Big( \big(1 - p_{M,k}(I_3) \big) \cdot \beta_{3\to R, k} + p_{M,k}(I_3) \cdot \beta_{3 \to M, k} \Big) \cdot I_{3,k}, \\ \dfrac{dM_{k}}{dt} & = & p_{M,k}(I_3) \cdot \beta_{3\to M,{k}} \cdot I_{3,k}, \\ \dfrac{dQ_{S,k}}{dt} & = & {\sum_{i = 1}^\infty q_{i,k} \cdot S_k(\tau_{i,k}^\text{init}) \cdot \Big(\delta(t - \tau_{i,k}^\text{init}) - \delta(t - \tau_{i,k}^\text{final}) \Big) } \label{eq:QSk}, \\ \dfrac{dQ_{E,k}}{dt} & = & - \gamma_{E,k} \cdot Q_{E, k} + {\sum_{i = 1}^\infty \bigg( q_{i,k} \cdot E_{k} \cdot \delta(t - \tau_{i,k}^\text{init}) - Q_{E,k} \cdot \delta(t - \tau_{i,k}^\text{final}) \bigg),} \\ \dfrac{dQ_{0,k}}{dt} & = & p_{I_0,k} \cdot \gamma_{E,k} \cdot Q_{E,k} - \beta_{0 \rightarrow R,k}\cdot Q_{0, k} \nonumber \\ && + {\sum_{i = 1}^\infty \bigg( q_{i,k} \cdot I_{0,k} \cdot \delta(t - \tau_{i,k}^\text{init}) - Q_{0,k} \cdot \delta(t - \tau_{i,k}^\text{final}) \bigg),} \\ \dfrac{dQ_{1,k}}{dt} & = & (1 - p_{I_0,k}) \!\cdot\! \gamma_{E,k} \!\cdot\! Q_{E,k} + p_{1,k} \!\cdot\! I_{1,k} - \big( p_{I_{2,k}} \!\cdot\! \beta_{1 \rightarrow 2, k} + (1 - p_{I_{2,k}}) \!\cdot\! \beta_{1 \rightarrow R,k} \big) \!\cdot\! Q_{1,k}, \end{eqnarray}\] donde \(j= 0, 1, 2, 3\) representa el estado de la infección y \(k, l = 1, 2, ...., m\) el grupo de interés. \(\Gamma_{j,k}\) representa las tasas de recrutamiento o migración y \(1/\gamma_{E,k}\) al tiempo promedio en estado de latencia. Se asume que todos los individuos en el estado \(I_{j,k}\) se recuperarán a una tasa constante \(\beta_{j\rightarrow R, k}\) y que \(p_{I_{j},k}\) es la proporción de infectados en cada k-compartimiento de \(I_{j,k}\) que pasan al compartimiento \(I_{j+1,k}\). La fuerza de infección de \(I_{j,l}\) en cada grupo \(k\) está representada por \(\gamma_{j, \ell}\) y se define como: \[ \gamma_{j, \ell} = \lambda_{j, \ell} \sum_{k = 1}^{m} c_{\ell}^{k} I_{j,k}, \] donde \(j =0, 1, 2, 3\), \(k = 1, ...., m\), la tasa de contagio es representada para cada infección del tipo \(j\) y para cada grupo \(l\), con \(\lambda_{j,l}\) como el riesgo de infección del grupo \(l\). La matriz de \(m\times m\) \(C= (c_{l}^{k})\) es la matriz de contacto donde cada entrada \(c_{l}^{k}\) representa la tasa de contacto entre grupos \(l\) y \(k\).

Para evaluar las condiciones de escasez hospitalaria se define: \[ p_{I_{M,k}}(I_3) = \begin{cases} q_{M,k}^1, & \textrm{ if } I_3 < \hat{I}_3, \\ q_{M,k}^2, & \textrm{ if } I_3 \geq \hat{I}_3, \end{cases}\] que denota el porcentaje de personas que pueden ser tratadas si \(\hat I_3\) representa el número de casos necesarios para saturar las unidades de cuidado intensivo del sistema de salud. Similarmente, se tiene que

\[ p_{I_{2,k}}(I_2) = \begin{cases} q_{3,k}^1, & \textrm{ if } I_2 < \hat{I}_2, \\ q_{3,k}^2, & \textrm{ if } I_2 \geq \hat{I}_2, \end{cases}\] donde \(\hat I_{2}\) representa el número de casos de infectados que saturan el sistema de salud para camas con conexión de oxígeno. Finalmente, \[ \beta_{2 \to 3, k}(I_2) = \begin{cases} \beta_{2 \to 3, k}^1, & \textrm{ if } I_2 < \hat{I}_2, \\ \beta_{2 \to 3, k}^2, & \textrm{ if } I_2 \geq \hat{I}_2, \\ \end{cases} \] representa el cambio en la duración si el sistema de salud tiene conexión a oxígeno disponible.

En este modelo los parámetros representan:

Variable	Interpretación
\(p_{I_0}\)	Proporción de asintomáticos
\(p_{I_2}\)	Proporción de infecciosos tipo 1 que se vuelven tipo 2
\(p_{I_3}\)	Proporción de infecciosos tipo 2 que se vuelven tipo 3
\(p_{M}\)	Proporción de infecciosos tipo 3 que se mueren
\(\tau_1\)	Día de implementación de la cuarentena \(Q_1\)
\(\tau_2\)	Día de implementación de la cuarentena \(Q_2\)
\(\Delta \tau_1\)	Duración de la cuarentena \(Q_1\)
\(\Delta \tau_2\)	Duración de la cuarentena \(Q_2\)
\(\lambda_{I_1\to I_2}\)	Promedio de días en los que infecciosos tipo 1 se vuelven tipo 2
\(\lambda_{I_2\to I_3}\)	Promedio de días en los que infecciosos tipo 2 se vuelven tipo 3
\(\gamma_E\)	Tasa de incubación
\(\lambda_{Q_2 \to I_2}\)	Tasa con la que infectados tipo 1 en cuarentena pueden volverse infectados tipo 2
\(\lambda_{I_1\to R}, \lambda_{I_2\to R}, \lambda_{I_3\to R}\)	Promedio de días en los que infecciosos tipo 1 se vuelven tipo 2
\(\lambda_{I_3\to M}\)	Promedio de días antes de morir para infectados tipo 3
\(\lambda_{R \to S}\)	Si hay remisión, número promedio de días en la que los recuperados se vuelven susceptibles de nuevo.
\(\gamma_{1}, \gamma_{2}, \gamma_{3}\)	Tasa de contagio para cada uno de los estadíos de infección sintomática \(I_1, I_2, I_3\)
\(\gamma_A\)	Tasa de contagio para los asintomáticos
\(q_1\)	Proporción de individuos en cuarentena sin síntomas (susceptibles y asintomáticos)
\(q_2\)	Proporción de individuos en cuarentena con síntomas de infección tipo \(1\)

Cuarentenas

Para mitigar el contagio se implementaron 3 tipos de cuarentena en el modelo:

Cuarentena tipo 1. Son las cuarentenas del tipo voluntario, donde un porcentaje \(q_{i,k}\) de la población adopta la cuarentena de un tiempo inicial \(\tau_{i,k}^{inicial}\) a un tiempo \(\tau_{i,k}^{final}\). Este tipo de cuarentenas pueden ser implementadas de distinta forma para cada grupo \(k\), indicado por el índice \(i\), para capturar este fenómeno de migración entre compartimientos se utilizaron deltas de Dirac. Por ejemplo, la solución de la ecuación \(\frac{d Q_{S,k}}{dt}\) es:

\[ Q_{S,k} = Q_{S,k}(0) + \sum_{i = 1}^\infty q_{i,k} \cdot S_k(\tau_{i,k}^\text{init}) \Big(\mathcal{H}(t - \tau_{i,k}^\text{init}) - \mathcal{H}(t - \tau_{i,k}^\text{final}) \Big), \] donde \(\mathcal{H}\) es una función Heaviside:

\[ \mathcal{H}(t) = \begin{cases} 0, & \text{if } t < 0 \\ 1, & \text{if } t > 0 \end{cases}. \]

Cuarentena tipo 2. Representan la cuarentena voluntaria del porcentaje \(p_{1,k}\) individuos que presentan síntomas de la enfermedad.
Cuarentena tipo 3. También llamadas toque de queda. Son controladas para cada grupo de la siguiente manera: \[ \alpha_k(t) \;=\; \begin{cases} 1, & t \in [0,\,6) \mod{7} \\ q_{S,k}, & t \in [6,\,7) \mod{7} \end{cases}. \] Cada \(q_{S,k}\) puede variar en el tiempo al rotar la circulación de individuos. (2020, @arXiv_v1, @Brauer, @BrauerCastillo, @Ferguson)

Condiciones iniciales

Las condiciones iniciales para el modelo se definen como \(M_k(0) = Q_{S,k}(0) = Q_{1,k}(0) = Q_{2,k}(0) =0\) para \(k=1, ..., m\) al inicio de la enfermedad, pues aún no habrán cuarentenas o muertes como consecuencia de la enfermedad. Las variables \(I_{j,k}(0)\) también pueden considerarse nulas. La población total inicial se divide entre el total inicial de susceptibles \(S_{k}(0)\), expuestos \(E_{k}(0)\) e infectados asintomáticos \(I_{0,k}(0)\) y debe satisfacer:

\[ \sum_{k = 1}^m \Big( S_{k}(0) + E_{k}(0) + I_{0,k}(0) \Big) \;=\; 1, \]

Programación

rm(list = ls())

#Required libraries to run the model
library(deSolve)
library(tidyverse)

## ── Attaching packages ────────────────────────────────────────────────────── tidyverse 1.3.0 ──

## ✓ ggplot2 3.3.0     ✓ purrr   0.3.3
## ✓ tibble  3.0.0     ✓ dplyr   0.8.5
## ✓ tidyr   1.0.2     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0

## ── Conflicts ───────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

library(sfsmisc)

## 
## Attaching package: 'sfsmisc'

## The following object is masked from 'package:dplyr':
## 
##     last

library(kableExtra)

## 
## Attaching package: 'kableExtra'

## The following object is masked from 'package:dplyr':
## 
##     group_rows

library(cowplot)

## 
## ********************************************************

## Note: As of version 1.0.0, cowplot does not change the

##   default ggplot2 theme anymore. To recover the previous

##   behavior, execute:
##   theme_set(theme_cowplot())

## ********************************************************

library(R0)

## Loading required package: MASS

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

library(plyr)

## ------------------------------------------------------------------------------

## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)

## ------------------------------------------------------------------------------

## 
## Attaching package: 'plyr'

## The following objects are masked from 'package:dplyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize

## The following object is masked from 'package:purrr':
## 
##     compact

library(mgcv)

## Loading required package: nlme

## 
## Attaching package: 'nlme'

## The following object is masked from 'package:dplyr':
## 
##     collapse

## This is mgcv 1.8-31. For overview type 'help("mgcv-package")'.

library(lhs)
library(minpack.lm)

#Files
params.file   <- "parameters.csv"
resource.file <- "resources.csv"

#Functions
source("modelo.edad.R")
source("read.parameters.R")
source("auxiliary.functions.R")
source("estimation.functions.R")
source("model.fit.R")

Datos

Los parámetros del modelo por edad:

	Category.1	Category.2	Category.3	Category.4
Category interpretation	< 20	20 to 60	60 to 80	80+
R0	2.23	2.23	2.23	2.23
Days on mild case before becoming severe	6	6	6	6
Days asymptomatic contagious before recovery	8	8	8	8
Days of mild case before recovery	8	8	8	8
Days of severe case before recovery with hospitalization	15	15	15	15
Days of severe case before recovery without hospitalization	15	15	15	15
Days of severe case before critical with hospitalization	1	1	1	1
Days of severe case before critical without hospitalization	1	1	1	1
Days of critical case before recovery	12	12	12	12
Days of critical case before death	7.5	7.5	7.5	7.5
Days of latency (with virus but not contagious)	4	4	4	4
Contact rate for mild cases	Estimate from R0	Estimate from R0	Estimate from R0	Estimate from R0
Contact rate for severe cases	0	0	0	0
Contact rate for critical cases	0	0	0	0
Contact rate for asymptomatic cases	Estimate from R0	Estimate from R0	Estimate from R0	Estimate from R0
Probability of severe given mild	0.16	0.16	0.16	0.16
Proportion of asymptomatic cases	0.37	0.37	0.37	0.37
Probability of death given critical with hospitalization	0.11	0.1771	0.435	0.86
Probability of death given critical without hospitalization	1	1	1	1
Probability of critical given severe with hospitalization	0.01	0.185	0.36	0.09
Probability of critical given severe without hospitalization	0.5	0.5	0.5	0.5
Initial susceptible cases	0.34672723	0.54168417	0.09583051	0.01575435
Initial mild cases	0.000000189	0.00000294	0.000000574	0.000000031
Initial latent cases	0	0	0	0
Initial severe cases	0	0	0	0
Initial critical cases	0	0	0	0
Initial asymptomatic cases	0	0	0	0

Los recursos que consideramos son:

	V2
Default limit of severe cases	2e+03
Default limit of critical cases	1e+02
Extended limit of severe cases	1e+03
Extended limit of critical cases	5e+00
Total population	1e+06

Modelo

Lectura de parámetros:

#Parameter list
initial.list <- read.parameters(params.file, resource.file = resource.file)
params       <- initial.list$params
state        <- initial.list$state

Ejecución simple

Para correr el modelo basta con ejecutar el siguiente código

#You can run model with run model.instance
ode.model <- run.model.continuous(params, state,  init.time = 0, end.time = 150)

Y graficar los casos infectados totales

ggplot.epidemiological.lines.infected(ode.model)

o bien seleccionando el caso:

ggplot.epidemiological.lines.infected(ode.model, infect.cats = c("I2","I3"))

El modelo es un modelo por categorías por tanto también se puede generar una gráfica de los estados infectados por categoría:

ggplot.epidemiological.lines.infected.cat(ode.model)

El tiempo para el pico se puede calcular mediante:

time.of.peak(ode.model)

## $Mild
##         Peak of I1 Time to peak of I1 
##          0.1696552         71.7000000 
## 
## $Severe
##         Peak of I2 Time to peak of I2 
##         0.02507798        78.00000000 
## 
## $Critical
##         Peak of I3 Time to peak of I3 
##         0.01921586        83.80000000

Variaciones de política

El modelo permite implementar distintas opciones de política en distintos sectores. En el siguiente ejemplo evaluamos:

La evolución usual de la enfermedad
1. Una política de cuarentena poblacional
2. Una política de cuarentena periódica “hoy no circula”
3. Una política de crecimiento usual de la enfermedad bajo un cambio en el \(R_0\)
4. Una política de cuarentena poblacional bajo un cambio en el \(R_0\)
5. Una política de cuarentena a adultos mayores de 60 años

#You can run model with run model.instance
model.1 <- run.model.continuous(params, state,  init.time = 0, end.time = 40)

#And then add quarantine
state   <- quarantine_all(model.1$state)
model.2 <- run.model.continuous(params, state, init.time = 40, end.time = 70)

#Then continue model wth periodic quarantine
state   <- unquarantine_all(model.2$state)
model.3 <- run.model.periodic(params, state, init.time = 70, end.time = 120,
                                  periodicity = 7, days = 2)

#Change R0
params  <- rnought.change.gamma.1(1.2, params)
model.4 <- run.model.continuous(params, model.3$state, init.time = 120, end.time = 150)

#And then add quarantine
state   <- quarantine_all(model.4$state)
params  <- rnought.change.gamma.1(2.1, params)
model.5 <- run.model.continuous(params, state, init.time = 150, end.time = 175)

#And then add quarantine for specific subgroups
state   <- quarantine_k(model.5$state, k = c(3,4))
model.6 <- run.model.continuous(params, state, init.time = 175, end.time = 200)

#Add quarantine only to infected
params$q2 <- rep(1/2, 4)
model.7   <- run.model.continuous(params, model.6$state, init.time = 200, 
                                  end.time = 250)

ggplot() +
  geom_line(aes(x = time, y = I1 + QI, 
                color = "Evolución usual"), data = model.1$dats) + 
  geom_line(aes(x = time, y = I1 + QI,
                color = "Cuarentena poblacional"), 
            data = model.2$dats) + 
  geom_line(aes(x = time, y = I1 + QI,
                color = "Cuarentena periódica"), data = model.3$dats) + 
  geom_line(aes(x = time, y = I1 + QI,
                color = "Evolución usual con R0 = 1.2"), data = model.4$dats) +
  geom_line(aes(x = time, y = I1 + QI,
                color = "Cuarentena poblacional  con R0 = 2.1"), 
            data = model.5$dats) +
  geom_line(aes(x = time, y = I1 + QI,
                color = "Cuarentena de mayores de 60 con R0 = 2.1"), 
            data = model.6$dats) +
  geom_line(aes(x = time, y = I1 + QI,
                color = "Sólo cuarentena del 50% de infectados"), 
                data = model.7$dats) +
  theme_classic() +
  ggtitle("Distintos esquemas de cuarentena del modelo") +
  ylab("Infectados (proporción)") +
  xlab("Tiempo")

Casos acumulados

Los casos acumulados pueden calcularse con la función cummulative.cases como se explica a continuación:

#Parameter list
initial.list    <- read.parameters(params.file)
ode.model       <- run.model.continuous(initial.list$params, initial.list$state,  
                                        init.time = 0, end.time = 100)
acumulados      <- cummulative.cases(ode.model)
ggplot.cummulative.lines.all(acumulados)

Desabasto hospitalario y efecto cascada

El modelo actual incluye un efecto de saturación del sistema “efecto cascada”. En los parámetros se puede establecer los niveles de saturación del sistema a partir de la base de recursos:

#Parameter list
initial.list <- read.parameters(params.file, estimate.saturation = TRUE, 
                                resource.file = resource.file)
params       <- initial.list$params
state        <- initial.list$state

También puede cambiarse de manera manual en los parámetros:

params$saturationI2 <- 0.001
params$saturationI3 <- 0.0001
saturated.model <- run.model.continuous(params, state,  init.time = 0, end.time = 100)

La función de tiempo para saturación estima la cantidad de días para saturación del sistema para hospitalizados y cuidado intensivo:

time.to.saturation(saturated.model)

##   I2   I3 
## 49.1 41.7

Si se desea correr el modelo sin saturación basta poner los parámetros de saturación en infinito:

params$saturationI2 <- Inf
params$saturationI3 <- Inf
unsaturated.model <- run.model.continuous(params, state,  init.time = 0, end.time = 100)

Y pueden compararse ambos escenarios:

ggplot() +
  geom_line(aes(x = time, y = M, color = "Saturated"), data = saturated.model$dats) +
  geom_line(aes(x = time, y = M, color = "Unaturated"), data = unsaturated.model$dats) +
  theme_bw() +
  ggtitle("Mortalidad por COVID-19")

Gráficas

Dado un modelo podemos producir las siguientes gráficas:

#Parameter list
initial.list <- read.parameters(params.file, resource.file = resource.file)
params       <- initial.list$params
state        <- initial.list$state
ode.model    <- run.model.continuous(params, state,  init.time = 0, 
                                     end.time = 100)

Gráficas de infectados ya sea agregados o por categoría:

ggplot.epidemiological.lines.infected(ode.model)

ggplot.epidemiological.lines.infected.cat(ode.model)

Gráficas de todas las curvas epidemiológicas agregadas o por edad:s

ggplot.epidemiological.lines.all(ode.model)

ggplot.epidemiological.lines.all.cat(ode.model)

Combinación de modelos

Para utilizar las funciones anteriores varios modelos pueden combinarse con la función model.bind como en el siguiente ejemplo:

#Parameter list
initial.list <- read.parameters(params.file)
params       <- initial.list$params
state        <- initial.list$state

#You can run model with run model.instance
model.1 <- run.model.continuous(params, state,  init.time = 0, end.time = 20)
model.2 <- run.model.continuous(params, model.1$state,  init.time = 20, end.time = 30)
model.3 <- run.model.continuous(params, model.2$state,  init.time = 30, end.time = 70)
model.4 <- run.model.continuous(params, model.3$state,  init.time = 70, end.time = 200)

model.list <- list("Escenario 1" = model.1, "Escenario 2" = model.2,
                   "Escenario 3" = model.3,
                   "Escenario 4" = model.4)

#Juntar todos los modelos
all.models <- model.bind(model.list)

Y pueden aplicarse las funciones de graficación previas al modelo:

ggplot.epidemiological.lines.infected(all.models)

O bien las de acumular casos:

acumulados <- cummulative.bind(model.list)
ggplot.cummulative.lines.all(acumulados)

También se puede usar la función geom_line.variable.from.list para poner líneas en objetos ggplot2:

plot <- geom_line.variable.from.list(model.list, varname = "I1") 
plot <- geom_line.variable.from.list(model.list, plot, "I2", linetype = "dashed")
plot <- geom_line.variable.from.list(model.list, plot, "I3", linetype = "dotted")
plot + theme_bw()

Estimación estadística

El \(R_0\) puede estimarse a partir de los casos acumulados y la fecha mediante:

#Casos acumulados de México:
total.cases <- c(2785, 2439, 2143, 1890, 1688, 1510, 1378, 1215, 1094, 
                       993, 848, 717, 585, 475, 405, 367, 316, 251, 203, 164, 
                       118, 93, 82, 53, 41, 26, 15, 11, 7)
date.cases        <- c("7/4/2020", "6/4/2020", "5/4/2020", "4/4/2020", "3/4/2020",
                       "2/4/2020", "1/4/2020", "31/3/2020", "30/3/2020", "29/3/2020",
                       "28/3/2020", "27/3/2020", "26/3/2020", "25/3/2020", "24/3/2020",
                       "23/3/2020", "22/3/2020", "21/3/2020", "20/3/2020", "19/3/2020",
                       "18/3/2020", "17/3/2020", "16/3/2020", "15/3/2020", "14/3/2020",
                       "13/3/2020", "12/3/2020", "11/3/2020", "10/3/2020")
date.cases        <- as.Date(date.cases, format = "%d/%m/%Y")

#Estimación del R0
R0 <- estimate.R0(total.cases, date.cases)

## Warning in est.R0.TD(epid = c(2785, 2439, 2143, 1890, 1688, 1510, 1378, :
## Simulations may take several minutes.

## Warning in est.R0.TD(epid = c(2785, 2439, 2143, 1890, 1688, 1510, 1378, : Using
## initial incidence as initial number of cases.

El cual puede verse en la siguiente tabla:

kable(tail(R0)) %>%  kable_styling(bootstrap_options = c("striped", "hover"))

	Fecha	R0	Lower bound R0	Upper bound R0
23	2020-04-01	1.174329	1.133294	1.216945
24	2020-04-02	1.185190	1.144658	1.225324
25	2020-04-03	1.193777	1.156106	1.231530
26	2020-04-04	1.207371	1.171411	1.243789
27	2020-04-05	1.215455	1.181736	1.249187
28	2020-04-06	1.220810	1.189417	1.252325

O bien puede graficarse el R0:

plot.R0(R0)

La estimación de casos diarios puede hacerse mediante un modelo Poisson:

dats <- data.frame(
  Fecha = rev(date.cases),
  n = c(7,4,4,11,15,12, 29, 11, 25, 26, 39, 48, 65, 51, 38, 70, 110, 132,
        131, 145, 101, 121, 163, 132, 178, 202, 253, 296, 346))

predicted.cases <- estimate.daily.cases(dats, dias.predict = 10, 
                                        confidence = 0.99, nsim = 1000)

kable(tail(predicted.cases)) %>%  kable_styling(bootstrap_options = c("striped", "hover"))

	Lower	Median	Upper	Fecha	n
34	505.960	612	753.075	2020-04-12	NA
35	548.970	698	845.115	2020-04-13	NA
36	618.985	785	976.160	2020-04-14	NA
37	703.985	889	1101.095	2020-04-15	NA
38	792.980	1008	1264.275	2020-04-16	NA
39	866.875	1142	1449.090	2020-04-17	NA

La estimación de casos acumulados (no recomendada) se puede hacer mediante la función estimate.cummulative.cases:

dats              <- dats %>% mutate(Cummulative = cumsum(n))
pred.cummulative.cases <- estimate.cummulative.cases(dats, nsim = 100)

kable(tail(pred.cummulative.cases)) %>%  kable_styling(bootstrap_options = c("striped", "hover"))

	Fecha	Intervalo bajo	Casos predichos	Intervalo superior	Observados
34	2020-04-12	4032.629	4344.035	4661.579	NA
35	2020-04-13	4326.328	4627.236	5010.264	NA
36	2020-04-14	4648.513	4929.536	5257.729	NA
37	2020-04-15	4902.355	5309.452	5556.290	NA
38	2020-04-16	5250.366	5591.110	5901.059	NA
39	2020-04-17	5600.416	5914.912	6216.977	NA

Existen funciones de graficación tanto para los casos diarios como los acumulados:

plot.cummulative.cases(pred.cummulative.cases)

Podemos graficar de igual forma los casos puntuales estimados:

plot.predicted.cases(predicted.cases)

El modelo puede ser ajustado a partir de los datos iniciales:

initial.list <- model.from.cummulative.cases(total.cases, 
                                             date.cases, params.file)

## Warning in est.R0.TD(epid = c(2785, 2439, 2143, 1890, 1688, 1510, 1378, :
## Simulations may take several minutes.

## Warning in est.R0.TD(epid = c(2785, 2439, 2143, 1890, 1688, 1510, 1378, : Using
## initial incidence as initial number of cases.

Modelo que se corre extactamente igual:

model.A <- run.model.continuous(initial.list$params, 
                                initial.list$state,  init.time = 0, end.time = 65)

Recordemos que se puede graficar y los objetos son cosas de ggplot:

numero.personas <- 120000000
ggplot.epidemiological.lines.infected(model.A, infect.cats = c("I2","I3"),
                                      scale = numero.personas,
                                      date.init = date.cases[length(date.cases)]) +
  geom_hline(aes(yintercept = 1000), linetype = "dotted") +
  annotate("text", label = "UCI Disponibles", x = date.cases[length(date.cases)],
           y = 1200, hjust = 0) + 
  ggtitle("Casos estimados por tipo de infección") +
  theme_bw()

Ajuste de modelo a datos

Consideremos casos observados (y escalados) a casos por cada 10,000 habitantes. El escalamiento puede que sea necesario para mejorar la optimización.

#Casos por cada 10,000 habitantes
observed.cases <- 10000*c(3181, 2785, 2439, 2143, 1890, 1688, 1510, 1378, 1215, 1094, 
                       993, 848, 717, 585, 475, 405, 367, 316, 251, 203, 164, 
                       118, 93, 82, 53, 41, 26, 15, 11, 7)/120000000
date.cases        <- c("8/4/2020", "7/4/2020", "6/4/2020", "5/4/2020", "4/4/2020", "3/4/2020",
                       "2/4/2020", "1/4/2020", "31/3/2020", "30/3/2020", "29/3/2020",
                       "28/3/2020", "27/3/2020", "26/3/2020", "25/3/2020", "24/3/2020",
                       "23/3/2020", "22/3/2020", "21/3/2020", "20/3/2020", "19/3/2020",
                       "18/3/2020", "17/3/2020", "16/3/2020", "15/3/2020", "14/3/2020",
                       "13/3/2020", "12/3/2020", "11/3/2020", "10/3/2020")
date.cases        <- as.Date(date.cases, format = "%d/%m/%Y")

#Creación de la base de datos necesaria para el escalamiento
dats              <- data.frame(I1 = rev(observed.cases), Fecha = rev(date.cases))
dats$time         <- dats$Fecha - dats$Fecha[1]
dats              <- dats %>% dplyr::select(-Fecha)

#Adjust model to parameters
initial.list <- read.parameters(params.file)
state        <- initial.list$state
params       <- initial.list$params

#Fit value to model
fitval       <- fit.model(dats, state)
model.info.1 <- run.fitted.model(fitval, state, params, end.time = 30)
model.info.2 <- run.fitted.model(fitval, state, params, end.time = 250)

#Graficación del ajuste
ggplot.data.to.fitted(model.info.1$cummulative, dats)

#Modelo ajustado completo:
ggplot.epidemiological.lines.infected(model.info.2$incidence)

Ejemplo

En este ejemplo se evalúan distintos escenarios a partir del código. En primer lugar se llaman los paquetes y funciones a utilizar:

#Required libraries to run the model
library(deSolve)
library(tidyverse)
library(sfsmisc)
library(kableExtra)
library(cowplot)
library(R0)
library(plyr)
library(mgcv)

#Files
params.file   <- "parameters.csv"
resource.file <- "resources.csv"

#Functions
source("modelo.edad.R")
source("read.parameters.R")
source("auxiliary.functions.R")
source("estimation.functions.R")
source("model.fit.R")

Luego se construye el modelo a partir de casos observados:

#Estimate cummulative cases for Mexico
total.cases <- c(2785, 2439, 2143, 1890, 1688, 1510, 1378, 1215, 1094, 
                       993, 848, 717, 585, 475, 405, 367, 316, 251, 203, 164, 
                       118, 93, 82, 53, 41, 26, 15, 11, 7)
date.cases        <- c("7/4/2020", "6/4/2020", "5/4/2020", "4/4/2020", "3/4/2020",
                       "2/4/2020", "1/4/2020", "31/3/2020", "30/3/2020", "29/3/2020",
                       "28/3/2020", "27/3/2020", "26/3/2020", "25/3/2020", "24/3/2020",
                       "23/3/2020", "22/3/2020", "21/3/2020", "20/3/2020", "19/3/2020",
                       "18/3/2020", "17/3/2020", "16/3/2020", "15/3/2020", "14/3/2020",
                       "13/3/2020", "12/3/2020", "11/3/2020", "10/3/2020")
date.cases        <- as.Date(date.cases, format = "%d/%m/%Y")

#Adjust model to parameters
initial.list <- model.from.cummulative.cases(total.cases, 
                                             date.cases, params.file)

## Warning in est.R0.TD(epid = c(2785, 2439, 2143, 1890, 1688, 1510, 1378, :
## Simulations may take several minutes.

## Warning in est.R0.TD(epid = c(2785, 2439, 2143, 1890, 1688, 1510, 1378, : Using
## initial incidence as initial number of cases.

Correremos distintos escenarios del modelo. En primer lugar, un modelo sin cuarentena:

#Run model without quarantine
model.caeteris.paribus <- run.model.continuous(initial.list$params, 
                                               initial.list$state,  
                                               init.time = 0, end.time = 1000)

Después un modelo con aislamiento de 50% de sintomáticos:

#Run model with isolation
params      <- initial.list$params
params$q2   <- rep(1, 4)/2
model.start <- run.model.continuous(initial.list$params, 
                                    initial.list$state,  
                                    init.time = 0, end.time = 20)
model.isolation <- run.model.continuous(params, 
                                        model.start$state,  
                                        init.time = 20, end.time = 1000)
model.isolation <- model.bind(list(model.start, model.isolation))

Un modelo con cuarentena periódica:

#Periodic quarantine
params$q2      <- rep(0, 4)
model.periodic <- run.model.periodic(params, 
                                     model.start$state,  
                                     periodicity = 7, days = 2, 
                                     quarantine.proportion = 1,
                                     init.time = 20, end.time = 1000)
model.periodic <- model.bind(list(model.start, model.periodic))

Un modelo con cuarentena poblacional:

#Cuarentena poblacional
state                     <- quarantine_all(model.start$state)
model.quarantine.all.next <- run.model.continuous(initial.list$params, 
                                                  state,  
                                                  init.time = 20, end.time = 110)

#Después de la cuarentena poblacional
state                     <- unquarantine_all(model.quarantine.all.next$state)
model.quarantine.all.last <- run.model.continuous(initial.list$params, 
                                                  state,  
                                                  init.time = 110, end.time = 1000)

#Juntamos los modelos
model.quarantine.all <- model.bind(list(model.start, 
                                        model.quarantine.all.next,
                                        model.quarantine.all.last))

Y un modelo con cuarentena para mayores de 60:

#Run model with quarantine for individuals > 60
state              <- quarantine_k(model.start$state, k = c(3,4))
model.quarantine.k <- run.model.continuous(initial.list$params, 
                                           state,  
                                           init.time = 20, end.time = 320)

#Run part where quarantine ends:
state                   <- unquarantine_k(model.quarantine.k$state, k = c(3,4))
model.quarantine.k.last <- run.model.continuous(initial.list$params, 
                                                  state,  
                                                  init.time = 320, end.time = 1000)

#Bind models
model.quarantine.k <- model.bind(list(model.start, 
                                      model.quarantine.k, 
                                      model.quarantine.k.last))

Finalmente, comparamos todos los modelos en una gráfica:

ggplot() +
  geom_line(aes(x = time/3, y = I2, 
                linetype = "Severe cases (hospitalized)", 
                color = "No Intervention"), 
            data = model.caeteris.paribus$dats) +
  geom_line(aes(x = time/3, y = I2, 
                linetype = "Severe cases (hospitalized)", 
                color = "Isolation of 50% symptomatic individuals"), 
            data = model.isolation$dats) +
  geom_line(aes(x = time/3, y = I2, 
                linetype = "Severe cases (hospitalized)", 
                color = "Population Quarantine (100% per 1 month)"), 
            data = model.quarantine.all$dats) +
  geom_line(aes(x = time/3, y = I2, 
                linetype = "Severe cases (hospitalized)", 
                color = "Periodic 50% population quarantine (2 days/week)"), 
            data = model.periodic$dats) +
  geom_line(aes(x = time/3, y = I2, 
                linetype = "Severe cases (hospitalized)", 
                color = "Isolation of individuals > 60 (100 days)"), 
            data = model.quarantine.k$dats) +
  geom_line(aes(x = time/3, y = I3, 
                linetype = "Critical cases (ICU)", 
                color = "No Intervention"), 
            data = model.caeteris.paribus$dats) +
  geom_line(aes(x = time/3, y = I3, 
                linetype = "Critical cases (ICU)", 
                color = "Isolation of 50% symptomatic individuals"), 
            data = model.isolation$dats) +
  geom_line(aes(x = time/3, y = I3, 
                linetype = "Critical cases (ICU)", 
                color = "Population Quarantine (100% per 1 month)"), 
            data = model.quarantine.all$dats) +
  geom_line(aes(x = time/3, y = I3, 
                linetype = "Critical cases (ICU)", 
                color = "Periodic 50% population quarantine (2 days/week)"), 
            data = model.periodic$dats) +
  geom_line(aes(x = time/3, y = I3, 
                linetype = "Critical cases (ICU)", 
                color = "Isolation of individuals > 60 (100 days)"), 
            data = model.quarantine.k$dats) +
  theme_bw() +
  ggtitle("Effects of distinct isolation and quarantine scenarios") +
  xlab("Days since the beginning of the pandemic") +
  ylab("Infected individuals (percent of total population)") +
  scale_y_continuous(labels = scales::percent) +
  scale_linetype_manual("Condition", 
                        values = c("Severe cases (hospitalized)" = "solid",
                                   "Critical cases (ICU)" = "dotted")) +
  scale_color_brewer("Intervention", palette = "Dark2")

TODO

Crear gráficas de demanda hospitalaria
Agregar fit estilo el paper de Marissa incluyendo la variable k de la ecuación 3
Agregar migración

Colaboración

Para colaborar en el modelo requieres conocimiento de epidemiología matemática, y/o epidemiología, sistemas dinámicos. Cualquier agregado está guiado por el comité científico.

Agregados de programación siempre y cuando se limiten a mejorar las gráficas, la evaluación del modelo pero no el modelo.

Traducciones bienvenidas.

Código

Ver código en Github

Referencias

Brandenburg, A. 2020. “Quadratic Growth During the 2019 Novel Coronavirus Epidemic.” arXiv:2002.03638v2.

Brauer, F. 2017. “Mathematical epidemiology: Past, present, and future.” Infectious Disease Modelling 2: 113–27. https://doi.org/{10.1016/j.idm.2017.02.001}.

Brauer, F, and Castillo-Chávez C. 2011. Mathematical Models in Population Biology and Epidemiology. Vol. 2. Springer. https://doi.org/{https://www.springer.com/gp/book/9781441931825}.

Ferguson, N, and et. al. 2020. “Impact of non-pharmaceutical interventions (NPIs) to reduce COVID-19 mortality and healthcare demand.” Imperial College Report 2. https://doi.org/{ https://doi.org/10.25561/77482}.

Nadim, SS, I Gosh, and J Chattopadhyay. 2019. “Short-Term Predictions and Prevention Strate-Gies for Covid-2019: A Model Based Study.” arXiv:2003.08150v1.

Modelo de M clases (edad/grupo socioeconómico/comorbilidad) con distintos escenarios de cuarentena y desabasto hospitalario para COVID-19

Equipo de modelaje del coronavirus

Last Update 2020-04-09 21:07:30

Licencias

Modelo

Variables

Ecuaciones

Cuarentenas

Condiciones iniciales

Programación

Datos

Modelo

Ejecución simple

Variaciones de política

Casos acumulados

Desabasto hospitalario y efecto cascada

Gráficas

Combinación de modelos

Estimación estadística

Ajuste de modelo a datos

Ejemplo

TODO

Colaboración

Código

Referencias