El Modelo SIR
Modelo SIR con demografía estable (tasa natalidad y mortalidad iguales) y alta inmunidad

El Modelo SIR

Un modelo clásico de epidemiología es el llamado modelo SIR de Kermack y McKendrick (Suceptibles-Infectados_Recuperados):

Modelo SIR básico

Sistemas de Ecuaciones Diferenciales

La finalidad de los modelos mecanísticos es poder determinar el estado de los compartimentos en algún momento del tiempo.
Para esto es necesario conocer la tasa de flujo de los individuos desde un compartimento hacia otro.
Las ecuaciones diferenciales son una herramienta matemática para plantear el problema.
Diversos procedimientos computacionales permiten la resolución de sistemas de ecuaciones diferenciales (cuando no hay una solución analítica conocida), y la obtención de las funciones para determinar el estado de los compartimentos en cualquier momento.
Algunos de estos procedimientos requieren conocer los parámetros de las ecuaciones (desde la teoría o empíricamente), y las condiciones iniciales del sistema.

Modelo SIR con tasas de flujo

Uso de deSolve en R

El paquete deSolve en R permite resolver sistemas de ecuaciones diferenciales ordinarias (ODE), asi como otras, cuando se conocen las condiciones iniciales (IVP).
Las ecuaciones básicas del modelo SIR (sin demografía, o nacimientos iguales a muertes), para cada compartimento, son las siguientes:

\[\frac {dS}{dt} = -\beta SI\] \[\frac {dI}{dt} = \beta SI - \gamma I\] \[\frac {dR}{dt} = \gamma I\]

Código R

library(deSolve)
#tamaño poblacional
N = 1
#estado inicial de los compartimentos
init <- c(S = 1-1e-6,
          I = 1e-6,
          R = 0)
#parámetros del modelo (coeficientes de las variables)
param <- c(beta = 1.4247,
           gamma = 0.14286)
#crear la función con las ODE
sir <- function(times, init, param) {
  with(as.list(c(init, param)), {
#ecuaciones diferenciales   
    dS <- -beta * S * I
    dI <-  beta * S * I - gamma * I
    dR <-                 gamma * I
#resultados de las tasas de cambio    
    return(list(c(dS, dI, dR)))
  })
}
#intervalo de tiempo y resolución
times <- seq(0, 70, by = 1)
#resolver el sistema de ecuaciones con función 'ode'
out <- ode(y = init, times = times, func = sir, parms = param)
#cambiar out a un data.frame
out <- as.data.frame(out*N) #aqui puede multiplicar 'out' por N
#eliminar la variable 'time' en out
out$time <- NULL
#mostrar 10 primeros datos
head(out, 10)

##            S            I            R
## 1  0.9999990 1.000000e-06 0.000000e+00
## 2  0.9999957 3.968277e-06 3.308132e-07
## 3  0.9999842 1.433888e-05 1.486622e-06
## 4  0.9999414 5.280414e-05 5.773735e-06
## 5  0.9997859 1.926859e-04 2.136590e-05
## 6  0.9992231 6.990818e-04 7.783497e-05
## 7  0.9972028 2.516429e-03 2.807825e-04
## 8  0.9899906 9.000750e-03 1.008659e-03
## 9  0.9648937 3.152281e-02 3.583502e-03
## 10 0.8846221 1.030844e-01 1.229344e-02

#gráfica
matplot(x = times, y = out, type = "l",
        xlab = "Tiempo", ylab = "S, I, R", main = "Modelo SIR básico",
        lwd = 1, lty = 1, bty = "l", col = 2:4)
#añadir leyenda de líneas
legend(40, 0.7, c("Susceptibles", "Infectados", "Recuperados"), 
       pch = 1, col = 2:4, bty = "n", cex = 1)

Uso de scipy.integrate en Python

El ‘ecosistema’ SciPy, basado en Python, contiene paquetes que permiten resolver sistemas de ecuaciones, integrarlas y graficar los resultados (incluye Matplotlib).
Usamos la función scipy.integrate.odeint para resolver (e integrar) el sistema de ecuaciones diferenciales del modelo SIR.

Código Python

## 3.6.5 |Anaconda, Inc.| (default, Apr 26 2018, 08:42:37) 
## [GCC 4.2.1 Compatible Clang 4.0.1 (tags/RELEASE_401/final)]

import scipy.integrate as spi
import numpy as np
import pylab as pl
'''tamaño poblacional'''
N=1
beta=1.4247
gamma=0.14286
'''time step'''
TS=1.0 
ND=70.0
S0=1-1e-6
I0=1e-6
INPUT = (S0, I0, 0.0)


def diff_eqs(INP,t):  
    Y=np.zeros((3))
    V = INP
    '''Las ecuaciones diferenciales'''
    Y[0] = - beta * V[0] * V[1]
    Y[1] = beta * V[0] * V[1] - gamma * V[1]
    Y[2] = gamma * V[1]
    return Y   # For odeint

t_start = 0.0; t_end = ND; t_inc = TS
t_range = np.arange(t_start, t_end+t_inc, t_inc)
RES = spi.odeint(diff_eqs,INPUT,t_range)

#Gráfica
pl.plot(RES[:,0]*N, '-g', label='Susceptibles')
pl.plot(RES[:,2]*N, '-k', label='Recuperados')
pl.plot(RES[:,1]*N, '-r', label='Infectados')
pl.legend(loc=0)
pl.title('Modelo SIR básico')
pl.xlabel('Tiempo')
pl.savefig('sirpy')

Gráfica Modelo SIR

Análisis del modelo SIR básico

Aunque el modelo SIR básico, sin demografía, es bastante simple en su estructura, sin embargo es muy útil para resaltar algunas propiedades de los sistemas epidemiológicos. A continuación exploramos algunas de esas propiedades, cambiando los parámetros del modelo SIR básico.

El umbral de la infección: \(\gamma / \beta\)

Examinaremos primero que efecto tiene disminuir el tamaño de la población inicial de susceptibles (S0), sin modificar los otros parámetros del modelo (\(\beta,\gamma\)). Para esto vamos a disminuir S0 de manera exponencial, y calcular el correspondiente valor máximo de I.

Ejemplo de cálculo de Imax

#usar el archivo de output (out) de cada corrida con los diferentes S0
head(out)

##           S            I            R
## 1 0.9999990 1.000000e-06 0.000000e+00
## 2 0.9999957 3.968277e-06 3.308132e-07
## 3 0.9999842 1.433888e-05 1.486622e-06
## 4 0.9999414 5.280414e-05 5.773735e-06
## 5 0.9997859 1.926859e-04 2.136590e-05
## 6 0.9992231 6.990818e-04 7.783497e-05

#obtener el máximo de I de cada corrida
imax <- max(out$I)
imax

## [1] 0.6603265

Gráfica de Imax vs logSo para determinar umbral de infección

#resultados de cálculos realizados con SIR básico
So <- c(1-1e-6, 1-1e-5, 1-1e-4, 1-1e-3, 1-1e-2, 1-1e-1, 
        0.8, 0.6, 0.4, 0.2, 0.1, 0.05)
Imax <- c(0.6675381,0.6675321,0.6674719,0.6668648,0.660271,
          0.5780076,0.4912547,0.3207578,0.1604105,0.0306711,
          1e-06, 1e-06)
#logaritmo de S0 para mejor visualización
logSo <- log(So)
#gráfica del logaritmo de S0 vs I máximo
plot(logSo,Imax)

Por debajo de cierto valor de S0 la infección no procede. Podemos obtener este valor de forma analítica rearreglando la ecuación del cambio en infectados: \[\frac{dI}{dt} = I(\beta S - \gamma)\] Si igualamos el paréntesis a 0: \[\beta S - \gamma = 0\] tenemos que para \(S \leq \frac{\gamma}{\beta}\), el cambio en los infectados es 0 o negativo, y la infección no procede. El valor de umbral de S0, \(\gamma / \beta\), es la tasa relativa de remoción de infectados, y su valor debe ser suficientemente pequeño para que se desarrolle la epidemia.

La vacunación tiene el efecto equivalente de disminuir S0, y de esta manera disminuir (y retrasar el pico de) los infectados. Sin embargo, los resultados indican que se necesita una alta proporción de vacunados (bajo S0) para controlar la epidemia.

Período infeccioso (\(1/\gamma\)) y \(R_0\)

Al inverso de la tasa relativa de remoción, se le conoce como la razón básica de reproducción o \(R_0\), y es uno de los valores más usados en epidemiología, y se define como: “the average number of secondary cases arising from an average primary case in an entirely susceptible population” (Keeling and Rohani, 2008). Con la misma argumentación anterior, en una población donde todos sean susceptibles (S0 = 1), para que un patógeno pueda invadir es necesario que \(R_0 > 1\).

El parámetro \(\gamma\) es conocido como la tasa de remoción o recuperación, aunque a menudo el interés es principalmente en su inverso, \(1/\gamma\), que determina el promedio del período infeccioso.

El período infeccioso de una enfermedad usualmente se puede determinar con datos epidemiológicos, y \(R_0\) puede calcularse contabilizando los casos nuevos producidos por un individuo infectado, y multiplicando por el período infeccioso (en una población con todos susceptibles inicialmente).

A continuación tenemos un ejemplo de valores de \(R_0\) y período infeccioso para algunas enfermedades infecciosas.

Tabla enfermedades infecciosas Ro

(tomada de Keeling and Rohani, 2008)

Tabla periodo infeccioso enfermedades

(tomada de Wearing y col., )

Modelo SIR con demografía estable (tasa natalidad y mortalidad iguales) y alta inmunidad

El modelo SIR básico anterior permite analizar el caso de enfermedades infecciosas que producen una epidemia relativamente rápido y en una población mayormente susceptible (“outbreak”). En este modelo el agente infeccioso eventualmente desaparece, por la falta de una población susceptible para su mantenimiento; esto significa que no se establece una enfermedad infecciosa endémica.

Para modelar a más largo plazo, debemos incorporar cambios demográficos, que incluyan la entrada de nuevos individuos susceptibles por medio de la natalidad; al mismo tiempo, para un modelo realista, debemos incluir la mortalidad en el modelo. Vamos a examinar ahora un modelo que incluye un parámetro demográfico (\(\mu\)), igual para natalidad y mortalidad (población estable), lo cual es muy cercano a lo que ocurre en países con cierta estabilidad económica y social.

El parámetro \(\mu\) se estima a partir de una esperanza de vida del hospedero igual a \(1/\mu\). Si en una población la esperanza de vida es 70 años (\(1/\mu\)), entonces:

\(\mu = \frac{1}{70 * 365}\) (en días)

Esta tasa de mortalidad y natalidad es independiente del patógeno, y se aplica por igual a todos los compartimentos del modelo SIR.

Las nuevas ecuaciones diferenciales son las siguientes:

\[\frac{dS}{dt} = \mu - \beta SI -\mu S \] \[\frac{dI}{dt} = \beta SI - \gamma I - \mu I\] \[\frac{dR}{dt} = \gamma I - \mu R\]

Modelo SIR demográfico en Python

import scipy.integrate as spi
import numpy as np
import pylab as pl
#parámetros
mu=1/(70*365.0)
beta=520/365.0
gamma=1/7.0
TS=1.0
ND=70*365
#condiciones iniciales
S0=0.1
I0=1e-4
R0=1-S0-I0
INPUT = (S0, I0, R0)
#ecuaciones diferenciales
def diff_eqs(INP,t):  
    '''The main set of equations'''
    Y=np.zeros((3))
    V = INP    
    Y[0] = mu - beta * V[0] * V[1] - mu * V[0]
    Y[1] = beta * V[0] * V[1] - gamma * V[1] - mu * V[1]
    Y[2] = gamma * V[1] - mu * V[2]
    return Y   # For odeint

t_start = 0.0; t_end = ND; t_inc = TS
t_range = np.arange(t_start, t_end+t_inc, t_inc)
RES = spi.odeint(diff_eqs,INPUT,t_range)

#Ploting
pl.subplot(311)
pl.plot(RES[:,0], '-g', label='Susceptibles')
pl.title('SIR demográfico')
pl.xlabel('Tiempo (días)')
pl.ylabel('Susceptibles')
pl.subplot(312)
pl.plot(RES[:,1], '-r', label='Infecciosos')
pl.xlabel('Tiempo (días)')
pl.ylabel('Infecciosos')
pl.subplot(313)
pl.plot(RES[:,2], '-k', label='Recuperados')
pl.xlabel('Tiempo (días)')
pl.ylabel('Recuperados')
pl.savefig('sirdempy')

Gráfica SIR demográfico

Algunas observaciones a partir del modelo demográfico

\(R_0\) se define como el número promedio de nuevos casos que provienen de un caso promedio de una población totalmente susceptible, y se puede calcular como \(R_0 = \frac{\beta}{\gamma}\), lo que es equivalente a decir que es la multiplicación de la tasa de contagio por contacto individual por el tiempo promedio infeccioso.

REFERENCIAS

Beckley, R., Weatherspoon, C., Alexander, M., Chandler, M., Johnson, A., Bhatt, G.S. n.d. Modeling epidemics with diﬀerential equations. Tennessee State University
Keeling, M.J., Rohani, P. 2008. Modeling infectious diseases in humans and animals. Princeton University Press, Princeton.
Krylova, O., Earn, D.J.D. 2013. Effects of the infectious period distribution on predicted transitions in childhood disease dynamics. J R Soc Interface 10. https://doi.org/10.1098/rsif.2013.0098
Wearing, H.J., Rohani, P., Keeling, M.J. 2005. Appropriate Models for the Management of Infectious Diseases. PLoS Med 2. https://doi.org/10.1371/journal.pmed.0020174

Modelos Mecanísticos Determinísticos con Ecuaciones Diferenciales

D. S. Fernández del Viso

24 Septiembre 2018