Coronavirus
Durante los últimos días hemos visto como el Coronavirus se ha comenzado a esparcir exponencialmente dentro de la población mundial. En el siguiente post buscamos simular bajo el modelo SIR, el efecto que tendría a lo largo del tiempo.
El modelo SIR
Una descripción simple de la capacidad de contagio de una enfermedad en una población está formulado en el modelo SIR, que divide la población N, entre tres compartimentos que cambian en función del tiempo, t.
- S(t) son los susceptibles pero aún no infectados con la enfermedad;
- I(t) son el número de individuos infectados
- R(t) son aquellos que se han recuperado del virus y ahora son inmunes.
El modelo describe el cambio en la población en cada uno de estos grupos en términos de dos parámetros, β y γ.
Gamma es la tasa de recuperación, dado esto, 1/gamma es la media del periodo donde un individuo puede contagiar a otro.
Las ecuaciones que describen el modelo, fueron derivadas por Kermack y McKendrick en el 27:
El siguiente código integra las ecuaciones para una enfermedad caracterizada con un B = 0.2, 1/gamma = 10 días, en una población de N = 1000. El modelo comienza con solamente una persona contagiada en el día 0: I(0) = 1. Las curvas representan los grupos a lo largo del tiempo.
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
# Poblacion total, N
N = 1000
# Numero inicial de individuos infectados y recuperados, I0 y R0.
I0, R0 = 1, 0
# Todo el resto es suceptible a la infección inicialmente.
S0 = N - I0 - R0
# Tasa de contacto, beta, tasa de recuperación y gamma en 1/días
beta, gamma = 0.2, 1./10
# Grilla de tiempo, en días.
t = np.linspace(0, 160, 160)
# Ecuaciones
def deriv(y, t, N, beta, gamma):
S, I, R = y
dSdt = -beta * S * I / N
dIdt = beta * S * I / N - gamma * I
dRdt = gamma * I
return dSdt, dIdt, dRdt
# Vector condiciones iniciales
y0 = S0, I0, R0
# Integrar las ecuaciones en el tiempo t.
ret = odeint(deriv, y0, t, args=(N, beta, gamma))
S, I, R = ret.T
# Plotear Data
fig = plt.figure(facecolor='w')
ax = fig.add_subplot(111, axisbelow=True)
ax.plot(t, S/N, 'b', alpha=0.5, lw=2, label='Suceptibles')
ax.plot(t, I/N, 'r', alpha=0.5, lw=2, label='Infectados')
ax.plot(t, R/N, 'g', alpha=0.5, lw=2, label='Recuperados con inmunidad')
ax.set_xlabel('Tiempo (días)')
ax.set_ylabel('% Población')
ax.set_ylim(0,1.2)
ax.yaxis.set_tick_params(length=0)
ax.xaxis.set_tick_params(length=0)
ax.grid(b=True, which='major', c='w', lw=2, ls='-')
legend = ax.legend()
legend.get_frame().set_alpha(0.5)
for spine in ('top', 'right', 'bottom', 'left'):
ax.spines[spine].set_visible(False)
plt.savefig('fig1.png', dpi = 300)
plt.show()
Resultados
Bajo el análisis de casi medio año, y con los parámetros definidos, encontraríamos el punto máximo de contagio en unos 65 días. Si consideramos el primer caso del Coronavirus el 31 de Diciembre, el máximo esperado de contagio estaría dado a mediados de Marzo, llegando en torno al 18% de la población – es decir – a 1,35 billones de personas. En perspectiva, la H1N1 llegó a más de 500 Millones de personas, por lo que nuestros números podrían sobreponderar la infección considerando que se ha visto un virus que demora entre 10-15 días en incubarse.
