Modelo SIR con Ecuaciones Diferenciales
D. S. Fernandez del Viso
10/2/2019
Ecuaciones diferenciales para el modelo SIR
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\]
Descripción y parámetros del modelo
- El modelo SIR asume el contacto (directo o mediante vector) entre individuos susceptibles (X) e infectados (Y), lo cual significa una transmisión por acción de masas, y por lo tanto es dependiente de la frecuencia de cada grupo actuando (susceptibles e infectados). Por esta razón, y para estandarizar el modelo para cualquier tamaño poblacional (N), utilizamos las proporciones de cada grupo:
\[S = X/N\] \[I = Y/N\] Igual lo hacemos para los recuperados en el modelo (R): \[R = Z/N\]
- El parámetro \(\beta\) es el producto de la tasa de contacto y la probabilidad de transmisión de la infección. También se utiliza la expresión \(\lambda = \beta I\) para indicar la fuerza de infección.
- El parámetro \(\gamma\) se denomina tasa de remoción o recuperación (usaremos recuperación); a partir de datos epidemiológicos se puede obtener el perÃodo infeccioso, que es el recÃproco \(1/\gamma\).
- Un parámetro derivado de los anteriores, es uno de los números más famosos de la epidemiologÃa es la razón reproductiva básica, \(R_0\), que se cuantifica como el número promedio de casos secundarios que se derivan de un caso primario promedio, en una población totalmente susceptible.
- Podemos demostrar que cuando la proporción de susceptibles iniciales, \(S(0)\), es menor que \(\gamma /\beta\), la tasa relativa de remoción, la infección no procede (esto se conoce como el umbral de infección).
- El inverso de \(\gamma /\beta\), por otra parte, expresa la razón reproductiva básica, \(R_0\).
Ejercicios con el modelo SIR y ecuaciones diferenciales en R
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.
Resolución de ecuaciones del modelo SIR
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 = 2.0,
gamma = 0.2)
#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, 20, 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
# ver resultados
head(out)## time S I R
## 1 0 0.9999990 1.000000e-06 0.000000e+00
## 2 1 0.9999915 7.786995e-06 7.541157e-07
## 3 2 0.9999469 4.793343e-05 5.214994e-06
## 4 3 0.9996718 2.955068e-04 3.272917e-05
## 5 4 0.9980070 1.793592e-03 1.993999e-04
## 6 5 0.9880578 1.074087e-02 1.201331e-03A continuación vamos a rearreglar los datos de los resultados (out), usando el paquete tidyverse. Esto facilita otros análisis y la gráficación de los resultados.
library("tidyverse")
new.out <- as.data.frame(out) %>% gather(key, value, -time)
head(new.out)## time key value
## 1 0 S 0.9999990
## 2 1 S 0.9999915
## 3 2 S 0.9999469
## 4 3 S 0.9996718
## 5 4 S 0.9980070
## 6 5 S 0.9880578En este nuevo objeto, new.out, la variable key contiene los grupos \(S\), \(I\), y \(R\), y value el respectivo valor del grupo en ese tiempo (time).
Ahora construiremos una gráfica usando el paquete ggplot2:
library("ggplot2")
ggplot(data=new.out,
aes(x = time,
y = value,
group = key,
col = key
)) +
ylab("Proporción de N") + xlab("Tiempo (dÃas)") +
geom_line(size = 1) +
scale_colour_manual(values = c("red", "green4", "blue"), name = "Grupo") +
scale_y_continuous(labels = waiver(), limits = c(0, 1))
Un dato importante es el valor máximo de infectados:
# valor máximo de I usando un pipe
new.out %>%
filter(key=="I") %>%
filter(value==max(value)) %>%
mutate(maxI = round(N * value, 2)) %>%
select(time, maxI)## time maxI
## 1 9 0.67Modelo SIR para sarampión en Humacao
El problema: si se declara una epidemia de sarampión en Humacao, es muy probable no contar con las camas suficientes para los enfermos. Hay que comprobarlo.
Datos a buscar:
- población de Humacao
- número de camas en hospitales
- parámetros del sarampión
Usar el modelo SIR en R.
Proponer soluciones.
library(deSolve)
#tamaño poblacional
N = 50000
#estado inicial de los compartimentos
init <- c(S = 1-1/50000,
I = 1/50000,
R = 0)
#parámetros del modelo (coeficientes de las variables)
param <- c(beta = 3.4,
gamma = 0.2)
#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, 20, 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)
head(out, 3)## time S I R
## 1 0 0.9999800 0.0000200000 0.0000000000
## 2 1 0.9994588 0.0005105764 0.0000306704
## 3 2 0.9868265 0.0123945561 0.0007789005#rearreglo del output
new.out <- as.data.frame(out) %>%
gather(key, value, -time) %>%
mutate(value=value*N) #de proporción a individuos
head(new.out, 3)## time key value
## 1 0 S 49999.00
## 2 1 S 49972.94
## 3 2 S 49341.33#gráfica conjunta
ggplot(data=new.out,
aes(x = time,
y = value,
group = key,
col = key
)) +
ylab("Individuos") + xlab("Tiempo (dÃas)") +
geom_line(size = 1) +
scale_colour_manual(values = c("red", "green4", "blue"), name = "Grupo") +
scale_y_continuous(labels = waiver(), limits = c(0, 60000))
# nivel máximo de infección
new.out %>%
filter(key=="I") %>%
filter(value==max(value)) %>%
mutate(maxI = round(value, 0)) %>%
select(time, maxI)## time maxI
## 1 4 37648Vacunación
library(deSolve)
#tamaño poblacional
N = 50000
pV = 0.95
effV = 0.95
#estado inicial de los compartimentos
init <- c(S = 1-1/50000-pV*effV,
I = 1/50000,
R = 0,
V = pV*effV)
#parámetros del modelo (coeficientes de las variables)
param <- c(beta = 3.4,
gamma = 0.2)
#crear la función con las ODE
sirv <- 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
dV <- 0
#resultados de las tasas de cambio
return(list(c(dS, dI, dR, dV)))
})
}
#intervalo de tiempo y resolución
times <- seq(0, 90, by = 1)
#resolver el sistema de ecuaciones con función 'ode'
out <- ode(y = init, times = times, func = sirv, parms = param)
#cambiar out a un data.frame
out <- as.data.frame(out)
#eliminar la variable 'V' en out
out$V <- NULL
#mostrar 10 primeros datos
head(out, 3)## time S I R
## 1 0 0.09748000 2.000000e-05 0.000000e+00
## 2 1 0.09747291 2.281319e-05 4.281246e-06
## 3 2 0.09746481 2.602130e-05 9.164461e-06#rearreglo del output
new.out <- as.data.frame(out) %>%
gather(key, value, -time) %>%
mutate(value=value*N) #de proporción a individuos
head(new.out, 3)## time key value
## 1 0 S 4874.000
## 2 1 S 4873.645
## 3 2 S 4873.241#gráfica conjunta
ggplot(data=new.out,
aes(x = time,
y = value,
group = key,
col = key
)) +
ylab("Individuos") + xlab("Tiempo (dÃas)") +
geom_line(size = 1) +
scale_colour_manual(values = c("red", "green4", "blue"), name = "Grupo") +
scale_y_continuous(labels = waiver(), limits = c(0, N-N*pV*effV))
# nivel máximo de infección
new.out %>%
filter(key=="I") %>%
filter(value==max(value)) %>%
mutate(maxI = round(value, 0)) %>%
select(time, maxI)## time maxI
## 1 59 448Modelo SIR con parámetro demográfico \(\mu\)
library(deSolve)
library(tidyverse)
#tamaño poblacional
N = 1
#estado inicial de los compartimentos
init <- c(S = 0.1,
I = 1e-6,
R = 0.9 - 1e-6)
#parámetros del modelo (coeficientes de las variables)
param <- c(beta = 2.4,
gamma = 0.2,
mu = 2e-4)
#crear la función con las ODE
sir <- function(times, init, param) {
with(as.list(c(init, param)), {
#ecuaciones diferenciales
dS <- mu - beta * S * I - mu * S
dI <- beta * S * I - gamma * I - mu * I
dR <- gamma * I - mu * R
#resultados de las tasas de cambio
return(list(c(dS, dI, dR)))
})
}
#intervalo de tiempo y resolución
times <- seq(0, 3650, 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)
#mostrar 3 primeros datos
head(out, 3)## time S I R
## 1 0 0.1000000 1.000000e-06 0.8999990
## 2 1 0.1001797 1.040837e-06 0.8998192
## 3 2 0.1003594 1.083810e-06 0.8996395#gráficas individuales
#infectados
ggplot(data=out,
aes(x = time,
y = I)) +
ylab("Proporción") + xlab("Tiempo (dÃas)") +
ggtitle("Infectados") +
geom_line(size = 0.5, col = "red") +
scale_y_continuous(labels = waiver())
#susceptibles
ggplot(data=out,
aes(x = time,
y = S)) +
ylab("Proporción") + xlab("Tiempo (dÃas)") +
ggtitle("Susceptibles") +
geom_line(size = 0.5, col = "blue") +
scale_y_continuous(labels = waiver())
#recuperados
ggplot(data=out,
aes(x = time,
y = R)) +
ylab("Proporción") + xlab("Tiempo (dÃas)") +
ggtitle("Recuperados") +
geom_line(size = 0.5, col = "green") +
scale_y_continuous(labels = waiver())
Modelo SIR demográfico y aleatorio
library(deSolve)
#tamaño poblacional
N = 1
#estado inicial de los compartimentos
init <- c(S = 0.1,
I = 1e-6,
R = 0.9 - 1e-6)
#parámetros del modelo (coeficientes de las variables)
param <- c(beta = 2.4,
gamma = 0.2,
mu = 2e-4)
#crear la función con las ODE
sir <- function(times, init, param) {
with(as.list(c(init, param)), {
#ecuaciones diferenciales
dS <- mu - beta * S * I * runif(1, 0.999, 1.001) - mu * S
dI <- beta * S * I - gamma * I - mu * I
dR <- gamma * I - mu * R
#resultados de las tasas de cambio
return(list(c(dS, dI, dR)))
})
}
#intervalo de tiempo y resolución
times <- seq(0, 3650, 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
#mostrar 6 primeros datos
head(out)## time S I R
## 1 0 0.1000000 1.000000e-06 0.8999990
## 2 1 0.1001797 1.041861e-06 0.8998192
## 3 2 0.1003594 1.086026e-06 0.8996395
## 4 3 0.1005389 1.132571e-06 0.8994599
## 5 4 0.1007183 1.181755e-06 0.8992802
## 6 5 0.1008976 1.232150e-06 0.8991007#gráficas individuales
#infectados
ggplot(data=out,
aes(x = time,
y = I)) +
ylab("Proporción") + xlab("Tiempo (dÃas)") +
ggtitle("Infectados") +
geom_line(size = 0.5, col = "red") +
scale_y_continuous(labels = waiver())
#susceptibles
ggplot(data=out,
aes(x = time,
y = S)) +
ylab("Proporción") + xlab("Tiempo (dÃas)") +
ggtitle("Susceptibles") +
geom_line(size = 0.5, col = "blue") +
scale_y_continuous(labels = waiver())
#recuperados
ggplot(data=out,
aes(x = time,
y = R)) +
ylab("Proporción") + xlab("Tiempo (dÃas)") +
ggtitle("Recuperados") +
geom_line(size = 0.5, col = "green") +
scale_y_continuous(labels = waiver())