La “gripe española” de 1918 (virus H1N1) mató a más de 300000 personas sólo en España y más de 40 millones en todo el mundo. Con la actual pandemia de coronavirus de 2020, aquella pandemia es un referente inevitable. ¿qué podemos aprender de ella?

La famosa “gripe española” de 1918 (que no se originó en España, aunque hizo estragos en el país, sino, posiblemente, en EEUU) causó estragos en todo el mundo. De ella podemos aprender muchas cosas y, en esta entrada, que voy a tratar de mantener lo más corta y sencilla posible, voy a mostrar a los lectores una de esas enseñanzas: cómo el número de infectados o de víctimas de una epidemia sigue un modelo matemático establecido. La intención en ésta serie de entradas es observar un aspecto muy concreto, en este caso, el modelo matemático que vamos a seguir en el análisis de COVID-19. No voy a entrar en detalles del fundamento matemático, ya que la intención es puramente divulgativa, de modo que el lector general pueda aprender a ver e interpretar éste tipo de análisis.

Los datos

Los datos que vamos a usar para este sencillo estudio me los ha proporcionado Antonio Quesada Ramos, quien ha publicado éste estudio sobre la epidemia de gripe de 1918 en Alcalá la Real (Jaén). Para ello, Antonio Quesada recopiló cuidadosamente los datos de fallecidos por gripe en dicha localidad. Ha tenido la amabilidad de enviarme sus datos para poder realizar algunos cálculos con ellos.

El modelo

Los fenómenos biológicos relacionados con el crecimiento, ya sea el crecimiento de un organismo (que en definitiva es el crecimiento de una población celular) o el crecimiento poblacional, tal como el número de individuos en una población o ecosistema (por ejemplo, el número de árboles en un bosque de área definida, el número de bacterias en un cultivo, el número de personas en una ciudad, el número de muertos en una guerra o el número de organismos infectados en una epidemia) se rigen por un conjunto de ecuaciones llamadas logísticas o sigmoidales, por la forma de S que presentan. Estas ecuaciones son soluciones particulares a la ecuación diferencial de Bernoulli, que describe el rate, como por ejemplo el número de fallecidos por unidad de tiempo debidos a una epidemia (y que tendrá una forma característica acampanada, como veremos). Una de estas soluciones para describir el comportamiento de mi población es la ecuación logística generalizada de cinco parámetros de Richards. Esta tiene la ventaja de que muestra asimetría y es flexible en el cálculo del punto de inflexión, por lo que puede describir bien un proceso epidémico; otras soluciones de crecimiento poblacional, como la ecuación logística de Gompertz, funcionan mejor con crecimientos como el número de árboles en un bosque.

La ecuación de Richards es:

en la que los cinco parámetros son

  • Lo : asíntota inicial (región próxima a cero pacientes)
  • L∞ : asíntota superior o “plateau”, que nos indicará que la fase epidémica ha frenado.
  • T: parámetro de forma, que define la asimetría y forma de la ecuación, dependiendo de dónde se ubique el punto de inflexión. 1/T es el factor de asimetría. Si 1/T=1, la curva sigmoidal es simétrica.
  • k: coeficiente de crecimiento, que depende de las características de la infección. Tendrá un valor entre 0 y 1 y determina a qué velocidad crecen los casos en la epidemia (en nuestro caso)
  • tm: tiempo al que se alcanzó el máximo crecimiento.

Una vez obtenidos por regresión no lineal los parámetros, se puede establecer una predicción de hacia dónde va a dirigirse la infección. Esta es una aproximación semiempírica: dados los parámetros obtenidos al acumular datos iniciales, tratamos de ver cuando va a frenar el crecimiento poblacional de individuos contagiados o muertos, lo cual tiene una utilidad obvia. Una limitación que nos vamos a encontrar es que, si bien la predicción será muy sólida en una curva monótona, si la infección es no monotónica, necesitaremos introducir flexibilidad en los parámetros. Otra limitación es que se acumulen variables, como por ejemplo, que haya un parámetro estacional o temporal en la infección. Para ello, se podría usar el modelo de von Bertalanffy ajustado a una variación periódica, o bien usar un modelo bifásico. Veremos en otra parte, cuando apliquemos este conocimiento a COVID-19, que sería interesante estudiar esto, pues probablemente estemos en un crecimiento bifásico. Nosotros aquí no vamos a introducir complicaciones, que superan el objetivo de una entrada de blog, y vamos a mantener la máxima simplicidad. Una vez realizado el ajuste no lineal, evaluaremos su calidad con algunos parámetros, incluyendo el gráfico de residuales, que nos puede aportar información interesante sobre el modelo, como identificar variaciones temporales superpuestas al crecimiento.

El análisis de los datos

Si representamos el número de fallecidos acumulados, es decir, el número total de fallecidos desde el inicio de la recogida de datos, obtenemos la siguiente gráfica:

El gráfico muestra los datos (puntos negros) y el ajuste al modelo de Richards del crecimiento (línea azul). Del ajuste obtenemos los parámetros:

  • Máximo: 262 fallecidos.
  • En 22 días desde que se inició la infección se alcanzó el punto de inflexión y el número de fallecidos/dia (rate) comienza a frenar.
  • Factor de asimetría 0.1
  • Coeficiente de crecimiento: 0.047. Esto es particularmente importante, por que la COVID-19 nos vamos a encontrar un coeficiente de crecimiento de 0.037 a 0.039.

Ahora, a partir de estos datos de fallecidos acumulados, podemos calcular el rate, o número de fallecidos/día y compararlo con la predicción teórica, obtenida derivando la curva azul de la gráfica de arriba. Obtenemos la característica curva de campana:

En ésta curva podemos comparar la predicción teórica (curva azul) según el modelo de Richards y, en negro, los datos reales. La predicción con la curva teórica sugería que el máximo se alcanzaría a los 18 días del inicio de la infección, aunque en los datos reales se alcanza a los 22. Esta discrepancia puede deberse a las irregularidades en la toma de los datos. En ésta curva, nos acercamos mucho a cero fallecidos por día a partir del día 70 desde el inicio. En una epidemia es esencial que ésta curva de campana se “aplaste” (algo que hemos oído mucho en la prensa). Esto es esencial para no saturar los servicios de emergencia y no provocar colapso social. Obsérvese que la predicción en azul está en buen acuerdo con los datos reales y ambas llegan a cero aproximadamente al mismo tiempo. Esto, que se comprueba con una prueba de predicción, es esencial para estimar cuándo podemos esperar bajar de un determinado número de casos por día.

Aún podemos ir más allá: podemos estudiar las variaciones diarias en el rate de casos. Esto es, como está acelerando o frenando la epidemia, o cómo están acelerando la cantidad de fallecidos por día, o cómo está frenando la cantidad de fallecidos por día y compararlo con la predicción teórica. Esto se obtiene mediante la derivada segunda del modelo de Richards. La derivada segunda es muy sensible a las variaciones del rate, y nos permite detectar pequeños cambios que podemos correlacionar con diferentes variables (datos ambientales, eventos durante el proceso de recogida de datos, etc. La gráfica de la derivada segunda tiene ésta forma:

Cuando la gráfica se encuentra en la zona positiva, la infección está acelerando. Cuando está en la zona negativa, está frenando. Como se ve, a partir del dia 18 en la predicción teórica o el 22 en los datos reales, la epidemia entró en la fase de frenada. En el dia 80 se llegan a los cero fallecidos y se puede considerar que la infección se ha detenido. Los datos reales están sujetos a fuertes variaciones, que pueden responder a muchas razones. La más evidente, son las variaciones en el proceso, desde que muere un paciente, se registra y se realiza el proceso burocrático y se contabiliza.

Evaluación del modelo

Hay varios modos de evaluar la bondad del ajuste de los datos al modelo no lineal teórico. El primero de ellos es el valor de R-cuadrado, que nos da un valor de R1=0.9980. Aquí, en ajustes a ecuaciones no lineales hay una diferencia fundamental con la regresión lineal. En ésta última, el valor de R-cuadrado nos informa sobre la linealidad de los datos (matemáticamente lineal, R2=1). Sin embargo, en regresión no lineal, un valor de R-cuadrado elevado no indica necesariamente que hay un buen ajuste al modelo. Por ello, éste parámetro debe usarse con precaución.

Otra forma sencilla de evaluar el modelo es mediante el análisis de los residuales. La gráfica de residuales (distancias de los puntos experimentales al modelo teórico propuesto) nos dará una distribución de puntos alrededor de una línea central. Si los datos experimentales se ajustaran exactamente, matemáticamente, al modelo, todos los puntos estarían en la línea central. Esto en datos reales es imposible, pues toda toma de datos experimentales tiene un rango de error y variabilidad. Así que, si el ajuste es bueno, esperamos que los residuales se distribuyan siguiendo una distribución gaussiana alrededor de la línea central. Visualmente, lo que veremos es una nube de puntos en torno a una línea central:

Esto no se hace sólo visualmente, sino que realizamos un test de hipótesis, el test de normalidad de D’Agostino y Pearson, que responde a ésta pregunta: Si el modelo que hemos usado para explicar los datos es correcto y simulamos datos que se han esparcido alrededor de la línea central siguiendo una distribución normal, ¿cual es la probabilidad de que estos datos simulados se desvíen tanto de la distribución normal como se desvían los datos reales?. Si esta probabilidad es elevada, quiere decir que podemos aceptar la hipótesis de que los datos experimentales se ajustan correctamente al modelo de Richards. En nuestro caso es así, ya que la probabilidad del test D’Agostino nos da un valor P=0.29, que es muy alto. Tan sólo hay que señalar que en la fase final, a partir del dia 60, tenemos un patrón, que puede ser fácilmente explicable por los motivos que expusimos antes.

Así, concluimos que con el modelo de Richards de cinco parámetros, de curva logística generalizada asimétrica, tenemos un buen modelo de crecimiento de la epidemia. Este modelo lo hemos aplicado a COVID-19, y en otra entrada explicaré los resultados hasta ahora. La aplicación de éste modelo, que he ido exponiendo en mi cuenta de Facebook, nos ha permitido predecir que a partir del 8-10 de mayo la velocidad de fallecimientos por COVID-19 bajará de los 50 casos/día y que hacia mediados de mayo entraremos en la asíntota de la fase final del proceso, salvo sorpresas. Recomiendo a quien tenga interés consultar en mi cuenta, para ir viendo cómo ha ido evolucionando el modelo y las interesantes observaciones que estamos viendo en éstos últimos días.

Referencias consultadas

  • https://besjournals.onlinelibrary.wiley.com/doi/pdf/10.1111/j.2041-210X.2012.00231.x
  • scihub.unblockit.one/10.1016/j.cam.2016.02.034
  • https://www.researchgate.net/publication/335704097_La_epidemia_de_gripe_de_1918_en_Alcala_la_Real?fbclid=IwAR2SskBVyzkpvm4sZmd78rCJ-4SBAs4blVlsAB0jG1qV1ot7_M1DOlw0Bo4
  • (1) Szparaga, A.; Kocira, S. Generalized Logistic Functions in Modelling Emergence of Brassica Napus L. PLoS One 2018, 13 (8). https://doi.org/10.1371/journal.pone.0201980.
  • (1) Román-Román, P.; Serrano-Pérez, J. J.; Torres-Ruiz, F. A Note on Estimation of Multi-Sigmoidal Gompertz Functions with Random Noise. Mathematics 2019, 7 (6), 1–18. https://doi.org/10.3390/MATH7060541.
  • (1) Aggrey, S. E. Comparison of Three Nonlinear and Spline Regression Models for Describing Chicken Growth Curves. Poult. Sci. 2002, 81 (12), 1782–1788. https://doi.org/10.1093/ps/81.12.1782.
  • (1) Kaplan, S.; Gürcan, E. K. Comparison of Growth Curves Using Non-Linear Regression Function in Japanese Quail. J. Appl. Anim. Res. 2018, 46 (1), 112–117. https://doi.org/10.1080/09712119.2016.1268965.
  • (1) Moler, C. Experiments with MATLAB. MathWorks, Co 2011, 278. https://doi.org/http://dx.doi.org/10.1017/CBO9780511813887.016.
  • (1) Szparaga, A.; Kocira, S. Generalized Logistic Functions in Modelling Emergence of Brassica Napus L. PLoS One 2018, 13 (8). https://doi.org/10.1371/journal.pone.0201980.
  • Akaike, H. 1974. A new look at the statistical model identification. IEEE Trans. Automat. Control, 19, 716-723
  • Janoschek, A. 1957. Das reaktionskinetische Grundgesetz und seine Beziehungen zum Wachstums- und Ertragsgesetz. Stat. Vjschr. 10, 25-37.
  • Sager, G. 1978 Zuwachsfunktionen vom Typ dw/dt=k·t^p(E-W)^n und ihre Integrale. Anat. Anz. 144, 366-374.
  • Sengul, T & Kiraz, S. 2005. Non-linear models for growth curves in large white turkeys. Turk J. Vet. Anim. Sci. 29, 331-337.
  • Somer, I.F. 1988. On a seasonally oscillating growth function. Fishbyte 6(1):8-11.
  • Winsor, C. 1932. The Gompertz curve as a growth equation. Proc. Nat. Acad. Sciences, 18, 1-8.
  • (1) Zeide, B. Analysis of Growth Equations. For. Sci. 1993, 39 (3), 594–616. https://doi.org/10.1093/forestscience/39.3.594.
  • (1) Martinez, A. S.; González, R. S.; Terçariol, C. A. S. Continuous Growth Models in Terms of Generalized Logarithm and Exponential Functions. Phys. A Stat. Mech. its Appl. 2008, 387 (23), 5679–5687. https://doi.org/10.1016/j.physa.2008.06.015.
  • DʼAgostino RB, Belanger A, DʼAgostino Jr RB. A Suggestion for Using Powerful and Informative Tests of Normality. The American Statistician. 1990; 44(4):316-21.
  • DʼAgostino R, Pearson ES. Tests for Departure from Normality. Empirical Results for the Distributions of b 2 and √b 1. Biometrika. 1973; 60(3):613-22.

Leave a Reply

Your email address will not be published.

sixteen − 11 =

This site uses Akismet to reduce spam. Learn how your comment data is processed.