Last active
July 16, 2020 02:01
-
-
Save N8python/2c13b413fe6fe5f27ad027a37a020b2b to your computer and use it in GitHub Desktop.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| library(deSolve) | |
| library(reshape2) | |
| library(ggplot2) | |
| N <- 1000 | |
| initial_state_values <- c( | |
| S = 0.999, | |
| E = 0, | |
| IMI = 0.001, | |
| IM = 0, | |
| A = 0, | |
| H = 0, | |
| C = 0, | |
| ICU = 0, | |
| R = 0, | |
| M = 0) | |
| parameters <- c( | |
| v = 0.0001, # Birth Rate | |
| u = 0.0001, # Death Rate | |
| bmi = 0.429, # Contact Rate for Mild Infections | |
| bm = 0.286, # Contact Rate for Moderate Infections | |
| bh = 0.0715, # Contact Rate for Hospitalized Infections | |
| ba = 0.143, # Contact Rate for Asymptomatic Infections. | |
| t = 0.166, # Proportion of the Exposed Compartment leaving each timestep. | |
| ke = 0.179, # Proportion of those leaving the exposed compartment who become asymptomatic. | |
| ya = 0.143, # Proportion of the Asymptomatic Compartment recovering each timestep. | |
| yimi = 0.143, # Proportion of the Mild Infection Compartment leaving each timestep. | |
| kimi = 0.5, # Proportion of those leaving the Mild Infection Compartment who recover. | |
| yim = 0.143, # Proportion of the Moderate Infection Compartment leaving each timestep. | |
| kim = 0.5, # Proportion of those leaving the Moderate Infection Compartment who recover. | |
| yh = 0.093, # Proportion of the Hospitalized Compartment leaving each timestep. | |
| kh = 0.25, # Proportion of those leaving the Hospitalized Compartment who progress to ARDS/Critical Condition. | |
| kyh = 0.7, # Proportion of those leaving the Hospitalized Compartment who recover. | |
| yc = 0.5, # Proportion of the Critical Condition Compartment leaving each timestep. | |
| yicu = 0.167, # Proportion of the ICU Compartment leaving each timestep. | |
| kicu = 0.4, # Proportion of those leaving the ICU compartment who die. | |
| p = 0.0027, # Proportion of the Recovered Compartment that loses immunity each timestep. | |
| beds = 10 # The Number of Beds in the ICU. | |
| ) | |
| startTime <- 0 # Time when the disease starts - 0 is winter, 180 is summer, 360 is winter, etc. | |
| r0_reduction <- 0.3 # Percent by which the r0 goes down during the summer - 0.3 means a 30% reduction at the peak of summer. | |
| days_simulated <- 7 * 52 * 10 # Set this variable to how many days you want to simulate. (by default, 10 years are simulated.) | |
| times <- seq(from = 0, to = days_simulated, by = 1) # | |
| beta_t <- function(time) { # Function to Calculate Seasonal Variance of r0 | |
| return((r0_reduction / 2)*sin(2*pi/365*(time + 91.25))+(1-r0_reduction / 2)) | |
| } | |
| cohort_model <- function(time, state, parameters) { | |
| with (as.list(c(state, parameters)), { | |
| lambda <- (bmi * IMI + bm * IM + bh * H + ba * A)*beta_t(time) | |
| dS <- v + p * R - lambda * S - u * S | |
| dE <- lambda * S - t * E - u * E | |
| dIMI <- t * (1 - ke) * E - yimi * IMI - u * IMI | |
| dIM <- yimi * (1-kimi) * IMI - yim * IM - u * IM | |
| dA <- t * ke * E - ya * A - u * A | |
| dH <- yim * (1 - kim) * IM - yh * H - u * H | |
| dC <- yh * kh * H - yc * C - u * C | |
| toICU <- min(yc * C, beds / N - ICU) | |
| dICU <- toICU - yicu * ICU - u * ICU | |
| dR <- yimi * kimi * IMI + yim * kim * IM + ya*A + yh * kyh * H + yicu * (1 - kicu) * ICU - p * R - u * R | |
| dM <- yh * (1 - kh - kyh) * H + (yc * C - toICU) + (yicu * kicu * ICU) | |
| return(list(c(dS, dE, dIMI, dIM, dA, dH, dC, dICU, dR, dM))) | |
| }) | |
| } | |
| output <- as.data.frame(ode(y = initial_state_values, times = times, func = cohort_model, parms = parameters)) | |
| output_long <- melt(as.data.frame(output), id = "time") | |
| ggplot(data = output_long, aes(x = time, y = value, colour = variable, group = variable)) + | |
| geom_line() + | |
| xlab("Time (days)") + # X Axis Label | |
| ylab("Number of people") + # Y Axis Label | |
| labs(title = "SEIRS Model", colour = "Compartments") # Graph Title & Legend Title |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment