##################################### ### voor op NOS-site en RIVM-site ### ##################################### library(deSolve) library(tidyverse) ################################ ### INPUT voor de simulaties ### ################################ ### 2017 population in 10 year age bands (till 80+) agedist <- function() { c(0.10546312, 0.11800504, 0.12623769, 0.12003092, 0.13795451, 0.14478945, 0.12200229, 0.08077414, 0.04474283) } popsize <- function(N = 17180000) return(N) ### contactmatrix of contact hours (Pienter3), and all control matrices. "99" is < 13 March (no contact reduction) # ctrlmats <- readRDS("data/Contactmatrices_19mrt2020.rds") # ## leisure * 80% # ctrlmats <- c(readRDS("data/ContactmatricesD3praktijk_leisure80_24mrt2020.rds"), # readRDS("data/Contactmatrices_19mrt2020.rds")) ## leisure * 70% ctrlmats <- c(readRDS("data/ContactmatricesD3praktijk_leisure70_24mrt2020.rds"), readRDS("data/ContactmatricesD3praktijk_midpoint_24mrt2020.rds")) contactmatrix <- function(scenario = "99", nr = 1, ...) { if(scenario %in% names(ctrlmats)) { return(ctrlmats[[scenario]]) } else if(nr == 0) { return(ctrlmats[["D3praktijk_leisure70"]]) } else { return(ctrlmats[[nr]]) } } ### relative infectiousness by age. "99" is < 13 March set.seed(25032020) rndcorrections <- runif(1000, -1, 1) ctrlinfs <- list( `00` = 1, `99` = 0.946,#0.77,#0.967, IsoMild = 0.98, #0.968*0.83,#0.96, HQMild = 0.86,#0.83*0.84, Elderly = c(1, 1, 0.98, 0.98, 0.98, 0.98, 0.98, 0.885, 0.895), IsoMildElderly = c(0.96, 0.96, 0.94, 0.94, 0.94, 0.94, 0.94, 0.85, 0.86), UK_CI = 0.825, UK_CI_HQ = 0.74, UK_SDO = 1, UK_SD = 1, UK_PC = 1 ) relinfectivity <- function(unequal = FALSE, scenario = "00", nr = 0, ...) { if(unequal) { baseline <- c(0.05, 0.1, 0.25, 0.55, 0.75, 0.9, 0.95, 0.95, 0.95) } else { baseline <- 1 } if(nr == 0) { rndcorrection <- 0 } else if(scenario == "IsoMild") { rndcorrection <- 0.01 * rndcorrections[nr] } else { rndcorrection <- 0.02 * rndcorrections[nr] } return(baseline * ctrlinfs[[scenario]] + rndcorrection) } # default is what we think MRC assumed, based on peak IC rate probsymptomatic <- function(MRC = TRUE, probsymp = 0.5, ...) { if(MRC) { return(rep(probsymp, 9)) } else { return(relinfectivity(TRUE)) } } # de default in deze en volgende is % in OSIRIS vroeg in epidemie probhospitalisation <- function(MRC = TRUE, probhosp = 0.15, ...) { if(MRC) { return(c(0.001, 0.003, 0.012, 0.032, 0.049, 0.102, 0.166, 0.243, 0.273)) } else { return(rep(probhosp, 9)) } } ### IN NL VEEL MINDER IC VOOR 80+ (NICE & OSIRIS) # MRC: (0.243*0.432)/(0.273*0.709) = 0.542 keer 'grotere' kans voor 70-80 naar IC dan voor 80+ # NICE: 149 in 70-80, 29 in 80+. Omgerekend naar demografie: (149/.0808)/(29/.0447) = 2.84 grotere kans. # OSIRIS: 0.092 in 70-80, 0.026 in 80+, dus 0.092/0.026 = 3.54 grotere kans # verhouding tussen MRC en NICE/OSIRIS is 2.84/.542 tot 3.54/.542 = 5.24 tot 6.53. Gebruik nu: 6 probICU <- function(MRC = TRUE, probICU = 0.15, ...) { if(MRC) { return(c(0.05, 0.05, 0.05, 0.05, 0.063, 0.122, 0.274, 0.432, 0.709 / 6)) } else { return(rep(probICU, 9)) } } # natural history parameters parmsNatHist <- function(Rstart = 2.2, SI = 5, ...) { return(list(betaA = Rstart / relinfectivity(scenario = "99") / SI, gammaA = 1 / SI)) } ### delays # modelled as negative binomial distributions, shapes give SDs of 4-5 with default means. shapeI2R <- function(shI2R = 8, ...) return(shI2R) shapeI2H <- function(shI2H = 8, ...) return(shI2H) shapeI2ICU <- function(shI2ICU = 8, ...) return(shI2ICU) # symptom to reporting from OSIRIS meanI2R <- function(meanS2R = 4, incubation = 5, ...) return(meanS2R + incubation) # symptom to hospitalisation from OSIRIS meanI2H <- function(meanS2H = 5, incubation = 5, ...) return(meanS2H + incubation) # symptom to ICU = hosp + 1 dag meanI2ICU <- function(meanS2ICU = 6, incubation = 5, ...) return(meanS2ICU + incubation) # infection to immune = 21 dagen, assumption, not relevant until serosurveys will be carried out meanI2Immune <- function(meanI2Imm = 21, ...) return(meanI2Imm) # durations # modelled as negative binomial distributions, shapes give highly variable stays shapestayhosp <- function(shDhosp = 2, ...) return(shDhosp) shapestayICU <- function(shDICU = 2, ...) return(shDICU) meanstayhosp <- function(meanDhosp = 16, ...) return(meanDhosp) meanstayICU <- function(meanDICU = 10, ...) return(meanDICU) ################# ### modelcode ### ################# # compartments: S, I, and cumulative incidence # I has letter A, to allow more I-compartments later, each with own letter ODEcompartments <- function() return(c("S", "IA", "cumI")) ### rate of new infections into different age classes covidlambda <- function(y, contacts, infectiousness, parms) { res <- rep(0, 9) for(i in ODEcompartments()) { if(substr(i, 1, 1) == "I") { res <- res + parms[[paste0("beta", substr(i, 2, 2))]] * colSums(y[paste0(i, 1:9)] * infectiousness * contacts) * y[paste0("S", 1:9)] } } names(res) <- paste0("la", 1:9) return(res) } ### model equations of SIR model: S, I, and cumulative incidence covidODEmodel <- function(t, y, parms) { la <- covidlambda(y, parms$contacts, parms$infectiousness, parms) res <- rep(0, length(y)) names(res) <- names(y) # dS/dt for(i in 1:9) { res[paste0("S", i)] <- -la[paste0("la", i)] } # dIA/dt for(i in 1:9) { res[paste0("IA", i)] <- la[paste0("la", i)] - parms[["gammaA"]] * y[paste0("IA", i)] } # dcumI/dt for(i in 1:9) { res[paste0("cumI", i)] <- la[paste0("la", i)] } # finished return(list(res)) } ### function for a single simulation, returning time series of all compartments covidsim <- function(contacts = c("99", "99"), infectivities = c("99", "99"), endtimes = c(30, 250), seeding = 10, ...) { ### scale transmission rates relbeta <- 1/eigen(t(contactmatrix(contacts[1], ...) * agedist()) * relinfectivity(...), only.values = T)[[1]][1] nathistparms <- parmsNatHist(...) for(i in names(nathistparms)) { if(substr(i, 1, 4) == "beta") nathistparms[[i]] <- nathistparms[[i]] * relbeta } ### define all compartments compartments <- t(outer(ODEcompartments(), 1:9, paste0)) ### initialize compartments initialstate <- rep(0, length(compartments)) names(initialstate) <- compartments imports <- rep(seeding/9, 9) / popsize() initialstate[paste0("S", 1:9)] <- agedist() - imports initialstate[paste0("IA", 1:9)] <- imports ### simulate no control ressimul <- ode( func = covidODEmodel, y = initialstate, times = seq(0, endtimes[1], 1), parms = c(list(contacts = contactmatrix(contacts[1], ...), infectiousness = relinfectivity(scenario = infectivities[1], ...)), nathistparms) ) ### simulate controls for(i in 2:length(contacts)) { ressimul <- rbind( head(ressimul, -1), ode( func = covidODEmodel, y = tail(ressimul, 1)[, compartments], times = seq(endtimes[i - 1], endtimes[i], 1), parms = c(list(contacts = contactmatrix(contacts[i], ...), infectiousness = relinfectivity(scenario = infectivities[i], ...)), nathistparms) ) ) } return(ressimul) } ### function for a single simulation, starting 12 Feb (day 30 = 13 March, start of control), # returning observed numbers of cases in age classes: # * symptomatic incidence (could be reported) # * hospitalisation incidence # * ICU incidence # * immunity prevalence # * hospitalisation prevalence # * ICU prevalence covidcontrolsim <- function(contacts = c("99", "99"), infectivities = c("99", "99"), endtimes = c(30, 250), seeding = 10, ...) { result <- covidsim(contacts = contacts, infectivities = infectivities, endtimes = endtimes, seeding = seeding, ...) # get daily incidence (new infections) cumincidencecurves <- result[, paste0("cumI", 1:9)] * popsize() + seeding/9 incidencecurves <- rbind(rep(seeding/9/popsize(), 9) , tail(cumincidencecurves, -1) - head(cumincidencecurves, -1)) # calculate numbers that will be observed (after delay) tobesymptomaticcurves <- t(incidencecurves) * probsymptomatic(...) tobehospitalisedcurves <- tobesymptomaticcurves * probhospitalisation(...) tobeICUcurves <- tobehospitalisedcurves * probICU(...) # calculate observation incidences (applying delays) incsymptomaticcurves <- sapply(1:(1 + tail(endtimes, 1)), function(x) colSums(t(tobesymptomaticcurves[, 1:x, drop = FALSE]) * dnbinom((x-1):0, shapeI2R(...), mu = meanI2R(...)))) inchospitalisedcurves <- sapply(1:(1 + tail(endtimes, 1)), function(x) colSums(t(tobehospitalisedcurves[, 1:x, drop = FALSE]) * dnbinom((x-1):0, shapeI2H(...), mu = meanI2H(...)))) incICUcurves <- sapply(1:(1 + tail(endtimes, 1)), function(x) colSums(t(tobeICUcurves[, 1:x, drop = FALSE]) * dnbinom((x-1):0, shapeI2ICU(...), mu = meanI2ICU(...)))) immunecurves <- sapply(1:(1 + tail(endtimes, 1)), function(x) colSums(cumincidencecurves[1:x, , drop = FALSE] * dpois((x-1):0, meanI2Immune()))) # calculate observation prevalences (applying durations) prevhospitalisedcurves <- sapply(1:(1 + tail(endtimes, 1)), function(x) colSums(t(inchospitalisedcurves[, max(1, x - meanstayhosp(...) + 1):x, drop = FALSE]))) prevICUcurves <- sapply(1:(1 + tail(endtimes, 1)), function(x) colSums(t(incICUcurves[, max(1, x - meanstayICU(...) + 1):x, drop = FALSE]))) # combine and add names toreturn <- as_tibble(t(incsymptomaticcurves)) %>% rename_all(~rownames(contactmatrix())) %>% mutate(time = 0:tail(endtimes, 1)) %>% pivot_longer(1:9, names_to = "ageclass", values_to = "incsymp") %>% full_join( as_tibble(t(inchospitalisedcurves)) %>% rename_all(~rownames(contactmatrix())) %>% mutate(time = 0:tail(endtimes, 1)) %>% pivot_longer(1:9, names_to = "ageclass", values_to = "inchosp"), by = c("time", "ageclass") ) %>% full_join( as_tibble(t(incICUcurves)) %>% rename_all(~rownames(contactmatrix())) %>% mutate(time = 0:tail(endtimes, 1)) %>% pivot_longer(1:9, names_to = "ageclass", values_to = "incICU"), by = c("time", "ageclass") ) %>% full_join( as_tibble(t(immunecurves)) %>% rename_all(~rownames(contactmatrix())) %>% mutate(time = 0:tail(endtimes, 1)) %>% pivot_longer(1:9, names_to = "ageclass", values_to = "immune"), by = c("time", "ageclass") ) %>% full_join( as_tibble(t(prevhospitalisedcurves)) %>% rename_all(~rownames(contactmatrix())) %>% mutate(time = 0:tail(endtimes, 1)) %>% pivot_longer(1:9, names_to = "ageclass", values_to = "prevhosp"), by = c("time", "ageclass") ) %>% full_join( as_tibble(t(prevICUcurves)) %>% rename_all(~rownames(contactmatrix())) %>% mutate(time = 0:tail(endtimes, 1)) %>% pivot_longer(1:9, names_to = "ageclass", values_to = "prevICU"), by = c("time", "ageclass") ) return(toreturn) } ### Calibration with control up to "today" # change 'seeding' until current ICU patients # are similar to observations in NICE # On 24 March: approx 560 known. ijksimulatie <- covidcontrolsim(contacts = c("99", "RND", "RND"), infectivities = c("99", "IsoMild", "HQMild"), endtimes = c(30, 41, 400), seeding = 61, unequal = F, MRC = T, Rstart = 2.2, meanDICU = 23, nr = 1) ijksimulatie <- covidcontrolsim(contacts = c("99", "99", "99"), infectivities = c("99", "99", "99"), endtimes = c(30, 41, 400), seeding = 61, unequal = F, MRC = T, Rstart = 2.2, meanDICU = 23, nr = 1) ijksimulatie %>% group_by(time) %>% summarise(totprevICU = sum(prevICU)) %>% filter(time == 41) ijksimulatie %>% filter(time <= 30) %>% summarise(totinchosp = sum(inchosp)) ijksimulatie %>% filter(time <= 30) %>% summarise(totincsymp = sum(incsymp)) panels_rijen <- c(NA) panels_kolommen <- c(NA) panels_kleuren <- c("geen bestrijding", "huidige pakket", "thuisquarantaine")#, "lockdown") contactkeus <- function(rij = 1, kolom = 1, kleur = 1) { case_when( kleur == 1 ~ c("99", "99", "99"), kleur == 2 ~ c("99", "RND", "RND"), kleur == 3 ~ c("99", "RND", "RND") ) } infectiviteitkeus <- function(rij = 1, kolom = 1, kleur = 1) { case_when( kleur == 1 ~ c("99", "99", "99"), kleur == 2 ~ c("99", "IsoMild", "IsoMild"), kleur == 3 ~ c("99", "IsoMild", "HQMild") ) } eindtijden <- function() c(30, 41, 400) allsimulresults <- list() set.seed(25032020) Rstarts <- c(2.2, runif(1000, 2, 2.4)) doublings <- c(2.9, runif(1000, 2.8, 3)) allfigdata <- list() for(repnr in 0:1000) { allsimulresults <- list() for(i in 1:1) { for(j in 1:1) { for(k in 1:3) { simul <- covidcontrolsim(contacts = contactkeus(kleur = k), infectivities = infectiviteitkeus(kleur = k), endtimes = eindtijden(), unequal = F, MRC = T, Rstart = Rstarts[repnr + 1], SI = (Rstarts[repnr + 1] - 1) * doublings[repnr + 1]/log(2), seeding = 61, nr = repnr, meanDICU = 23) allsimulresults <- c(allsimulresults, list(simul %>% mutate(rij = panels_rijen[i], kolom = panels_kolommen[j], kleur = panels_kleuren[k])) ) } } } allfigdata <- c(allfigdata, list(do.call(rbind, allsimulresults) %>% group_by(time, rij, kolom, kleur) %>% summarise(ICUbezetting = sum(prevICU)) %>% ungroup() %>% mutate(tijd = as.Date("2020-02-12") + time) %>% select(time, tijd, ICUbezetting, kleur) %>% mutate(replicate = repnr) )) print(repnr) } saveRDS(allfigdata, "results/allfigdata_25032020") # saveRDS(allfigdata, "results/allfigdata_leisure70") # saveRDS(allfigdata, "results/allfigdata_leisure80") allfigdata <- readRDS("results/allfigdata_leisure70_IC20") figuredata <- do.call(rbind, allfigdata) %>% rename(repl = `replicate`) figuredata %>% filter(time >= 25 & time < 250) %>% mutate(repl = repl + 1000 * (kleur == "huidige pakket") + 2000 * (kleur == "thuisquarantaine")) %>% # filter(kleur == "geen bestrijding") %>% # group_by(tijd, kleur) %>% ggplot(aes(x = tijd, y = ICUbezetting, group = repl, color = kleur)) + geom_line(size = 0.1) + theme_light() + labs(x = "Dag", y = "Aantal bezette IC-plaatsen") + # scale_x_date(date_labels = "%d-%m") + # scale_x_date(date_breaks = "1 week", date_labels = "%d-%m") + scale_y_continuous(limits = c(0,5000)) + scale_color_discrete(name = "scenario") + guides(color = guide_legend(override.aes = list(size = 1))) + # facet_wrap(~ maatregelen, ncol = 1) + geom_hline(aes(yintercept = 1208), color = "black", linetype = "dashed") + annotate(geom = "text", x = as.Date("2020-03-07"), y = 1300, label = "Nederlandse IC-capaciteit", hjust = 0, size = 3) # ggsave("results/sims_spaghettileisure70_ICU20.tiff", units = "mm", width = 180, height = 120, dpi = 150) # eigenvalues <- sapply(1:1000, # function(x) eigen(t(ctrlmats[[x]] * agedist()) * 0.9525 * relinfectivity(F)/.946, only.values = T)[[1]][1]) # # figuredata %>% # filter(time >= 25 & time < 250 & kleur == "huidige pakket") %>% # mutate( # R0 = Rstarts[repl], # control = eigenvalues[repl], # lambda = doublings[repl]) %>% # # filter(kleur == "geen bestrijding") %>% # # group_by(tijd, kleur) %>% # ggplot(aes(x = tijd, y = ICUbezetting, group = repl, color = lambda)) + # geom_line(size = 0.1) + # theme_light() + # labs(x = "Dag", y = "Aantal bezette IC-plaatsen") + # # scale_x_date(date_labels = "%d-%m") + # # scale_x_date(date_breaks = "1 week", date_labels = "%d-%m") + # scale_y_continuous(limits = c(0,5000)) + # scale_color_gradient2(name = "doubling\ntime", midpoint = 2.9, low = "red", high = "blue") + # # guides(color = guide_legend(override.aes = list(size = 1))) + # # facet_wrap(~ maatregelen, ncol = 1) + # geom_hline(aes(yintercept = 1208), color = "black", linetype = "dashed") + # annotate(geom = "text", x = as.Date("2020-03-07"), y = 1300, # label = "Nederlandse IC-capaciteit", hjust = 0, size = 3) # ggsave("results/sims_spaghettileisure70_ICU20_doubling.tiff", units = "mm", width = 180, height = 120, dpi = 150) # figuredata %>% # filter(time >= 25 & time < 250 & kleur == "huidige pakket") %>% # mutate( # R0 = Rstarts[repl], # control = eigenvalues[repl], # lambda = doublings[repl]) %>% # # filter(kleur == "geen bestrijding") %>% # # group_by(tijd, kleur) %>% # ggplot(aes(x = tijd, y = ICUbezetting, group = repl, color = R0)) + # geom_line(size = 0.1) + # theme_light() + # labs(x = "Dag", y = "Aantal bezette IC-plaatsen") + # # scale_x_date(date_labels = "%d-%m") + # # scale_x_date(date_breaks = "1 week", date_labels = "%d-%m") + # scale_y_continuous(limits = c(0,5000)) + # scale_color_gradient2(name = "Reproduction\nnumber", midpoint = 2.2, high = "red", low = "blue") + # # guides(color = guide_legend(override.aes = list(size = 1))) + # # facet_wrap(~ maatregelen, ncol = 1) + # geom_hline(aes(yintercept = 1208), color = "black", linetype = "dashed") + # annotate(geom = "text", x = as.Date("2020-03-07"), y = 1300, # label = "Nederlandse IC-capaciteit", hjust = 0, size = 3) # ggsave("results/sims_spaghettileisure70_ICU20_reprratio.tiff", units = "mm", width = 180, height = 120, dpi = 150) # figuredata %>% # filter(time >= 25 & time < 250 & kleur == "huidige pakket") %>% # mutate( # R0 = Rstarts[repl], # control = eigenvalues[repl], # lambda = doublings[repl]) %>% # # filter(kleur == "geen bestrijding") %>% # # group_by(tijd, kleur) %>% # ggplot(aes(x = tijd, y = ICUbezetting, group = repl, color = control)) + # geom_line(size = 0.1) + # theme_light() + # labs(x = "Dag", y = "Aantal bezette IC-plaatsen") + # # scale_x_date(date_labels = "%d-%m") + # # scale_x_date(date_breaks = "1 week", date_labels = "%d-%m") + # scale_y_continuous(limits = c(0,5000)) + # scale_color_gradient2(name = "Contact\nreduction", midpoint = .450, low = "blue", high = "red") + # # guides(color = guide_legend(override.aes = list(size = 1))) + # # facet_wrap(~ maatregelen, ncol = 1) + # geom_hline(aes(yintercept = 1208), color = "black", linetype = "dashed") + # annotate(geom = "text", x = as.Date("2020-03-07"), y = 1300, # label = "Nederlandse IC-capaciteit", hjust = 0, size = 3) # ggsave("results/sims_spaghettileisure70_ICU20_control.tiff", units = "mm", width = 180, height = 120, dpi = 150) median(eigenvalues) distsdata <- allsimulsummaries[[1]] distsdata %>% # filter(scenario %in% c("zonder bestrijding", "doorgaan met huidige pakket")) %>% ggplot(aes(x = agecat, y = immuniteit, fill = scenario)) + geom_col(position = position_dodge()) ggsave("results/sims_packhqlock_immunity.tiff", units = "mm", width = 180, height = 120, dpi = 150) distsdata %>% # filter(scenario %in% c("zonder bestrijding", "doorgaan met huidige pakket")) %>% ggplot(aes(x = agecat, y = ICopnames, fill = scenario)) + geom_col(position = position_dodge()) ggsave("results/sims_packhqlock_ICU.tiff", units = "mm", width = 180, height = 120, dpi = 150) figuredata %>% filter(time >= 25 & time < 200 & kleur != "huidige pakket") %>% mutate(kleur = if_else(kleur == "thuisquarantaine", "met gevolgde\n maatregelen", kleur)) %>% group_by(tijd, kleur) %>% summarise( lowerbound = quantile(ICUbezetting, .025), upperbound = quantile(ICUbezetting, .975), midbound = median(ICUbezetting) ) %>% ungroup() %>% ggplot(aes(x = tijd, y = midbound, color = kleur)) + geom_line(size = 2) + geom_ribbon(aes(ymin = lowerbound, ymax = upperbound, fill = kleur, color = NULL), alpha = 0.3) + theme_light() + labs(x = "Dag", y = "Aantal bezette IC-plaatsen") + # scale_x_date(date_labels = "%d-%m") + scale_color_discrete(name = "scenario") + scale_fill_discrete(name = "scenario") + scale_x_date(date_breaks = "1 month", date_labels = "1 %b") + coord_cartesian(ylim = c(0,3000)) + # guides(color = guide_legend(override.aes = list(size = 1))) + # facet_wrap(~ maatregelen, ncol = 1) + geom_hline(aes(yintercept = 1208), color = "black", linetype = "dashed") + annotate(geom = "text", x = as.Date("2020-07-15"), y = 1280, label = "Nederlandse IC-capaciteit", hjust = 0, size = 3) ggsave("results/sims_25032020.tiff", units = "mm", width = 200, height = 120, dpi = 150) ggsave("results/sims_25032020.pdf", units = "mm", width = 200, height = 120, dpi = 150) ggsave("results/sims_25032020.jpg", units = "mm", width = 200, height = 120, dpi = 300) # hoogte van de piek figuredata %>% filter(kleur == "thuisquarantaine") %>% group_by(repl) %>% summarise(maxIC = max(ICUbezetting)) %>% pull(maxIC) %>% quantile(probs = c(.025,.5,.975)) # locatie van de piek figuredata %>% filter(kleur == "thuisquarantaine") %>% group_by(repl) %>% summarise(peaktime = time[ICUbezetting == max(ICUbezetting)]) %>% pull(peaktime) %>% quantile(probs = c(.025,.5,.975)) %>% `+`(as.Date("2020-02-12"))