Dado que hoy en día estamos en el medio de una pandemia me pareció adecuado intentar entender el ABC de la epidemiología matemática. Quizá en otro universo nosotros los estudiantes de matemática estaríamos trabajando en estos modelos respondiendo a la necesidad de la sociedad, pero no es así. En esta entrada quería intentar entender el modelo más básico para modelizar la cantidad de infectados por una pandemia en función del tiempo. Para eso voy a tener que desempolvar mi (escaso) conocimiento de ecuaciones diferenciales y un poco de código para poder experimentar tranquilamente.

Las ecuaciones del modelo.

Nuestro objetivo es intentar modelizar la cantidad de infectados qué va a haber por día sabiendo los que hay hoy y usando lo que se conoce del COVID-19. ¿Por qué nos interesa esto? Porque nuestro frágil sistema de salud tiene una capacidad muy limitada de atender gente. Es decir que si tenemos más de N para N chico infectados en un día no va a haber lugar donde tratarlos y la situación rápidamente podría pasar a ser como en Lombardía o Guayaquil (obviando que hay otras causas para las situaciones en esos lugares). Ahora que sabemos porqué esto es malo pasemos a la instancia de modelizar.

Nuestro modelo de lo más naif es el SIR (Susceptible Infectado Recuperado) que particiona a la población en esas tres clases. Una suposición que estamos haciendo viene dada de que esta pandemia es causada por un virus por lo tanto la gente recuperada no se vuelve a contagiar en el rango de tiempo que nosotros modelizamos (esto no quiere decir que la inmunidad dure para siempre pero para lo rápido que se esparce el virus podemos pensarlo que este es el caso). A su vez como la epidemia se esparce tan rápido no hace falta considerar los nacimientos y muertes de la población así que la podemos asumir constante $N$. Las ecuaciones diferenciales nos relacionan estas 3 particiones de la población de la manera esperada. La primer ecuación es la siguiente

\[\dfrac{dS}{dt} = - \dfrac{\beta I S}{N}\]

La derivada es negativa porque la cantidad de susceptibles naturalmente va bajando porque se transforman en infectados. La constante $\beta$ depende del virus y de las medidas que se toman para evitar el contagio (sea una cuarentena obligatoria o aislamiento de infectados, etc).
La segunda ecuación nos dice el cambio en los infectados,

\[\dfrac{dI}{dt} = \dfrac{\beta I S}{N} - \gamma I\]

donde aparecen todos los susceptibles que se infectan y le restamos los infectados que se recuperan. Esta segunda constante que aparece $\gamma$ tiene que ver con el periódo de incubación de la enfermedad y está en unidades $1/t$ donde $t$ es la unidad de tiempo, que generalmente suele ser días. En el caso del COVID-19 se necesitan aproximadamente 14 días para saber si alguien se va a recuperar de la enfermedad (puede llegar a fallecer en el camino a cumplir esta espera). Lo importante es que durante esos 14 días puede seguir contagiando otras personas. Finalmente los recuparados siguen esta ecuación,

\[\dfrac{dR}{dt} = \gamma I\]

que tiene que ver con la cantidad de infectados que se recuperan. Notemos que no aparece en ningún momento la mortalidad esto se debe a qué para nuestro modelo no nos estamos preocupando por eso (se estima que la mortalidad es un >4% para los infectados que son atendidos) sino que nos estamos preocupando por la cantidad de gente que se puede llegar a infectar al mismo tiempo para evitar el colapsamiento del sistema de salud.

Transición entre los diferentes compartimientos de la población en una epidemia.

El número de reproducción.

Si consideramos las constantes que antes aparecieron en el modelo SIR podemos armarnos el siguiente número, $R_0 = \beta/\gamma$ llamado número de reproducción. Moralmente este número nos dice la cantidad de gente que un infectado infecta en promedio (decimos en promedio porque existen casos como los super-spreaders). Este número es fundamental(!) ya que decide,

  1. Si hay una epidemia o no la hay. Si $R_0$ es mayor que 1 entonces cada vez va a haber más infectados y vamos a llegar a tener una epidemia. En cambio si es menor resulta que el virus va a ir siendo cada vez menos común hasta desaparecer.
  2. Determina la cantidad de infectados que va a haber para el final de la epidemia y de esta manera las casualidades del virus. A mayor $R_0$ mayor cantidad de infectados
  3. Determina la taza de crecimiento de la epidemia. Esto sale de las ecuaciones, que no solo a mayor valor de $R_0$ van a llegar a mayor cantidad de infectados sino que va a llegar más rápido (!).

Entonces la pregunta natural es cómo manejamos que este número sea lo más chico posible. Para eso notemos que esta constante $R_0$ depende de $\beta$ que tiene que ver con la capacidad de infección de la enfermedad en las circustancias dadas (en caso de cuarentena obligatoria para una misma enfermedad es menor que en caso de no estar haciendo absolutamente nada para frenarla). Mientras que $\gamma$ está fija y depende intrínsecamente del virus, por lo tanto no queda más que modificar $\beta$.

Herd immunity.

Otro temita que se está escuchando mucho sobretodo para los que siguen al Reino Unido es el tema de la herd immunity. Básicamente es una estrategia que querían tomar para que la población quede expuesta al virus y desarollo la inmunidad. Con el tiempo la epidemia debería dejar de existir debido a que la mayoría sería inmune a ella. La estrategia resulto ser un fracaso en parte porque mucha más gente necesita hospitalización de lo que era prevista y se pasó a una cuarentena flexible dentro de todo.

La pregunta que está de fondo a esta política es la siguiente. ¿Cuánta gente se necesita inmunizar, ya sea por una vacuna que por el momento es inexistente o por contagiarse y recuperarse, para lograr calmar la epidemia? Consideremos que un porcentaje $0 < p < 1$ de la población tiene inmunidad. Veamos que tan grande tiene que ser $p$ para evitar la epidemia o en otras palabras, como ya sabemos, que $R_0 \le 1.$

Notemos que $R_0$ como nombre de una constante es bastante curioso por el $0$ de sub-índice. Este nos señala que en realidad es una función del tiempo y la estamos mirando en el momenot $t=0$ que es donde estamos parados en el modelo. ¿Qué pasaría si el 90% de la población está infectada? En ese caso es claro que la cantidad de gente que un infectado logra infectar va a ser menor que si el 10% lo estuviese. Esto nos dice en cierta manera que el $R$ depende de la cantidad de gente que está infectada en el momento. Para eso notemos que si $s =1-p$ es el porcentaje de gente infectada en un momento dado entonces \(R = sR_0\) y recordemos que nos gustaría que $R < 1$ entonces despejando de esta desigualdad obtenemos el siguiente valor de $p$ mínimo que garantiza que no haya una epidemia, \(p= 1 - \dfrac{1}{R_0}.\) Esto nos dice que a su vez, $R_0$ nos define qué porcentaje de la población es necesario inmunizar para garantizar que se caiga la epidemia. En el caso que nos compete que es el del covid se estima que $R_0$ es aproximadamente 2,4 acá en Argentina por lo que $p>0.65$ para tener la codiciada herd immunity.

Código.

Acá pongo un código simpático en python que encontré en interné para poder ver distintos escenarios en nuestra máquina.

#Población total, N.
N = 500000
# Infectados iniciales I0 y recuperados iniciales R0.
I0, R0 = 10, 0
# El resto, S0 es susceptible.
S0 = N - I0 - R0
### Las constantes del modelo
beta, gamma = 0.18, 1./14 
### Cuanto tiempo se extiende el modelo (en días)
t = np.linspace(0, 365, 365)

### Las ecuaciones diferenciales del SIR.
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 de condiciones iniciales.
y0 = S0, I0, R0
### Resolvedor (mágico) de las ecuaciones diferenciales.
ret = odeint(deriv, y0, t, args=(N, beta, gamma))
S, I, R = ret.T

### Ploteo de tres curvas S(t), I(t), R(t).
fig = plt.figure(facecolor='w')
ax = fig.add_subplot(111, axis_bgcolor='#dddddd', axisbelow=True)
ax.plot(t, S/1000, 'b', alpha=0.5, lw=2, label='Susceptible')
ax.plot(t, I/1000, 'r', alpha=0.5, lw=2, label='Infected')
ax.plot(t, R/1000, 'g', alpha=0.5, lw=2, label='Recovered with immunity')
ax.set_xlabel('Time /days')
ax.set_ylabel('Number (1000s)')
ax.set_ylim(0,5000)
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.show()

Conclusión.

Este modelo si bien es simpático y fácil de entender dentro de todo, tiene algunas desventajas como asumir que $\beta$ es realmente constante o que el virus se propaga de forma uniforme en la población pero al ser suficientemente sencillo y general encuentra su uso en varios casos en particular. Algunas mejores de este sistema son las que se usan para predecir el día que va a haber un pico de infectados en distintos países como EEUU, Argentina, Alemania (no encuentro la cita por el momento pero la había visto en twitter). El modelo ya tiene casi 100 años y se nota por la simplicidad que nos remite al modelo de Lotka-Volterra de una población de depredadores y sus presas.

Referencia.

  • Emilia Vynnycky, Richard White - An Introduction to Infectious Disease Modelling - OUP Oxford (2010)