library(deSolve) # Paquete deSolve para resolver las ecuaciones diferenciales
library(tidyverse) # Paquetes ggplot2 y dplyr de tidyverse
library(cowplot) # Paquete gridExtra para unir gráficos.
#> Warning: package 'cowplot' was built under R version 4.4.1
#>
#> Attaching package: 'cowplot'
#> The following object is masked from 'package:lubridate':
#>
#> stamp
# desafio 1 ---------------------------------------------------------------
# Parámetros
ph <- 0.7 # Probabilidad de transmisión del vector al hospedador dada una picadura por un mosquito infeccioso a un humano susceptible
pv <- 0.7 # Probabilidad de transmisión del hospedador al vector dada una picadura por un mosquito susceptible a un humano infeccioso
Lv <- 10 # Esperanza de vida de los mosquitos (en días)
Lh <- 50 * 365 # Esperanza de vida de los humanos (en días)
Iph <- 7 # Periodo infeccioso en humanos (en días)
IP <- 6 # Periodo infeccioso en vectores (en días)
EIP <- 8.4 # Período extrínseco de incubación en mosquitos adultos
muv <- 1/Lv # Tasa per capita de mortalidad del vector (1/Lv)
muh <- 1/Lh # Tasa per capita de mortalidad del hospedador (1/Lh)
alphav <- muv # Tasa per capita de natalidad del vector
alphah <- muh # Tasa per capita de natalidad del hospedador
gamma <- 1/Iph # Tasa de recuperación en humanos
delta <- 1/EIP # Tasa extrínseca de incubación
# Tamaño de la población
Nh <- 100000 # Número de humanos
m <- 2 # Proporción vector a humano
Nv <- m * Nh # Número de vectores
R0 <- 3 # Número reproductivo
b <- sqrt((R0 * muv*(muv+delta) * (muh+gamma)) /
(m * ph * pv * delta)) # tasa de picadura
betah <- ph*b # Coeficiente de transmisión del mosquito al humano
betav <- pv*b # Coeficiente de transmisión del humano al mosquito
TIME <- 100 # Número de años a simular
# modelo en R -------------------------------------------------------------
# ode()
zika_model <- function(tiempo, variable_estado, parametros) {
with(as.list(c(variable_estado, parametros)), # entorno local para evaluar derivados
{
# Humanos
dSh <- (alphah * Nh) - (betah * (Iv/Nh) * Sh) - (muh * Sh)
dIh <- (betah * (Iv/Nh) * Sh) - ((gamma + muh) * Ih)
dRh <- (gamma * Ih) - (muh * Rh)
# Vectores
dSv <- (alphav * Nv) - (betav * (Ih/Nh) * Sv) - (muv * Sv)
dEv <- (betav * (Ih/Nh) * Sv) - ((delta + muv)* Ev)
dIv <- (delta * Ev) - (muv * Iv)
list(c(dSh, dIh, dRh, dSv, dEv, dIv))
}
)
}
# ode(func = zika_model)
# desafio 3 ---------------------------------------------------------------
#Tiempo
tiempo_simulacion <- seq(1, 365 * TIME , by = 1)
# Especifique los parámetros
parametros_simulacion <- c(
muv = muv,
muh = muh,
alphav = alphav,
alphah = alphah,
gamma = gamma,
delta = delta,
betav = betav,
betah = betah,
Nh = Nh,
Nv = Nv
)
# Condiciones iniciales del sistema
condiciones_inicio <- c(
Sh = Nh , # Número inicial de Sh en T0, total de la poblacion
Ih = 0, # Número inicial de Ih en T0, ningun infectado
Rh = 0, # Número inicial de Rh en T0, ningun recuperado
Sv = Nv-1, # Número inicial de Sv en T0, habra un mosquito infectado
Ev = 0, # Número inicial de Ev en T0, ningun expuesto
Iv = 1 # Número inicial de Iv en TO, unico mosquito infectado
)
# Resuelva las ecuaciones
out <- ode(
y = condiciones_inicio, # Condiciones iniciales
times = tiempo_simulacion, # Tiempo
func = zika_model, # Modelo
parms = parametros_simulacion # Parámetros
) %>%
as.data.frame() # Convertir a data frame
# out
out %>% as_tibble()
#> # A tibble: 36,500 × 7
#> time Sh Ih Rh Sv Ev Iv
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 100000 0 0 199999 0 1
#> 2 2 100000. 0.176 0.0131 199999. 0.0338 0.906
#> 3 3 100000. 0.313 0.0484 199999. 0.116 0.828
#> 4 4 99999. 0.418 0.101 199999. 0.225 0.769
#> 5 5 99999. 0.501 0.167 199999. 0.346 0.728
#> 6 6 99999. 0.566 0.243 199999. 0.469 0.705
#> 7 7 99999. 0.620 0.328 199999. 0.589 0.698
#> 8 8 99999. 0.667 0.420 199999. 0.703 0.705
#> 9 9 99999. 0.711 0.518 199998. 0.811 0.724
#> 10 10 99999. 0.753 0.623 199998. 0.913 0.753
#> # ℹ 36,490 more rows
# desafio 4 ---------------------------------------------------------------
out$years <- out$time / 365
out$weeks <- out$time / 7
out %>% as_tibble()
#> # A tibble: 36,500 × 9
#> time Sh Ih Rh Sv Ev Iv years weeks
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 100000 0 0 199999 0 1 0.00274 0.143
#> 2 2 100000. 0.176 0.0131 199999. 0.0338 0.906 0.00548 0.286
#> 3 3 100000. 0.313 0.0484 199999. 0.116 0.828 0.00822 0.429
#> 4 4 99999. 0.418 0.101 199999. 0.225 0.769 0.0110 0.571
#> 5 5 99999. 0.501 0.167 199999. 0.346 0.728 0.0137 0.714
#> 6 6 99999. 0.566 0.243 199999. 0.469 0.705 0.0164 0.857
#> 7 7 99999. 0.620 0.328 199999. 0.589 0.698 0.0192 1
#> 8 8 99999. 0.667 0.420 199999. 0.703 0.705 0.0219 1.14
#> 9 9 99999. 0.711 0.518 199998. 0.811 0.724 0.0247 1.29
#> 10 10 99999. 0.753 0.623 199998. 0.913 0.753 0.0274 1.43
#> # ℹ 36,490 more rows
out %>% glimpse()
#> Rows: 36,500
#> Columns: 9
#> $ time <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 1…
#> $ Sh <dbl> 100000.00, 99999.81, 99999.64, 99999.48, 99999.33, 99999.19, 999…
#> $ Ih <dbl> 0.0000000, 0.1759705, 0.3125805, 0.4182923, 0.5006964, 0.5662925…
#> $ Rh <dbl> 0.00000000, 0.01309050, 0.04840205, 0.10092344, 0.16679457, 0.24…
#> $ Sv <dbl> 199999.0, 199999.1, 199999.1, 199999.0, 199998.9, 199998.8, 1999…
#> $ Ev <dbl> 0.00000000, 0.03383092, 0.11623920, 0.22521428, 0.34594459, 0.46…
#> $ Iv <dbl> 1.0000000, 0.9061977, 0.8282000, 0.7686612, 0.7279212, 0.7049395…
#> $ years <dbl> 0.002739726, 0.005479452, 0.008219178, 0.010958904, 0.013698630,…
#> $ weeks <dbl> 0.1428571, 0.2857143, 0.4285714, 0.5714286, 0.7142857, 0.8571429…
# parte 11 ----------------------------------------------------------------
# humanos -----------------------------------------------------------------
# Revise el comportamiento general del modelo para 100 años
p1h <- ggplot(data = out, aes(y = (Rh + Ih + Sh), x = years)) +
geom_line(color = 'grey68', size = 1) +
ggtitle('Población humana total') +
theme_bw() + ylab('Número') + xlab('Años')
#> Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
#> ℹ Please use `linewidth` instead.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.
p2h <- ggplot(data = out, aes(y = Sh, x = years)) +
geom_line(color = 'royalblue', size = 1) +
ggtitle('Población humana susceptible') +
theme_bw() + ylab('Número') + xlab('Años')
p3h <- ggplot(data = out, aes(y = Ih, x = years)) +
geom_line(color = 'firebrick', size = 1) +
ggtitle('Población humana infecciosa') +
theme_bw() + ylab('Número') + xlab('Años')
p4h <- ggplot(data = out, aes(y = Rh, x = years)) +
geom_line(color = 'olivedrab', size = 1) +
ggtitle('Población humana recuperada') +
theme_bw() + ylab('Número') + xlab('Años')
plot_grid(p1h, p2h, p3h, p4h, ncol = 2)
# vectores ----------------------------------------------------------------
# Revise el comportamiento general del modelo
p1v <- ggplot(data = out, aes(y = (Sv + Ev + Iv), x = years)) +
geom_line(color = 'grey68', size = 1) +
ggtitle('Población total de mosquitos') +
theme_bw() + ylab('Número') + xlab('Años')
p2v <- ggplot(data = out, aes(y = Sv, x = years)) +
geom_line(color = 'royalblue', size = 1) +
ggtitle('Población susceptible de mosquitos') +
theme_bw() + ylab('Número') + xlab('Años')
p3v <- ggplot(data = out, aes(y = Ev, x = years)) +
geom_line(color = 'orchid', size = 1) +
ggtitle('Población expuesta de mosquitos') +
theme_bw() + ylab('Número') + xlab('Años')
p4v <- ggplot(data = out, aes(y = Iv, x = years)) +
geom_line(color = 'firebrick', size = 1) +
ggtitle('Población infecciosa de mosquitos') +
theme_bw() + ylab('Número') + xlab('Años')
plot_grid(p1v, p2v, p3v, p4v, ncol = 2)
# proporcion --------------------------------------------------------------
p1 <- ggplot(data = out, aes(y = Sh/(Sh+Ih+Rh), x = years)) +
geom_line(color = 'royalblue', size = 1) +
ggtitle('Población humana susceptible') +
theme_bw() + ylab('Proporción') + xlab('Años') +
coord_cartesian(ylim = c(0,1))
p2 <- ggplot(data = out, aes(y = Ih/(Sh+Ih+Rh), x = years)) +
geom_line(color = 'firebrick', size = 1) +
ggtitle('Población humana infecciosa') +
theme_bw() + ylab('Proporción') + xlab('Años') +
coord_cartesian(ylim = c(0,1))
p3 <- ggplot(data = out, aes(y = Sv/(Sv + Ev + Iv), x = years)) +
geom_line(color = 'olivedrab', size = 1) +
ggtitle('Población mosquitos susceptible') +
theme_bw() + ylab('Proporción') + xlab('Años') +
coord_cartesian(ylim = c(0,1))
p4 <- ggplot(data = out, aes(y = Iv/(Sv + Ev + Iv), x = years)) +
geom_line(color = 'darkgreen', size = 1) +
ggtitle('Población mosquitos infectada') +
theme_bw() + ylab('Proporción') + xlab('Años') +
coord_cartesian(ylim = c(0,1))
plot_grid(p1, p2, p3, p4, ncol = 2) 
# la primera epidemia -----------------------------------------------------
# Revise la primera epidemia
dat <- out %>% as_tibble() %>% filter(weeks < 54)
p1e <- ggplot(dat, aes(y = Ih, x = weeks)) +
geom_line(color = 'firebrick', size = 1) +
ggtitle('Población de humanos infecciosos') +
theme_bw() + ylab('Número') + xlab('Semanas')
p2e <- ggplot(dat, aes(y = Rh, x = weeks)) +
geom_line(color = 'olivedrab', size = 1) +
ggtitle('Población humana recuperada') +
theme_bw() + ylab('Número') + xlab('Semanas')
p3e <- ggplot(dat, aes(y = Sv, x = weeks)) +
geom_line(color = 'orchid', size = 1) +
ggtitle('Población de mosquitos susceptibles') +
theme_bw() + ylab('Número') + xlab('Semanas')
p4e <- ggplot(dat, aes(y = Iv, x = weeks)) +
geom_line(color = 'gray68', size = 1) +
ggtitle('Población mosquitos infectados') +
theme_bw() + ylab('Número') + xlab('Semanas')
plot_grid(p1e, p2e)
plot_grid(p1e, p2e,p4e, p3e)
Created on 2024-08-02 with reprex v2.1.1