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.
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:
\[\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}\]
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\) |
Para mitigar el contagio se implementaron 3 tipos de cuarentena en el modelo:
\[ 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}. \]
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, \]
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")
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 |
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 |
Lectura de parámetros:
#Parameter list
initial.list <- read.parameters(params.file, resource.file = resource.file)
params <- initial.list$params
state <- initial.list$state
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)