Resolución de EDOs para Modelar Epidemias con Python

En la clase pasada vimos que un sistema de EDOs para describir la propagación de una epidemia está dado, en su forma más simple por:

Donde β es una constante cuyo significado revelaremos a continuación. Para poder resolver estas ecuaciones y obtener las funciones que nos permitirán predecir cómo evolucionará la epidemia usaremos métodos muy populares en matemáticas. Puedes profundizar en este curso de ecuaciones diferenciales aquí. En nuestro caso particular usaremos trucos poco comunes pero necesarios para esta situación específica.

Revisemos la solución en un paso a paso detallado:

1. Suponemos que la población total N es constante por lo que siempre se cumple que S(t) + I(t)=N. De esto despejamos así S(t) = N - I(t) y lo usamos para obtener una sola ecuación con una sola incógnita.

2. De lo cual la segunda ecuación diferencial se convierte en:

3. La ecuación que acabamos de obtener pertenece a una clase especial de ecuaciones denominada ecuaciones diferenciales de Bernoulli, y para estas podemos definir una nueva variable así:

4. Donde el cambio para la derivada se traduce en y'(t) = -I'(t)/I(t)². Ahora, reemplazando todo lo anterior en la ecuación [3] nos permite obtener la siguiente ecuación lineal:

5. La cual se resuelve fácilmente con el método de separación de variables para obtener una solución con la forma:

6. Donde, A es una constante que se determina por la condición inicial del problema. En nuestro caso la condición inicial problema es la cantidad de personas infectadas al inicio de nuestro sistema. Lo usual es definir que al principio siempre existe al menos 1 persona infectada, y esto lo podemos definir así I(t) = 1

7. Al reemplazar la condición inicial en la ecuación [5] obtenemos ahora:

8. Finalmente, de la ecuación [6] podemos obtener también la función para el número de individuos susceptibles con la relación:

Y así, hemos obtenido funciones que nos permiten saber cómo evolucionan las variables del sistema dinámico para cualquier valor del tiempo. Ahora, igual que hicimos con el ejercicio del péndulo, podemos escribir código en Python que nos permita visualizar las predicciones que se derivan de estas funciones.

Primero importamos las librerías necesarias y definimos los parámetros y funciones de nuestro modelo:

Así, hemos definido la población total con 100 personas y un solo individuo inicialmente infectado, mientras que el valor de β lo hemos definido de forma arbitraria, por ahora. Luego creamos listas de elementos que constituyen las tablas de datos sobre las cuales graficamos nuestras funciones:

En el ejemplo vemos que estamos tomando el tiempo desde 0 hasta 49 días y así mismo podemos ver la otra lista con el número de infectados total, día a día. Aquí lo extraño es ver que un día tenga 1.80 infectados y estos números no deben preocuparnos pues son simplemente aproximaciones que se derivan en que estamos usando una variable continua para aproximar un fenómeno que es discreto. Así por ejemplo, que en el día 2 haya 1.80 infectados y luego en el día 3 haya 2.42 infectados nos indica que aproximadamente le toma a la epidemia entre dos y tres días alcanzar un total de dos personas infectadas. Ahora, una visualización completa del comportamiento en el tiempo la obtenemos así:

Donde el eje horizontal es el tiempo en días y el eje vertical es el número de individuos. Varios hechos observamos de este gráfico:

1. Después de aproximadamente 30 días, toda la población resulta afectada.
2. El número de infectados siempre aumenta porque no estamos asumiendo dentro del modelo la posibilidad de recuperación de los individuos.

Como ves, las predicciones que se derivan de este modelo son algo obvias pero totalmente coherentes a la luz de las simplificaciones que incluimos en nuestro modelo.

Ahora te propongo un ejercicio adicional, calculando las soluciones de nuestro sistema para diferentes valores del parámetro β como se ve a continuación:

En este caso solo usamos las curvas de población de infectados. Lo que vemos es que a medida que el parámetro β crece , la curva alcanza el máximo en menos tiempo, esto indica que el parámetro β es una medida de la rapidez de propagación de la epidemia y ese es el modo en que debemos interpretar este número.

El hecho de que podamos predecir el valor de una variable para cualquier instante de tiempo futuro o pasado es lo que entendemos por un sistema determinístico, ya que su evolución está completamente determinada por una regla claramente definida. Es así, como el lenguaje de las ecuaciones diferenciales, aplicadas a estas situaciones, nos permite codificar el concepto de determinismo en números concretos. En este link puedes ver el notebook completo con el código.