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

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)