Modelo SIR, modelo epidemiológico, con Python

Hoy vamos a hablar de ecuaciones diferenciales (y de modelos epidemiológicos). Imaginad que existe una pandemia mundial de una nueva enfermedad que llamaremos X. Vamos a establecer una serie de ideas para poder trabajar en un modelo de simulación de estas pandemias.

El modelo SIR

Voy a empezar hablando un poco del modelo SIR. Según la Wikipedia:

El modelo SIR es uno de los modelos epidemiológicos más simples capaces de capturar muchas de las características típicas de los brotes epidémicos. El nombre del modelo proviene de las iniciales S (población susceptible), I (población infectada) y R (población recuperada). El modelo relaciona las variaciones de las tres poblaciones (Susceptible, Infectada y Recuperada) a través de la tasa de infección y el período infeccioso promedio.

Este modelo lo formularon, inicialmente, William Ogilvy Kermack y A. G. McKendrick y lo publicaron en 1927 en los proceedings of the Royal Society.

La población susceptible, \(S\), es aquella que no ha pasado el virus y no tiene inmunidad al agente infeccioso. Esta población se puede infectar si está en contacto con el agente infeccioso.

La población infectada, \(I\), es aquella que ya tiene el agente infeccioso en su cuerpo y pueden, por tanto, infectar a miembros de la población susceptible.

La población recuperada, \(R\), es la población que ha sufrido el agente infeccioso pero ya no lo va a sufrir más porque se ha recuperado y ha generado inmunidad y no lo transmite o porque ha muerto. Estos ya no pueden infectar a otros.

Definición del modelo matemáticamente

Lo que tenemos es lo siguiente:

\(S \Rightarrow I \Rightarrow R\)

Es decir, un miembro de la población \(S\) puede pasar a ser miembro de la población \(I\), si está en contacto con el agente infeccioso, y, de ahí, pasará a ser miembro de la población \(R\).

El modelo dispone de varios parámetros que controlan cómo se infectan y recuperan los miembros en las distintas poblaciones:

  • \(\beta\): El parámetro beta nos indica cómo de sencillo es que un miembro de la población \(I\) infecte a un miembro de la población \(S\).

Como vamos a usar probabilidades usaremos el parámetro beta durante un intervalo de tiempo como la probabilidad de que un miembro de \(I\) infecte a un miembro de \(S\) en un intervalo de tiempo \(t\): \( \beta \Delta t \)

  • \( \gamma \): el parámetro gamma nos indica la facilidad o dificultad de que un miembro de la población \(I\) se recupere, es decir, pase a la población \(R\) (con inmunidad o sin vida).

Al igual que antes, la probabilidad de que un miembro de \(I\) se recupere, pase a \(R\), en un intervalo de tiempo \(t\) será: \( \gamma \Delta t \)

Puedes ver que tienes tres clases, tres posibilidades para que un miembro esté en una población. Las tres clases conforman la población total, \(N = S + I + R\). Por otro lado, tienes dos parámetros que te indican las probabilidades de que un individuo se mueva de una clase a otra. Con esto vamos a ver si derivamos una serie de ecuaciones para resolver el problema.

Los susceptibles en un momento de tiempo posterior, \( S(t + \Delta t) \), serán igual a los susceptibles en este momento, \( S(t) \), menos (salen del grupo \(S\)) los posibles susceptibles que han sido infectados entre este momento y el momento de tiempo posterior, \( \beta \Delta t \frac{S(t)}{N} I(t) \):

\( S(t + \Delta t) = S(t) – \beta \Delta t \frac{S(t)}{N} I(t) \)

Los infectados en un momento de tiempo posterior, \( I(t + \Delta t) \), serán igual a los infectados en este momento, \( I(t) \), más (entran en el grupo \(I\)) los susceptibles que han sido infectados entre este momento y el momento de tiempo posterior, \( \beta \Delta t \frac{S(t)}{N} I(t) \), menos (salen del grupo \(I\)) los infectados que se recuperan entre este momento y el momento de tiempo posterior, \( \gamma \Delta t I(t) \):

\( I(t + \Delta t) = I(t) + \beta \Delta t \frac{S(t)}{N} I(t) – \gamma \Delta t I(t) \)

Por último, los recuperados en un momento de tiempo posterior, \( R(t + \Delta t) \), serán igual a los recuperados en este momento, \( R(t) \), más (entran en el grupo \(R\)) los infectados que se recuperan entre este momento y el momento de tiempo posterior, \( \gamma \Delta t I(t) \):

\( R(t + \Delta t) = R(t) + \gamma \Delta t I(t) \)

Estas ecuaciones que acabamos de derivar las vamos a reorganizar para que sean más familiares a lo que es una derivada:

\( \frac{S(t + \Delta t) – S(t)}{\Delta t} = – \beta \frac{S(t)}{N} I(t) \)

\( \frac{I(t + \Delta t) – I(t)}{\Delta t} = \beta \frac{S(t)}{N}) I(t) – \gamma I(t) \)

\( \frac{R(t + \Delta t) – R(t)}{\Delta t} = \gamma I(t) \)

Si ahora hacemos que los pasos temporales sean muy pequeños, \( \lim\limits_{\Delta t \to 0} \frac{f(t + \Delta t) – f(t)}{\Delta t}\), llegamos al concepto de derivada que mide la variación de la función (en este caso \(S\), \(I\) o \(R\)) cuando hay una pequeña variación de \(t\). Las ecuaciones diferenciales que nos quedan son:

\( \frac{dS}{dt} = – \frac{\beta S I}{N} \)

\( \frac{dI}{dt} = \frac{\beta S I}{N} – \gamma I \)

\( \frac{dR}{dt} = \gamma I \)

Con esas ecuaciones diferenciales nos queda un modelo relativamente simple para poder resolver el problema de conocer la evolución de la pandemia.

Condiciones iniciales

El modelo necesitará unas condiciones iniciales, condiciones de \(S\), \(I\) o \(R\) en el momento \(t=0\), desde las que habrá que hacerlo evolucionar. Las condiciones iniciales podrían ser algo como lo siguiente:

La población susceptible inicial será toda la población, \(N\), menos una pequeña parte de individuos que poseen el agente infeccioso, \(n\), menos una serie de individuos que se han recuperado (empezamos con 0):

\( S(t=0) = N – n – 0 \)

La población de infectados iniciales son los que son capaces de desencadenar la pandemia y empieza por una parte muy modesta de la población total (\( n << N \)):

\( I(t=0) = n \)

Como estamos al inicio de una posible pandemia el número inicial de gente recuperada es cero en este momento, \( t = 0 \):

\( R(t=0) = 0 \)

Además, asumimos, para simplificar, que no hay nacimientos, es decir, que la población no crece:

\( \frac{dS}{dt} + \frac{dI}{dt} + \frac{dR}{dt} = 0 \)

De lo anterior se deduce que la población la consideramos constante:

\( S(t) + I(t) + R(t) = cte = N \)

Tasa básica de reproducción

En epidemiología existe un valor que se conoce como tasa básica o ratio básico de reproducción. Se expresa como \(R_0\). Se puede ver como el número promedio de casos nuevos que genera un caso dado a lo largo de un período infeccioso considerando que toda la población es susceptible.

Vamos a intentar derivar \(R_0\) para el caso del modelo SIR usando una serie de suposiciones y la ecuación diferencial para \(I\). Partiendo de:

\( \frac{dI}{dt} = \frac{\beta S I}{N} – \gamma I \)

Y considerando que \(S >> I\) y que \(S \approx N\) la ecuación anterior podría quedar similar a:

\( \frac{dI}{dt} = \frac{\beta N I}{N} – \gamma I = (\beta – \gamma) I\)

Si resolvemos la anterior ecuación diferencial obtenemos:

\( \int_{I_0}^I \frac{dI}{I} = \int_0^t (\beta – \gamma) dt \)

\( I(t) = I_0 e^{(\beta – \gamma)t} \)

La ecuación anterior nos vendría a decir como evoluciona la población de infectados. Si la exponencia es positiva la función sería creciente mientras que si es negativa sería decreciente:

\( \beta – \gamma > 0, \nearrow \)

\( \beta – \gamma < 0, \searrow \)

Si dividimos lo anterior por \(\gamma\) y despejamos llegamos a:

\( \frac{\beta}{\gamma} > 1, \nearrow \)

\( \frac{\beta}{\gamma} < 1, \searrow \)

El ratio que encontramos entre los parámetros \(\beta\) y \(\gamma\) sería el \(R_0\). Si el \(R_0\) está por encima de 1 es mala cosa porque están aumentando los casos. Si, por el contrario, está por debajo de 1 indica que se está aplanando la curva (estas palabras se han estado escuchando bastante durante los últimos meses).

Programando

Como aquí solemos hablar de programación y de resolver problemas científicos usando programas vamos a meternos en harina para resolver este problema.

En el pasado se ha hablado bastante de ecuaciones diferenciales:

Por lo que no voy a entrar en muchos detalles. Podéis echar un ojo a los anteriores artículos para saber más.

Todo el código que figura a continuación está basado, con pequeñas adaptaciones, en lo que puedes encontrar en este enlace.

Como siempre, primero importamos una serie de cosas que vamos a usar:

Definimos la función donde resolveremos el sistema de ecuaciones diferenciales del modelo SIR:

Definimos una función para pintar los resultados:

Definimos algunas cosas:

Con el siguiente código vamos a usar todo lo anterior descrito hasta ahora y veremos los resultados. Para ello usamos odeint. Puedes ver información de las posibilidades de odeint en la documentación oficial. Básicamente, le pasamos nuestro sistema de ecuaciones diferenciales, deriv, las condiciones iniciales, y0, los puntos donde resolver el sistema de ecuaciones diferenciales, t, y otras condiciones que definen el sistema de ecuaciones diferenciales, N, beta, gamma.

Puedes tocar todos estos parámetros y ver cómo afecta cada uno de ellos al resultado final: N, I0, R0, beta, gamma.

Modificando los recuperados

Vamos a añadir una pequeña modificación para mostrar la gente que moriría en función de una tasa que definamos. Imaginemos que muere el 5% de los individuos que se infectan. Voy a modificar la función plot para reflejar esto:

Notas finales

Si le ves alguna utilidad a todo esto que hemos visto hoy nos puedes mandar un comentario o, incluso, puedes publicar tus propios artículos en pybonacci si no dispones de otro espacio para ello.

Nota de descargo: no soy biólogo, médico ni nada por el estilo por lo que si he metido alguna gamba y me lo hacéis saber os lo agradezco.

Saludos.

Deja una respuesta

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *

ninety four − 86 =