Modelado de la epidemia del coronavirus

Wilfried Coenen (http://fluidosuc3m.es/people/wcoenen/)
Ana Medina Palomo

Primera publicación: 18 de marzo de 2020
Última actualización: 2 de abril de 2020


English version (Sorry, under construction)

 

Disclaimer 1: No soy epidemiólogo. Este proyecto surge de un par de tardes aburridas en confinamiento con mi mujer en Leganés, y sus resultados deberán ser interpretados como tal. Por favor, no vayáis reservando vuestras vacaciones en base a las predicciones de aquí abajo, porque seguramente os vais a arruinar.

Disclaimer 2: Mientras estaba escribiendo esto, he visto que hay un grupo en la Universidad Politécnica de Valencia que también se llevan aburriendo varios días. Ellos implementan el SIR estándar, pero parecen acoplarlo con herramientas de machine learning, me imagino para ir ajustando la tasa de transmisión de manera automática según van saliendo datos. Sus resultados están aquí, con informes diarios muy bien presentados.

Disclaimer 3: Repito que es la primera vez que miro modelos de propagación de epidemias. Estaré encantado de recibir feedback de gente que realmente sabe (correo / twitter @_w_coenen_).

 

1. Introducción

Este trabajo es un intento de modelado — muy muy simplista — de la epidemia de coronavirus 2020, con comparaciones con los datos publicados para España e Italia. Se basa en el modelo compartimental SEIR, pero con una ligera modificación en forma de un tiempo de retraso para tomar en cuenta la duración de la enfermedad. La explicación detallada del modelo se puede leer a continuación. Después se discute como tomar en cuenta la diferencia entre el número de infectados detectados y el número real de infectados, la mortalidad del virus, y como modelar los efectos de la cuarentena. Finalmente, se presentan resultados para España e Italia.

 

2. El modelo

2.1 Ecuaciones de evolución

En el modelo SEIRD, una población (N) se divide en cinco grupos: los susceptibles (S), los expuestos (E), los infectados (I), los recuperados (R), y los fallecidos (D). Cada persona sólo puede estar en un grupo a la vez. Los infectados pueden pasar el virus a los miembros del grupo de susceptibles, convirtiéndose en expuestos (E). Después de un tiempo de incubación \tau_i, los expuestos pasan al grupo de infectados. Los infectados o bien se recuperan después de un tiempo \tau_r, o bien fallecen después de un tiempo \tau_d. La fracción que se muere es \mu, y por tanto la fracción que se recupera es 1-\mu. Los recuperados no pueden volver a contraer la enfermedad. El ritmo al cual se generan expuestos es proporcional al número de infectados elevado a una potencia p entre 0 y 1, I^p, al porcentaje de susceptibles que hay en la población, S/N, y al parámetro de transmisión \beta.

El sistema de ecuaciones diferenciales que gobierna la evolución temporal del número de miembros en cada grupo es

 \begin{align*} \frac{dS}{dt} & = - \frac{\beta S(t) I(t)}{N}, \\ \frac{dE}{dt} & = \frac{\beta S(t) I(t)}{N} - E(t - \tau_i), \\ \frac{dI}{dt} & = E(t - \tau_i) - \mu E(t - \tau_d) - (1-\mu) E(t - \tau_r), \\ \frac{dD}{dt} & = \mu E(t - \tau_d), \\ \frac{dR}{dt} & = (1-\mu) E(t - \tau_r) \end{align*}

con condiciones iniciales adecuadas.

La diferencia del modelo usado en este trabajo con el modelo SIR estándar [3] es el uso de los tiempos de retraso \tau_i, \tau_r y \tau_d, por ejemplo cuando modelamos la pérdida de infectados por fallecimiento E(t-\tau_d). En el modelo SIR estándar el ritmo al cual perdemos infectados es proporcional al número instantáneo de infectados, -\gamma I — no me parece muy correcto. La introducción del tiempo de retraso convierte las ecuaciones de evolución en ecuaciones diferenciales con retraso (delay differential equations, dde). Los resultados numéricos del modelo se han obtenidos con la rutina dde23 de MATLAB.

En los cálculos, se ha tomado el valor \tau_i = 1 para el tiempo de incubación. El tiempo de enfermedad se ha inferido de los datos; para España \tau_r = 10 y \tau_d = 16. El ratio de mortalidad (ver sección 5) es \mu = 0.007.

2.2 Parámetros de transmisión y de crecimiento

En el modelo de arriba, el parámetro de transmisión, \beta, es proporcional al número de encuentros que hay entre miembros de la población y a la probabilidad de contagio que hay en cada encuentro. El distanciamiento social reduce el número de contactos y por tanto reduce el valor de \beta.

La potencia p modela el grado de mezclado en la población. El caso p = 1 corresponde con un mezclado homogéneo, y resultará en la fase inicial (cuando S/N \simeq 1) en un crecimiento exponencial del número de infectados. En la realidad, los procesos de mezclado de la población raramente son homogéneos, y por tanto la mayoría de las epidemias muestran un crecimiento sub-exponencial [3, 4]. La primera figura de abajo muestra un árbol de contagio epidémico con 12 generaciones de transmisión de un caso con crecimiento exponencial (izquierda) y de un caso con crecimiento sub-exponencial (derecha). La segunda figura ilustra la influencia del parámetro p sobre las curvas de crecimiento de una epidemia. Es de esperar que las medidas de distanciamiento social también causan una reducción del valor de p.

Al contrario de los países Europeos, donde el crecimiento del coronavirus hasta ahora es casi exponencial, los países Asiáticos se caracterizan por un crecimiento muy sub-exponencial [2]. En los resultados para España, se ha tomado el valor p = 0.85, un valor parecido al valor estimado para la pandemia de la gripe Española en 1918 [4 (fig. 6)].


Fuente: [3].


Fuente: [3].

 

3. El número real de infectados y el tiempo de retraso en la publicación de los casos

Debido a la limitada capacidad de hacer pruebas de COVID-19, el número real de infectados es mucho mayor que el número de infectados detectados [5-8]. Sin embargo, el cociente exacto entre ambos números no es conocido en este momento. Algunos autores [8] estiman hasta un factor 95 de diferencia para España, y un factor 63 para Italia. Otros [5] estiman un factor 20 para España y Italia. Este último es el valor que se ha empleado en el presente trabajo. Ademas del factor de proporcionalidad, hay un desfase temporal entre la curva real del número de infectado y la curva de los infectados detectados. Teniendo en cuenta la siguiente cadena de sucesos, podemos estimar que si una persona se expone al virus el día 0, el día 1 se vuelve infeccioso, alrededor del día 5 tiene síntomas, el día 8 se le hace el test, el día 9 salen los resultados del test positivo, y el día 11 aparece como parte de los números publicados [9]. Por otro lado, podemos asumir que el número de los fallecidos publicados sí es igual al número de fallecidos reales instantáneos—aunque aquí también hay polémica [10]. De las misma manera, se asume que el número publicado de recuperados es instantáneo, pero que se relaciona con el número real de casos por un factor de proporcionalidad de 20. La siguiente figura ilustra la relación entre las curvas de casos reales y detectados.

De manera matemática, si C_p, R_p y D_p indican los números publicados de infectados acumulados, recuperados y fallecidos, los números reales de infectados acumulados, recuperados y fallecidos serán

 \begin{align*} C_r(t) &= r_p C_p(t-\tau_p), \\ R_r(t) &= r_p R_p(t), \\ D_r(t) &= D_p(t), \end{align*}

donde r_p es el factor de proporcionalidad (20 para España e Italia), y \tau_p es el tiempo entre la infección de una persona y la publicación de su caso (11 días para España e Italia).

 

4. Casos activos

Hay cierta confusión sobre los casos activos. Visto el retraso considerable en la publicación de las infecciones, no se puede simplemente restar los recuperados publicados y muertes publicados de los casos acumulados publicados para obtener los casos activos conforme con los datos publicados. Es mejor trabajar con los casos reales, dónde sí podemos hacer esa resta, i.e. el número real de casos activos, I_r, es

 I_r(t) = C_r(t) - R_r(t) - D_r(t).

A posteriori, podemos ir al sistema de referencia de los casos publicados como

 I_p(t) = I_r(t+\tau_p) / r_p

para estimar el número de casos activos conforme con los datos publicados. La siguiente figura para Hubei, China, muestra la relación entre todos las variables reales y publicados (usando, para ser consistente con lo anterior, \tau_p = 11 y r_p = 20, los valores adecuados para Hubei pueden ser diferentes).

 

5. Mortalidad

La mortalidad real del virus se caracteriza por el número de fallecidos dividido por el número real de casos (infection fatality rate (IFR). Sin embargo, en los medios normalmente hablan de la mortalidad como el número de fallecidos dividido por el número de casos detectados (case fatality rate (CFR). La relación entre el IFR y el CFR es igual el factor de detección, descrito en el anterior párrafo (aquí 20 para España e Italia). Como se puede apreciar en el ejemplo de la siguiente figura, ambos coeficientes pueden cambiar en el tiempo, tal que solo tiene sentido medirlos en el estado pseudo-estacionario que se alcanza al final de la epidemia, o bien en la fase inicial si esa es puramente exponencial. En el modelo matemático, la mortalidad real entra a través del parámetro \mu. Para los resultados numéricos de España, se ha usado \mu = 0.007. Esto significa que de cada 1000 personas que se infectan de coronavirus, 7 fallecen.

 

6. Como modelar los efectos de la cuarentena

Las medidas de distanciamiento social reducen el número de contactos por unidad de tiempo entre los miembros de la población. Por tanto, en el modelo actúan a través de una disminución de la tasa de transmisión \beta y posiblemente también del nivel de crecimiento subexponential p. Sin embargo, las medidas de cuarentena estricta, como la que hay en España, dividen la población en dos de una manera fundamental: los que hacen la cuarentena, y los que, por su profesión, no lo hacen. Esta división está en desacuerdo con la hipótesis básica de los modelos tipo SIR, i.e. que hay un mezclado más o menos homogéneo de la población. En la cuarentena, los dos grupos tienen una dinámica muy diferente: el grupo que está en casa únicamente sale de vez en cuando a comprar al Mercadona, donde tiene una muy pequeña probabilidad de infectarse, mientras el grupo que sigue yendo a trabajar se sigue mezclando con mucha agilidad, por ejemplo en el transporte público o en el trabajo. Por eso tiene sentido modelar la evolución de los infectados en todo la población del país como la suma de las evoluciones de dichos dos grupos, cada una con su propia dinámica. En realidad, como dos semanas después del inicio de la cuarentena se enunciaron medidas adicionales, conviene separar la población en tres grupos: A, B y C.

Antes del inicio de la cuarentena, los tres grupos tienen la misma tasa de transmisión. Cuando empieza la cuarentena, el grupo A, que se tiene que quedar en casa, baja su tasa de transmisión de manera considerable, por ejemplo con un 90%. Los grupos B y C, que siguen yendo a trabajar, un asustados por la situación, bajan también su tasa de transmisión, pero en menor medida, por ejemplo con un 50%. Con las medidas adicionales dos semanas después, el grupo B se tiene que quedar en casa, mientras el grupo C sigue yendo a trabajar. De esa manera, la tasa de transmisión del grupo B disminuye considerablemente hasta alcanzar el nivel del grupo A, mientras la tasa de transmisión del grupo C sigue igual. Es difícil estimar el tamaño relativo de los tres grupos. Aquí se asume que el grupo A, i.e. la parte grande de la población que hace la cuarentena, es el 80% de la población, mientras que los grupos B y C constituyen cada uno el 10% de la población. Si se resuelve numéricamente el modelo de evolución en cada grupo y se suman los resultados parciales, obtenemos la siguiente gráfica. Se puede observar que los efectos de la cuarentena se ven en el grupo A unos 10-15 días después su inicio, un tiempo que es del orden del tiempo \tau_p entre infección y publicación. De la misma manera, los efectos de las medidas adicionales se empiezan a ver en el grupo B un tiempo similar después de su inicio. Finalmente el grupo C llega a su pico un mes más tarde. Obviamente los resultados son muy susceptibles a los porcentajes de bajada de las tasas de transmisión, y a los tamaños relativos de los tres grupos.

 

7. Resultados

España

Ahora se presenta los resultados del modelo para España, y se compara con los datos publicados. Nótese que, poder poder hacer dicha comparación, presentamos las curvas  C_p, I_p, R_p, D_p . Hay que tener en cuenta que el número real de casos reales es mucho más alto. La relación entre los dos se explica en las secciones 3 y 4.

Según las predicciones, alrededor del 7 de abril los casos activas llegarán a su máximo, y empezarán a bajar. La asíntota de los casos acumulados es del orden de 300000, y la asíntota del número total de fallecidos del orden de 40000. Sobre el 15 de junio bajaremos otra vez de los 10000 casos activos. Hay que destacar que las predicciones a largo plazo son muy difíciles de hacer, y habrá que revisar los parámetros cada semana para ajustar la predicción.

Italia

Los resultados para Italia muestran que está en su pico de infectados activos. Los números para Italia son muy parecidos a los de España, con unos 300000 infectados acumulados finales y unos 40000 fallecidos finales.

Efectos de la cuarentena

Los efectos de la cuarentena son claros. Para los non-believers, a continuación se muestran los resultados numéricos del modelo, sin tomar en cuenta la bajada de la tasa de transmisión, i.e. como si no hubiese habido cuarentena. En ese caso, se habría alcanzado unos 300000 muertos, un orden de magnitud mayor que con la cuarentena. Conclusión, la cuarentena funciona y hay que quedarse en casa!

 

Referencias

[1] Link a los datos

[2] http://nrg.cs.ucl.ac.uk/mjh/covid19/

[3] Chowell, G., Sattenspiel, L., Bansal, S., & Viboud, C. (2016). Mathematical models to characterize early epidemic growth: A review. Physics of life reviews, 18, 66-97.

[4] Chowell, G., Viboud, C., Simonsen, L., & Moghadas, S. M. (2016). Characterizing the reproduction number of epidemics with early subexponential growth dynamics. Journal of The Royal Society Interface, 13(123), 20160659.

[5] Jombart, T., van Zandvoort, K., Russell, T., Jarvis, C., Gimma, A., Abbott, S., ... & Pearson, C. (2020). Inferring the number of COVID-19 cases from recently reported deaths. medRxiv.

[6] Russell, T. W., Hellewell, J., Jarvis, C. I., van-Zandvoort, K., Abbott, S., Ratnayake, R., ... & CMMID nCov working group. (2020). Estimating the infection and case fatality ratio for COVID-19 using age-adjusted data from the outbreak on the Diamond Princess cruise ship. medRxiv

[7] https://institucional.us.es/blogimus/2020/03/como-estimar-el-numero-de-infectados-reales-por-covid-19-el-caso-de-andalucia-e-italia/

[8] Seth Flaxman et al. (2020) Estimating the number of infections and the impact of non- pharmaceutical interventions on COVID-19 in 11 European countries

[9] https://ourworldindata.org/coronavirus

[10] https://elpais.com/sociedad/2020-03-29/cada-pais-cuenta-los-muertos-a-su-manera-y-ninguno-lo-hace-bien.html