Las matemáticas de la epidemia, Parte I

En un año normal, cualquiera diría que hablar de epidemias y su difusión no es exactamente un tema de "vida cotidiana", por lo menos en la rica y medicalizada Europa (lo es, lamentablemente, en muchos países de Africa, Asia o América, por poco que hablemos de ello). Pero el 2020 que acabaos de despedir ha sido diferente: ha sido el año en que las epidemias han entrado prepotentemente en nuestra vida, el año en que todos hemos hablado de IA, de R0, de curvas ascendentes, de doblar las curvas, etc.

La pandemia de covid ha afectado no sólo, y de manera trágica, la vida de muchos, sino también el discurso público. De repente, España (y el resto del mundo) se ha llenado de personas que querían saber más sobre las epidemias, sus causas y su manera de extenderse. España se ha llenado también, y lamentablemente, de personas que sentenciaban mucho y sabían muy poco, pero este no es el lugar para hablar de ellos.

Las matemáticas han estudiado los fenómenos de difusión viral desde hace más de 40 años. Al principio estos estudios se limitaban a la difusión de virus y enfermedades pero luego, como es a menudo el caso de los modelos matemáticos, se ha notado que muchos fenómenos aparentemente alejado siguen las mismas leyes y hoy en día estos modelo se usan también para estudiar, por ejemplo, la difusión de opiniones o la difusión de bulos en las redes sociales. Aquí nos interesa la aplicación de estos modelos a la difusión de virus, y nos limitaremos a los más sencillos.

Antes de empezar, una recomendación: los modelos matemáticos no son la realidad sino, en el mejor de los casos, una versión muy simplificada de la realidad.Un buen modelo matemáticos nos puede dar indicación preciosas sobre el comportamiento del fenómeno real, pero siempre hay que tener cuidado a la hora de tomar decisiones o expresar opiniones basándose en un modelo matemáticos. Antes de hacerlo, es necesario por lo menos saber las simplificaciones que se han hecho, su justificación y los posibles errores que se pueden introducir.

Y ahora el modelo. Nos interesa la evolución de la epidemia a lo largo del tiempo. Mediremos el tiempo en días, y usaremos la variable t para indicar el día en que nos encontramos. La variable asume valores $t=0,1,2,...,100$, es decir, seguiremos el comportamiento de la epidemia durante 100 días. Por cada día seguiremos cuatro magnitudes:

• $S(t)$: el número de personas "susceptibles" el día $t$. Los susceptibles son las personas no contagiadas o que se han recuperado pero no tienen anticuerpos. Es decir, las personas que se pueden contagiar del virus.
• $I(t)$: el número de personas "infectadas" el día $t$. Las personas enfermas.
• $R(t)$: el número de personas "recuperadas" el día $t$. Se trata de las personas que han superado la enfermedad y que todavía tienen anticuerpos, es decir, las personas que el día t no están enfermas y no se pueden infectar.
• $V(t)$: el número de víctimas, es decir, el número de personas que en este día han sucumbido a la enfermedad.

También consideraremos un número total N de personas fijo, que no cambia durante nuestro análisis (esto quiere decir que el intervalo en que analizamos la epidemia es lo suficientemente corto para que el nacimiento de niños y la muerte de personas por causas diferentes de la covid no son relevantes). En cualquier momento una persona puede ser susceptible, infectada, recuperada o ser una víctima: cada persona se encuentra en una de estas categoría, y ninguna está en dos categoría, por tanto cada día t: \[ S(t) + I(t) + R(t) + V(t) = N \]

La base del modelo es la llamada "predicción a un paso". Supongamos que conocemos los valores $S(t)$, $I(t)$, $R(t)$, $V(t)$, el día $t$. ¿Cuáles son los valores $S(t+1)$, $I(t+1)$, $R(t+1)$, $V(t+1)$ para el día siguiente? Empecemos con $S(t)$. Podemos escribir:

\[ S(t+1) = S(t) - A + B \]

En esta ecuación $A$ es el número de personas susceptibles que el día $t$ se han infectado (y que por tanto el día $t+1$ ya no son susceptibles), $B$ es el número de personas que el día $t$ eran recuperados pero que al día $t+1$ han perdido los anticuerpos y por tanto están susceptibles.

El número $A$ depende de los contactos entre personas susceptibles e infectadas. No podemos saber cuantos son exactamente estos contactos, pero los podemos estimar estadísticamente. Un modelo muy utilizado es el que los físicos llaman "de moléculas de gas". Si asumimos que las personas se mueven mucho y cada persona puede encontrar a cualquier otra persona, entonces el número $A$ es proporcional a $S(t)\cdot{I(t)}$. Introducimos un parámetro $\beta$, que representa la probabilidad que un encuentro resulte en un contagio. Con esto podemos escribir

\begin{equation} A = \beta \frac{S(t)I(t)}{N}\rule{10em}{0pt}(1) \end{equation}

(El factor $1/N$ es un factor de normalización para que la ecuación describa una probabilidad de encuentros.) El valor $B$ depende de cuantas personas se han recuperado en ese día. Si hay $R(t)$ personas recuperadas y llamamos $\gamma$ la fracción de personas recuperadas que pierden los anticuerpos cada día, tenemos

\begin{equation} B = \gamma R(t)\rule{10em}{0pt}(2) \end{equation}

Con estas definiciones podemos escribir la ecuación que regula la evolución del número de personas susceptibles:

\begin{equation} S(t+1) = S(t) - \beta \frac{S(t)I(t)}{N} + \gamma R(t) \end{equation}

Consideremos ahora las personas infectadas. En este caso podemos escribir

\begin{equation} I(t+1) = I(t) + A - C - D \end{equation}

La variable $A$ son las personas que se han infectado ese día, y ya hemos calculado su valor en la ecuación (1). La variable $C$ son las personas que ese día se han recuperado, y $D$ son las personas que han muerto ese día. Así como hemos hecho antes, podemos asumir que cada día una fracción dada de las personas infectadas se recuperan, y otra fracción dada muere, por tanto escribiremos

\begin{equation} \begin{array}{rcll} C &= &\rho I(t) &\rule{10em}{0pt}(3) \\ D &= &\zeta I(t) &\rule{10em}{0pt}(4) \end{array} \end{equation}

Con esto obtenemos la ecuación

\begin{equation} I(t+1) = I(t) + \beta \frac{S(t)I(t)}{N} - \rho I(t) - \zeta I(t) \end{equation}

Tenemos ahora las personas que se recuperan. En este caso escribimos

\begin{equation} R(t+1) = R(t) - B + C \end{equation}

La variable $B$, definida en la ecuación (2) son las personas que ese día han vuelto a ser susceptibles, mientras $C$, definida en (3) son las personas infectadas que ese día se han recuperado. Por tanto podemos escribir

\begin{equation} R(t+1) = R(t) - \gamma R(t) + \rho I(t) \end{equation}

Finalmente, para las víctimas, escribimos:

\begin{equation} V(t+1) = V(t) + \zeta I(t) \end{equation}

Las víctimas de hoy son las víctima de ayer más las personas que han muerto entre ayer y hoy. Tenemos ahora nuestro modelo completo:

\begin{equation} \begin{array}{rcl} S(t+1) &= &S(t) - \beta \frac{S(t)I(t)}{N} + \gamma R(t) \\ I(t+1) &= &I(t) + \beta \frac{S(t)I(t)}{N} - \rho I(t) - \zeta I(t) \\ R(t+1) &= &R(t) - \gamma R(t) + \rho I(t) \\ V(t+1) &= &V(t) + \zeta I(t) \end{array} \end{equation}

Estas ecuaciones tienen una característica que nos ayuda a tener confianza de su validez. Se puede demostrar que

\begin{equation} S(t+1)+I(t+1)+R(t+1)+V(t+1) = S(t)+I(t)+R(t)+V(t) \end{equation}

Es decir, las ecuaciones nos confirman que el número total de personas (es decir, N) no cambia: las personas pasan de una categoría a otra sin alterar su número total.

Para completar el modelo es necesario establecer el valor de las constantes β, γ, ρ, ζ. Para esto se pueden usar los datos efectivos de la epidemias (tasa de crecimiento, duración media de la enfermedad, mortalidad, etc.). El procedimiento para conseguir estos parámetros no es muy complicado pero en este momento nos llevaría demasiado lejos. Si usamos como datos los relativos a la primera fase de la epidemia en España, conseguimos los valores siguientes:

\begin{equation} \begin{array}{rcl} \beta &= &0.48 \\ \gamma &= &0 \\ \rho &= &0.05 \\ \zeta &= &0.06 \end{array} \end{equation}

Notamos que se ha puesto γ=0. ¿Qué quiere decir esto? Si nos fijamos en las ecuaciones, γ=0 quiere decir que ninguna persona recuperada volverá a ser susceptible. Esto no es completamente cierto, pero en la primera fase de la epidemia (de Marzo a Junio) no se confirmó en España ningún caso de reinfección, y esto resulta en nuestra elección de γ=0. ¿Qué resultados nos dan nuestras ecuaciones? No es posible calcular soluciones exactas, pero es muy fácil introducirle en un ordenador y efectuar los cálculos para ver que curvas de contagio nos dan. Nos interesan sobre todo dos cantidades: el número de casos nuevos diarios C(t)=I(t)-I(t-1) y el número de víctimas diarias M(t)=V(t)-V(t-1). Para arrancar el cálculo, asumimos que el primer día (t=0) todo el mundo es susceptible excepto una persona que se ha infectado (el "paciente cero"). Es decir:

\begin{equation} \begin{array}{rcl} S(0) &= &N-1 \\ I(0) &= &1 \\ R(0) &= &0 \\ V(0) &= &0 \end{array} \end{equation}

Si insertamos nuestras formulas en un sencillo programa de ordenador poniendo, por ejemplo, $N=10.000$, los resultados son los de Figura 1.a, donde vemos en rojo el número de infectados cada día y en negro el número de víctimas diarias.

El número de infectados (lo que hemos llamado $I(t)$) no es lo mismo que el número de nuevos casos diarios; esto lo tenemos en la Figura 2. El comportamiento es muy parecido a lo que observamos en los datos reales, sobre todo en la fase ascendente (el descenso es más lento en nuestro caso que en el caso real, en cuanto nuestro modelo, por el momento, no tiene en cuenta un elemento muy importante: el confinamiento del 14 de Marzo). A continuación me centraré sobre todo en el número de casos (Figura 1), dado que es el que determina la saturación de los hospitales.

Se trata claramente de una versión muy idealizada: los datos reales no tienen esta regularidad, principalmente porque el parámetro $\beta$, que depende de los contactos entre las personas, no es exactamente constante: de un día para otro los contactos entre las personas varian de manera imprevisible.

Los modelos matemáticos no nos dan una imagen completa de la difusión de la epidemia, pero son muy útiles para estudiar sus características generales y ver, a grandes líneas, el efecto de varias estrategias y comportamiento. Para esto,es mejor evitar el ruido que hace las curvas más realistas y trabajar con la curva idealizada, que nos muestra mejor los rasgos generales sin confundirnos con los detalles que pueden variar imprevisiblemente. Hablaré de ello en el próximo artículo.