diff --git a/app/R/custom_graphics.R b/app/R/custom_graphics.R
index a58bab9..0c3d5b0 100644
--- a/app/R/custom_graphics.R
+++ b/app/R/custom_graphics.R
@@ -60,3 +60,7 @@ cmmid_pal <- function(n){
   colorRampPalette(colors = c("#388a68","#0D5257"))(n)
   
 }
+
+my_palette <- c(Unreported = scales::alpha(cmmid_color, 0.5),
+                Reported = cmmid_color,
+                Projected = "#00AEC7")
diff --git a/app/R/dhlabel.R b/app/R/dhlabel.R
new file mode 100644
index 0000000..b81dec1
--- /dev/null
+++ b/app/R/dhlabel.R
@@ -0,0 +1,7 @@
+dhlabel <- function(x){
+    if (x == 1){
+        "doubling"
+    } else {
+        "halving"
+    }
+}
\ No newline at end of file
diff --git a/app/R/gamma_q2shapescale.R b/app/R/gamma_q2shapescale.R
new file mode 100644
index 0000000..7ce457b
--- /dev/null
+++ b/app/R/gamma_q2shapescale.R
@@ -0,0 +1,256 @@
+gamma_q2shapescale <- function(q, p=c(0.025,0.975),
+                                       precision=0.001, derivative.epsilon=1e-3, start.with.normal.approx=T, start=c(1.1, 1.1), plot=F, plot.xlim=numeric(0))
+{
+    # Version 1.0.2 (December 2012)
+    #
+    # Function developed by 
+    # Lawrence Joseph and Patrick Bélisle
+    # Division of Clinical Epidemiology
+    # Montreal General Hospital
+    # Montreal, Qc, Can
+    #
+    # patrick.belisle@rimuhc.ca
+    # http://www.medicine.mcgill.ca/epidemiology/Joseph/PBelisle/GammaParmsFromQuantiles.html
+    #
+    # Please refer to our webpage for details on each argument.
+    
+    f <- function(x, theta){dgamma(x, shape=theta[1], rate=theta[2])}
+    F.inv <- function(x, theta){qgamma(x, shape=theta[1], rate=theta[2])}
+    f.cum <- function(x, theta){pgamma(x, shape=theta[1], rate=theta[2])}
+    f.mode <- function(theta){shape <- theta[1]; rate <- theta[2]; mode <- ifelse(shape>1,(shape-1)/rate,NA); mode}
+    theta.from.moments <- function(m, v){shape <- m*m/v; rate <- m/v; c(shape, rate)}
+    dens.label <- 'dgamma'
+    parms.names <- c('shape', 'rate')
+    
+    
+    if (length(p) != 2) stop("Vector of probabilities p must be of length 2.")
+    if (length(q) != 2) stop("Vector of quantiles q must be of length 2.")
+    p <- sort(p); q <- sort(q)
+    
+    #_____________________________________________________________________________________________________
+    
+    print.area.text <- function(p, p.check, q, f, f.cum, F.inv, theta, mode, cex, plot.xlim, M=30, M0=50)
+    {
+        par.usr <- par('usr')
+        par.din <- par('din')
+        
+        p.string <- as.character(round(c(0,1) + c(1,-1)*p.check, digits=4))
+        str.width <- strwidth(p.string, cex=cex)
+        str.height <- strheight("0", cex=cex)
+        
+        J <- matrix(1, nrow=M0, ncol=1)
+        
+        x.units.1in <- diff(par.usr[c(1,2)])/par.din[1]
+        y.units.1in <- diff(par.usr[c(3,4)])/par.din[2]
+        aspect.ratio <- y.units.1in/x.units.1in
+        
+        # --- left area  -----------------------------------------------------------
+        
+        scatter.xlim <- c(max(plot.xlim[1], par.usr[1]), q[1])
+        scatter.ylim <- c(0, par.usr[4])
+        x <- seq(from=scatter.xlim[1], to=scatter.xlim[2], length=M)
+        y <- seq(from=scatter.ylim[1], to=scatter.ylim[2], length=M)
+        x.grid.index <- rep(seq(M), M)
+        y.grid.index <- rep(seq(M), rep(M, M))
+        
+        grid.df <- f(x, theta)
+        
+        # Estimate mass center
+        tmp.p <- seq(from=0, to=p[1], length=M0)
+        tmp.x <- F.inv(tmp.p, theta)
+        h <- f(tmp.x, theta)
+        mass.center <- c(mean(tmp.x), sum(h[-1]*diff(tmp.x))/diff(range(tmp.x)))
+        
+        # Identify points under the curve
+        # (to eliminate them from the list of candidates)
+        gridpoint.under.the.curve <- y[y.grid.index] <= grid.df[x.grid.index]
+        w <- which(gridpoint.under.the.curve)
+        x <- x[x.grid.index]; y <- y[y.grid.index]
+        if (length(w)){x <- x[-w]; y <- y[-w]}
+        
+        # Eliminate points to the right of the mode, if any
+        w <- which(x>mode)
+        if (length(w)){x <- x[-w]; y <- y[-w]}
+        
+        # Eliminate points for which the text would fall out of the plot area
+        w <- which((par.usr[1]+str.width[1]) <= x & (y + str.height) <= par.usr[4])
+        x <- x[w]; y <- y[w]
+        
+        # For each height, eliminate the closest point to the curve
+        # (we want to stay away from the curve to preserve readability)
+        w <- which(!duplicated(y, fromLast=T))
+        if (length(w)){x <- x[-w]; y <- y[-w]}
+        
+        # For each point, compute distance from mass center and pick the closest point
+        d <- ((x-mass.center[1])^2) + ((y-mass.center[2])/aspect.ratio)^2
+        w <- which.min(d)
+        x <- x[w]; y <- y[w]
+        
+        if (length(x))
+        {
+            text(x, y, labels=p.string[1], adj=c(1,0), col='gray', cex=cex)
+        }
+        else
+        {
+            text(plot.xlim[1], mean(par.usr[c(3,4)]), labels=p.string[1], col='gray', cex=cex, srt=90, adj=c(1,0))
+        }
+        
+        # --- right area  ----------------------------------------------------------
+        
+        scatter.xlim <- c(q[2], plot.xlim[2])
+        scatter.ylim <- c(0, par.usr[4])
+        x <- seq(from=scatter.xlim[1], to=scatter.xlim[2], length=M)
+        y <- seq(from=scatter.ylim[1], to=scatter.ylim[2], length=M)
+        x.grid.index <- rep(seq(M), M)
+        y.grid.index <- rep(seq(M), rep(M, M))
+        grid.df <- f(x, theta)
+        
+        # Estimate mass center
+        tmp.p <- seq(from=p[2], to=f.cum(plot.xlim[2], theta), length=M0)
+        tmp.x <- F.inv(tmp.p, theta)
+        h <- f(tmp.x, theta)
+        mass.center <- c(mean(tmp.x), sum(h[-length(h)]*diff(tmp.x))/diff(range(tmp.x)))
+        
+        # Identify points under the curve
+        # (to eliminate them from the list of candidates)
+        gridpoint.under.the.curve <- y[y.grid.index] <= grid.df[x.grid.index]
+        w <- which(gridpoint.under.the.curve)
+        x <- x[x.grid.index]; y <- y[y.grid.index]
+        if (length(w)){x <- x[-w]; y <- y[-w]}
+        
+        # Eliminate points to the left of the mode, if any
+        w <- which(x<mode)
+        if (length(w)){x <- x[-w]; y <- y[-w]}
+        
+        # Eliminate points for which the text would fall out of the plot area
+        w <- which((par.usr[2]-str.width[2]) >= x & (y + str.height) <= par.usr[4])
+        x <- x[w]; y <- y[w]
+        
+        # For each height, eliminate the closest point to the curve
+        # (we want to stay away from the curve to preserve readability)
+        w <- which(!duplicated(y))
+        if (length(w)){x <- x[-w]; y <- y[-w]}
+        
+        # For each point, compute distance from mass center and pick the closest point
+        d <- ((x-mass.center[1])^2) + ((y-mass.center[2])/aspect.ratio)^2
+        w <- which.min(d)
+        x <- x[w]; y <- y[w]
+        
+        if (length(x))
+        {
+            text(x, y, labels=p.string[2], adj=c(0,0), col='gray', cex=cex)
+        }
+        else
+        {
+            text(plot.xlim[2], mean(par.usr[c(3,4)]), labels=p.string[2], col='gray', cex=cex, srt=-90, adj=c(1,0))
+        }
+    }
+    
+    # ......................................................................................................................................
+    
+    Newton.Raphson <- function(derivative.epsilon, precision, f.cum, p, q, theta.from.moments, start.with.normal.approx, start)
+    {
+        Hessian <- matrix(NA, 2, 2)
+        
+        if (start.with.normal.approx)
+        {
+            # Probably not a very good universal choice, but proved good in most cases in practice
+            m <-  diff(q)/diff(p)*(0.5-p[1]) + q[1]
+            v <- (diff(q)/diff(qnorm(p)))^2
+            theta <- theta.from.moments(m, v)
+        }
+        else theta <- start
+        
+        
+        change <- precision + 1
+        niter <- 0
+        # Newton-Raphson multivariate algorithm
+        while (max(abs(change)) > precision)
+        {
+            Hessian[,1] <- (f.cum(q, theta) - f.cum(q, theta - c(derivative.epsilon, 0))) / derivative.epsilon
+            Hessian[,2] <- (f.cum(q, theta) - f.cum(q, theta - c(0, derivative.epsilon))) / derivative.epsilon
+            
+            f <- f.cum(q, theta) - p
+            change <- solve(Hessian) %*% f
+            last.theta <- theta
+            theta <- last.theta - change
+            
+            # If we step out of limits, reduce change
+            
+            if (any(theta<0))
+            {
+                k <- min(last.theta/change)
+                theta <- last.theta - k/2*change
+            }
+            
+            niter <- niter + 1
+        }
+        
+        list(theta=as.vector(theta), niter=niter, last.change=as.vector(change))
+    }
+    
+    # ...............................................................................................................
+    
+    plot.density <- function(p, q, f, f.cum, F.inv, mode, theta, plot.xlim, dens.label, parms.names, cex)
+    {
+        if (length(plot.xlim) == 0)
+        {
+            plot.xlim <- F.inv(c(0, 1), theta)
+            
+            if (is.infinite(plot.xlim[1]))
+            {
+                tmp <- min(c(0.001, p[1]/10))
+                plot.xlim[1] <- F.inv(tmp, theta)
+            }  
+            
+            if (is.infinite(plot.xlim[2]))
+            {
+                tmp <- max(c(0.999, 1 - (1-p[2])/10))
+                plot.xlim[2] <- F.inv(tmp, theta)
+            }
+        }
+        plot.xlim <- sort(plot.xlim)
+        
+        
+        x <- seq(from=min(plot.xlim), to=max(plot.xlim), length=1000)
+        h <- f(x, theta)
+        x0 <- x; f0 <- h
+        ylab <- paste(c(dens.label, '(x, ', parms.names[1], ' = ', round(theta[1], digits=5), ', ', parms.names[2], ' = ', round(theta[2], digits=5), ')'), collapse='')
+        plot(x, h, type='l', ylab=ylab)
+        
+        # fill in area on the left side of the distribution
+        x <- seq(from=plot.xlim[1], to=q[1], length=1000)
+        y <- f(x, theta)
+        x <- c(x, q[1], plot.xlim[1]); y <- c(y, 0, 0)
+        polygon(x, y, col='lightgrey', border='lightgray')
+        # fill in area on the right side of the distribution
+        x <- seq(from=max(plot.xlim), to=q[2], length=1000)
+        y <- f(x, theta)
+        x <- c(x, q[2], plot.xlim[2]); y <- c(y, 0, 0)
+        polygon(x, y, col='lightgrey', border='lightgray')
+        # draw distrn again
+        points(x0, f0, type='l')
+        h <- f(q, theta)
+        points(rep(q[1], 2), c(0, h[1]), type='l', col='orange')
+        points(rep(q[2], 2), c(0, h[2]), type='l', col='orange')
+        # place text on both ends areas
+        print.area.text(p, p.check, q, f, f.cum, F.inv, theta, mode, cex, plot.xlim)  
+        
+        xaxp <- par("xaxp")
+        x.ticks <- seq(from=xaxp[1], to=xaxp[2], length=xaxp[3]+1)
+        q2print <- as.double(setdiff(as.character(q), as.character(x.ticks)))
+        
+        mtext(q2print, side=1, col='orange', at=q2print, cex=0.6, line=2.1)
+        points(q, rep(par('usr')[3]+0.15*par('cxy')[2], 2), pch=17, col='orange')
+    }
+    
+    #________________________________________________________________________________________________________________
+    
+    
+    parms <- Newton.Raphson(derivative.epsilon, precision, f.cum, p, q, theta.from.moments, start.with.normal.approx, start=start)
+    p.check <- f.cum(q, parms$theta)
+    
+    if (plot) plot.density(p, q, f, f.cum, F.inv, f.mode(parms$theta), parms$theta, plot.xlim, dens.label, parms.names, 0.8)
+    
+    list(shape=parms$theta[1], rate=parms$theta[2], scale=1/parms$theta[2], last.change=parms$last.change, niter=parms$niter, q=q, p=p, p.check=p.check)
+}
diff --git a/app/R/los_distributions.R b/app/R/los_distributions.R
index 957d5e9..195aef0 100644
--- a/app/R/los_distributions.R
+++ b/app/R/los_distributions.R
@@ -29,44 +29,72 @@ weibull_mucv2shapescale <- function(mean, cv){
 ## we want to avoid having 0 days LoS, so what we do is actually:
 ## - generate a discretised Gamma with (mean - 1)
 ## - add 1 to simulated LoS
-los_dist <- function(distribution = "gamma", mean, cv) {
+los_dist <- function(distribution = "gamma", q) {
   
   if (distribution == "gamma"){
-    params <- epitrix::gamma_mucv2shapescale(mean - 1, cv)
+    if (all(q == min(q))){
+      params <- epitrix::gamma_mucv2shapescale(mu = min(q), cv = 1e-08)
+    } else {
+      params <- gamma_q2shapescale(q, p = c(0.25, 0.75))
+    }
     auxil <- distcrete::distcrete("gamma",
                                   shape = params$shape,
                                   scale = params$scale,
                                   w = 0.5, interval = 1)
-    r <- function(n) auxil$r(n) + 1
+    r <- function(n) auxil$r(n)
     d <- function(x) {
-      out <- auxil$d(x - 1)
-      out[x < 1] <- 0
+      out <- auxil$d(x)
+      #out[x < 1] <- 0
       out
     }
     q <- function(x) {
-      1 + auxil$q(x)
+      auxil$q(x)
     }
   }
   
   if (distribution == "weibull"){
     
+    if (all(q == min(q))){
+      params <- list(shape = Inf, scale = min(q))
+    } else {
+      params <- weibull_q2shapescale(q, p = c(0.25, 0.75))  
+    }
     
-    
-    params <- weibull_mucv2shapescale(mean, cv)
     auxil <- distcrete::distcrete(distribution,
                                   shape = params$shape,
                                   scale = params$scale,
                                   w = 0.5, interval = 1)
     
-    r <- function(n) auxil$r(n)
+    r <- function(n) auxil$r(n) 
     d <- function(x) {
-      auxil$d(x)
+      out <- auxil$d(x)
+      #out[x < 1] <- 0
+      out
     }
     q <- function(x) {
-      auxil$q(x)
+      auxil$q(x) 
     }
     
   }
   
   list(r = r, d = d, q = q)
 }
+
+los_params <- function(distribution = "gamma", q) {
+  
+  if (all(q == min(q))){
+    params <- epitrix::gamma_mucv2shapescale(mu = min(q), cv = 1e-08)
+  } else {
+    params <- gamma_q2shapescale(q, p = c(0.25, 0.75))
+  }
+  
+  if (distribution == "weibull"){
+    if (all(q == min(q))){
+      params <- list(shape = Inf, scale = min(q))
+    } else {
+      params <- weibull_q2shapescale(q, p = c(0.25, 0.75))
+    }
+  }
+  
+  as.list(params)
+}
diff --git a/app/R/make_los_parameters.R b/app/R/make_los_parameters.R
index aad0d5a..0059bc4 100644
--- a/app/R/make_los_parameters.R
+++ b/app/R/make_los_parameters.R
@@ -1,19 +1,31 @@
+
+
 los_parameters <- do.call(
     what = rbind,
     args = list(
-        `Custom` = 
-            data.frame(mean_los =  7,
-                       cv_los   =  0.1,
-                       los_dist = "gamma"),
+        `Rees and Nightingale et al. non-critical care` = 
+            data.frame(los_25 = 3,
+                       los_75 = 9,
+                       los_dist = "weibull"),
+        `Rees and Nightingale et al. critical care` =
+            data.frame(los_25 = 4,
+                       los_75 = 11,
+                       los_dist = "weibull"),
         `Zhou et al. non-critical care` = 
-            data.frame(mean_los = 11.52095,
-                       cv_los   =  0.5227232,
+            data.frame(los_25 = 7,
+                       los_75 = 14,
                        los_dist = "weibull"),
         `Zhou et al. critical care` =
-            data.frame(mean_los =  8.862269,
-                       cv_los   =  0.5227232,
-                       los_dist = "weibull")
+            data.frame(los_25 = 4,
+                       los_75 = 12,
+                       los_dist = "weibull"),
+        `Custom` = 
+            data.frame(los_25   =  5,
+                       los_75   =  13,
+                       los_dist = "gamma")
     ))
 
-los_parameters$name <- rownames(los_parameters)
+
+
+los_parameters$name     <- rownames(los_parameters)
 los_parameters$los_dist <- as.character(los_parameters$los_dist)
diff --git a/app/R/make_r0.R b/app/R/make_r0.R
new file mode 100644
index 0000000..cc4f408
--- /dev/null
+++ b/app/R/make_r0.R
@@ -0,0 +1,14 @@
+make_r0 <- function(n, mean, cv) {
+    if (cv == 0){
+        R0 <- rep(mean, n)
+        return(list(R0 = R0))
+    } else {
+        R0_parms <- epitrix::gamma_mucv2shapescale(mean, cv)
+        R0 <- rgamma(n     = n,
+                     shape = R0_parms$shape,
+                     scale = R0_parms$scale)
+        return(list(R0 = R0, R0_parms = R0_parms))
+    }
+    
+    
+}
diff --git a/app/R/make_si.R b/app/R/make_si.R
new file mode 100644
index 0000000..da44a8a
--- /dev/null
+++ b/app/R/make_si.R
@@ -0,0 +1,8 @@
+make_si <- function(mean, cv){
+    
+    parms <- epitrix::gamma_mucv2shapescale(mean, cv)
+    si <- distcrete::distcrete("gamma", shape = parms$shape, scale = parms$scale, interval = 1, w = 0)
+    
+    si
+    
+}
diff --git a/app/R/plot_admissions.R b/app/R/plot_admissions.R
new file mode 100644
index 0000000..7dc9722
--- /dev/null
+++ b/app/R/plot_admissions.R
@@ -0,0 +1,64 @@
+#' Function to plot discretised duration distribution
+#'
+#' This function provides custom plots for the output of `distcrete`.
+#'
+#' @param los a length of stay distribution of class `distcrete`
+#' @param title a string to be added to the resulting ggplot
+#' 
+#' @author Sam Clifford
+#' 
+
+plot_admissions <- function(data,
+                            projections, 
+                            reporting = 100, 
+                            title = NULL) {
+    
+    
+    
+    projections <- summarise_beds(projections[-c(1:nrow(data)),])
+    
+    projections$Status <- "Projected"
+    projections$n_admissions <- NA
+    projections$date <- projections$Date
+    
+    # my_palette <- c(Reported = cmmid_color,
+    #                 Unreported = alpha(cmmid_color, 0.5),
+    #                 Projected = "#00AEC7")
+    
+    palette_to_use <- c("Reported", 
+                        "Unreported",
+                        "Projected")[
+                            c(TRUE,
+                              reporting < 100,
+                              TRUE)]
+    
+    out_plot <- plot_data(data = data, 
+                          reporting = reporting,
+                          title = title)
+    out_plot +
+        ggplot2::geom_col(data = projections,
+                          width = 0.8,
+                          aes(fill = Status,
+                              y = Median)) +
+        ggplot2::geom_segment(data = projections,
+                              aes(y = `lower 95%`,
+                                  xend = date,
+                                  yend = `upper 95%`),
+                              size = 1
+        ) +
+        # ggplot2::geom_segment(data = projections,
+        #                       aes(y = `lower 50%`,
+        #                           xend = date,
+        #                           yend = `upper 50%`),
+        #                       size = 2,
+        #                       color = cmmid_pal(1)
+        # ) +
+        ggplot2::scale_fill_manual(values = my_palette[palette_to_use],
+                                   breaks = palette_to_use,
+                                   name = "Reporting status") +
+        ggplot2::theme(legend.position = "bottom") +
+        ggplot2::scale_y_continuous(limits = c(0, NA), breaks = int_breaks) +
+        ggplot2::scale_x_date(date_label = "%d %b %y")
+    
+    
+}
diff --git a/app/R/plot_branching.R b/app/R/plot_branching.R
new file mode 100644
index 0000000..634a4f6
--- /dev/null
+++ b/app/R/plot_branching.R
@@ -0,0 +1,64 @@
+#' Functions for plotting branching process parameters
+#'
+#' 
+#' @param R A `list` containing `R0`, a vector of sampled basic reproduction
+#'   numbers, and `R0_parms`, a list containing the `shape` and `scale` parameters
+#'   used to generate these samples
+#'
+#' 
+#' @author Sam Clifford
+
+
+plot_r0 <- function(R){
+    
+    
+    
+    if (sd(R$R0) == 0){
+        r0_plot <- ggplot2::ggplot(data = data.frame(R = unique(R$R0),
+                                                     y = 1)) +
+            ggplot2::geom_segment(aes(x = R,
+                                      xend = R,
+                                      y = 0*y,
+                                      yend = y),
+                              size = 4,
+                              color = cmmid_color)
+    } else {
+        q <- qgamma(p = c(0.0001, 0.9999),
+                    shape = R$R0_parms$shape,
+                    scale = R$R0_parms$scale)
+        
+        r0_plot <- ggplot2::ggplot(data = data.frame(R = q),
+                                   aes(x = R)) +
+            ggplot2::stat_function(geom = "area",
+                                   n = 301,
+                                   fill = cmmid_color,
+                                   fun = dgamma, 
+                                   args = list(shape = R$R0_parms$shape,
+                                               scale = R$R0_parms$scale)) 
+    }
+    
+    r0_plot +
+        ggplot2::xlab(expression(R[0])) +
+        ggplot2::ylab("Density") + 
+        ggplot2::theme_bw() +
+        ggplot2::ggtitle("Basic reproduction number") +
+        large_txt 
+    
+}
+
+plot_secondary <- function(R, dispersion){
+    
+    df <- data.frame(R = R$R0)
+    df$y <- rnbinom(n = nrow(df), size = dispersion, mu = df$R)
+    q <- qnbinom(p = 0.99, size = dispersion, mu = df$R)
+    
+    ggplot2::ggplot(data = df, aes(x=y)) +
+        ggplot2::geom_bar(width = 0.8, fill = cmmid_color, color = NA,
+                          aes(y = ..prop..)) +
+        ggplot2::theme_bw() +
+        ggplot2::xlab("Number of cases") +
+        ggplot2::ylab("Density") +
+        ggplot2::scale_x_continuous(limits = c(NA, q)) +
+        ggplot2::ggtitle("Secondary cases generated") +
+        large_txt 
+}
diff --git a/app/R/plot_data.R b/app/R/plot_data.R
index 9fe7b0f..07b407c 100644
--- a/app/R/plot_data.R
+++ b/app/R/plot_data.R
@@ -8,11 +8,9 @@
 #' @author Sam Clifford
 #' 
 
-plot_data <- function(data, reporting = 100, simulation_duration = 1,
+plot_data <- function(data, reporting = 100, 
                       title = NULL) {
     data$date <- as.Date(data$date)
-    data <- rbind(data, data.frame(date = tail(data$date, 1) + simulation_duration,
-                                   n_admissions = 0))
     data$Status <- "Reported"
     
     if (reporting < 100){
@@ -24,20 +22,35 @@ plot_data <- function(data, reporting = 100, simulation_duration = 1,
         data$Status <- factor(data$Status, levels = c("Unreported", "Reported"))
     }
     
-    ggplot2::ggplot(data = data,
+    my_palette <- my_palette[names(my_palette) %in% unique(data$Status)]
+    
+    out_plot <- ggplot2::ggplot(data = data,
                     ggplot2::aes(x = date, y = n_admissions)) +
-        ggplot2::geom_col(fill = cmmid_color, width = 0.8,
-                          aes(alpha = Status)) +
+        ggplot2::geom_col(width = 0.8,
+                          aes(fill = Status)) +
         ggplot2::xlab("Date") +
-        ggplot2::scale_alpha_manual(values = c("Unreported" = 0.5,
-                                               "Reported" = 1),
-                                    name = "Reporting status") +
+        ggplot2::scale_fill_manual(values = my_palette,
+                                   name = "Reporting status",
+                                   limits = names(my_palette)) +
         ggplot2::ylab("Number of admissions") +
         ggplot2::ggtitle(title) +
         ggplot2::theme_bw() +
         ggplot2::scale_x_date(date_label = "%d %b %y") +
         ggplot2::scale_y_continuous(breaks = int_breaks, limits = c(0, NA)) + 
         rotate_x +
-        large_txt +
-        ggplot2::theme(legend.position = "bottom")
+        large_txt 
+    
+    if (length(unique(data$Status)) > 1){
+        out_plot <- out_plot + ggplot2::theme(legend.position = "bottom")
+    } else {
+        out_plot <- out_plot + ggplot2::theme(legend.position = "none")
+    }
+    
+    if (max(data$n_admissions) == 1){
+        out_plot <- out_plot + 
+            ggplot2::scale_y_continuous(breaks = c(0, 1),
+                                        limits = c(0, NA))  
+    }
+    
+    out_plot
 }
diff --git a/app/R/plot_doubling_distribution.R b/app/R/plot_doubling_distribution.R
index c6c7a2a..3c75ea8 100644
--- a/app/R/plot_doubling_distribution.R
+++ b/app/R/plot_doubling_distribution.R
@@ -8,16 +8,21 @@
 #' @author Sam Clifford
 #' 
 
-plot_doubling_distribution <- function(x, title = NULL) {
-  dat <- data.frame(x = x)  
+plot_doubling_distribution <- function(v, trim = TRUE, ...) {
+  
+  if (trim){
+    q <- quantile(v, c(0.001, 0.999))
+    v <- v[v < q[2] & v > q[1]]
+  }
+  
+  dat <- data.frame(x = v)  
   ggplot2::ggplot(data = dat,
                   ggplot2::aes(x = x)) +
     ggplot2::geom_histogram(fill = cmmid_color,
                             col = "white",
-                            bins = 30) +
-    ggplot2::xlab("Doubling time (days)") +
-    ggplot2::ylab("Frequency") +
-    ggplot2::ggtitle(title) +
+                            bins = 30,
+                            aes(y = ..density..)) +
+    ggplot2::labs(...) +
     ggplot2::theme_bw() +
     large_txt
 }
diff --git a/app/R/plot_los_distribution.R b/app/R/plot_los_distribution.R
index bab29ce..0a86b5b 100644
--- a/app/R/plot_los_distribution.R
+++ b/app/R/plot_los_distribution.R
@@ -8,7 +8,7 @@
 #' @author Sam Clifford
 #' 
 
-plot_los_distribution <- function(los, title = NULL) {
+plot_los_distribution <- function(los, ...) {
   min_days <- max(1, los$q(0.001) - 1)
   max_days <- max(1, los$q(.999) + 2)
   days     <- min_days:max_days
@@ -18,10 +18,8 @@ plot_los_distribution <- function(los, title = NULL) {
   ggplot2::ggplot(data = dat,
                   ggplot2::aes(x = days, y = y)) +
     ggplot2::geom_col(fill = cmmid_color, width = 0.8) +
-      ggplot2::xlab("Days in hospital") +
-      ggplot2::ylab("Probability") +
-      ggplot2::ggtitle(title) +
-      ggplot2::theme_bw() +
-      ggplot2::scale_x_continuous(breaks = int_breaks) +
-  large_txt 
+    ggplot2::labs(...) +
+    ggplot2::theme_bw() +
+    ggplot2::scale_x_continuous(breaks = int_breaks) +
+    large_txt 
 }
diff --git a/app/R/plot_r_distribution.R b/app/R/plot_r_distribution.R
index d65ac34..d1c8a0d 100644
--- a/app/R/plot_r_distribution.R
+++ b/app/R/plot_r_distribution.R
@@ -19,7 +19,7 @@ plot_r_distribution <- function(los, title = NULL) {
                   ggplot2::aes(x = days, y = y)) +
     ggplot2::geom_col(fill = cmmid_color, width = 0.8) +
       ggplot2::xlab("Days in hospital") +
-      ggplot2::ylab("Probability") +
+      ggplot2::ylab("Density") +
       ggplot2::ggtitle(title) +
       ggplot2::theme_bw() +
       ggplot2::scale_x_continuous(breaks = int_breaks) +
diff --git a/app/R/plot_results.R b/app/R/plot_results.R
new file mode 100644
index 0000000..30ae4c4
--- /dev/null
+++ b/app/R/plot_results.R
@@ -0,0 +1,89 @@
+#' Function to plot outputs of predict_beds
+#'
+#' This function provides custom plots for the output of `predict_beds`.
+#' 
+#' It exists only to wrap `plot_beds` and `plot_admissions`
+#'
+#' @param x the output of `predict_beds`
+#'
+#' 
+#' @author Sam Clifford
+#' 
+#' @examples
+
+plot_results <- function(results,
+                         reporting = 100){
+    
+   
+    n_data <- nrow(results$data)
+    
+    beds  <- summarise_beds(results$beds)
+    beds$var <- "Bed occupancy"
+    beds$Status <- "Projected"
+    
+    results$data$Status <- "Reported"
+    
+    if (reporting < 100){
+        unreported <- results$data
+        unreported$Status <- "Unreported"
+        unreported$n_admissions <- round(results$data$n_admissions*100/reporting - results$data$n_admissions)
+        
+        results$data <- rbind(results$data, unreported)
+        results$data$Status <- factor(results$data$Status, levels = c("Unreported", "Reported"))
+    }
+    
+    results$data$var <- "Cases per day"
+    results$data$Date <- as.Date(results$data$date)
+    
+    cases        <- summarise_beds(results$admissions[-n_data,])
+    cases$var    <- "Cases per day"
+    cases$Status <- "Projected"
+    
+    palette_to_use <- c("Reported", 
+                        "Unreported",
+                        "Projected")[
+                            c(TRUE,
+                              reporting < 100,
+                              TRUE)]
+    
+    
+    ggplot2::ggplot(
+        mapping = aes(x = Date)) +
+        ggplot2::geom_col(data = cases,
+                          aes(fill = Status,
+                              y = Median),
+                          width = 0.8) + 
+        ggplot2::geom_col(data = results$data,
+                          aes(fill = Status,
+                              y = n_admissions),
+                          width = 0.8) + 
+        ggplot2::geom_linerange(data = cases,
+                                aes(ymin = `lower 95%`,
+                                    ymax = `upper 95%`)) + 
+        ggplot2::geom_ribbon(data = beds,
+                             aes(ymin = `lower 95%`,
+                                 ymax = `upper 95%`,
+                                 fill = Status),
+                             color = NA,
+                             alpha = 0.25) + 
+        ggplot2::geom_line(data = beds, aes(y = Median,
+                                            color = Status)) +
+        ggplot2::facet_wrap( ~ var, ncol = 1, scales = "free_y") +
+        ggplot2::theme_bw() + 
+        ggplot2::scale_fill_manual(values = my_palette[palette_to_use],
+                                   breaks = palette_to_use,
+                                   name = "Reporting status") +
+        ggplot2::scale_color_manual(values = my_palette[palette_to_use],
+                                    breaks = palette_to_use,
+                                    name = "Reporting status", 
+                                    guide = FALSE) +
+        ggplot2::theme(legend.position = "bottom",
+                       axis.title = element_blank(),
+                       strip.background = element_blank(),
+                       panel.border = element_rect(colour = "black"),
+                       strip.text = element_text(size = 18)) +
+        ggplot2::scale_y_continuous(limits = c(0, NA), breaks = int_breaks) +
+        ggplot2::scale_x_date(date_label = "%d %b %y") +
+        large_txt + rotate_x 
+   
+}
\ No newline at end of file
diff --git a/app/R/predict_admissions.R b/app/R/predict_admissions.R
index 34a28d7..cf58626 100644
--- a/app/R/predict_admissions.R
+++ b/app/R/predict_admissions.R
@@ -29,50 +29,87 @@
 #' x
 #' 
 
-predict_admissions <- function(date_start,
-                               n_start,
+predict_admissions <- function(dates,
+                               n_admissions,
                                doubling,
+                               R,
+                               si,
+                               dispersion,
                                duration,
                                reporting = 1) {
-
+  
   ## Sanity checks
-  if (length(date_start) != 1L) stop("`date_start` must contain exactly one number")
-  if (length(n_start) != 1L) stop("`n_start` must contain exactly one number")
-  if (!all(is.finite(n_start))) stop("`n_start` is not a number")
-  if (any(n_start < 1)) stop("`n_start` must be >= 1")
-
-  if (!all(is.finite(doubling))) stop("`doubling` is not a number")
-
+  if (length(dates) < 1L) stop("`dates` must contain at least one number")
+  if (length(n_admissions) < 1L) stop("`n_admissions` must contain at least one number")
+  if (!all(is.finite(n_admissions))) stop("`n_admissions` are not all numeric and finite")
+  #if (any(n_start < 1)) stop("`n_admissions` must be >= 1") # i'm not sure this is necessary
+  
+  if (!is.null(doubling) & !all(is.finite(doubling))) stop("`doubling` is not a number")
+  
   if (!is.finite(duration)) stop("`duration` is not a number")
   if (duration < 1) stop("`duration` must be >= 1")
-
+  
   if (length(reporting) != 1L) stop("`reporting` must contain exactly one value")
   if (!is.finite(reporting)) stop("`reporting` is not a number")
   if (reporting <= 0) stop("`reporting` must be > 0")
   if (reporting > 1) stop("`reporting` must be <= 1")
   
-
   ## Outline:
-
+  
   ## This function calculates future admissions using an exponential model. The
   ## growth rate is calculated from the doubling time, using: r = log(2) / d
   
   ## future dates and initial conditions
-  future_dates <- seq(date_start, length.out = duration, by = 1L)
-  initial_admissions <- round(n_start / reporting)
-
-  ## calculate growth rate from doubling times
-  r_values <- log(2) / doubling
+  future_dates       <- seq(tail(dates,1), length.out = duration, by = 1L)
+  all_admissions     <- round(n_admissions / reporting)
+  initial_admissions <- tail(all_admissions, 1) # this is for the doubling process
 
-  ## calculate future admissions
-  future_admissions <- lapply(r_values,
-                              function(r)
-    round(initial_admissions * exp(r * (seq_len(duration) - 1))))
-
-  ## build output
-  future_admissions <- matrix(unlist(future_admissions), ncol = length(doubling))
-  out <- projections::build_projections(x = future_admissions,
-                                        date = future_dates)
+  
+  if (!is.null(doubling) & length(doubling) > 0){
+    ## calculate growth rate from doubling times
+    r_values <- log(2) / doubling
+    dates_num <- as.numeric(dates)
+    tail_date <- tail(dates_num,1)
+    
+    ## calculate future admissions
+    future_admissions <- lapply(r_values,
+                                function(r){
+                                  # fit a model for the doubling/halving rate to recent data
+                                  # use this to get expected number of cases on final day of admissions
+                                  model <- stats::glm(all_admissions ~ 1, 
+                                                      offset = r*dates_num,
+                                                      family = "poisson")
+                                  admissions0 <- as.numeric(
+                                    stats::predict.glm(object = model,
+                                                       newdata = data.frame(
+                                                         dates_num = tail_date,
+                                                         r = r), 
+                                                       type = "response"))
+                                  # end fit a model
+                                  round(admissions0 * exp(r * (seq_len(duration) - 1)))
+                                })
+    
+    ## build output
+    future_admissions <- matrix(unlist(future_admissions), ncol = length(doubling))
+    out <- projections::build_projections(x = future_admissions,
+                                          date = future_dates)
+  } else {
+    
+    # use branching process
+    current_incidence <- incidence::incidence(dates = rep(dates, n_admissions))
+    out <- projections::project(x = current_incidence, 
+                                R = R,
+                                si = si,
+                                n_sim = length(R), 
+                                n_days = duration,
+                                R_fix_within = TRUE,
+                                model = "negbin",
+                                size = dispersion)
+  } 
+  
+  
+  
+  
   out
 }
 
diff --git a/app/R/predict_beds.R b/app/R/predict_beds.R
index a2e9c9c..c945b76 100644
--- a/app/R/predict_beds.R
+++ b/app/R/predict_beds.R
@@ -29,7 +29,7 @@
 ##                          doubling = c(3.5,5, 6.1, 4.6),
 ##                          duration = 14) 
 ##  x
- 
+
 ##  ## make toy duration of hospitalisation (exponential distribution)
 ##  r_duration <- function(n = 1) rexp(n, .2)
 
@@ -43,24 +43,31 @@
 ## plot(beds)
 
 predict_beds <- function(n_admissions, dates, r_los, n_sim = 10) {
-
+  
   ## sanity checks
   if (!length(dates)) stop("`dates` is empty")
-
+  
   if (!is.finite(n_admissions[1])) stop("`n_admissions` is not a number")
-  if (n_admissions[1] < 1) stop("`n_admissions` must be >= 1")
 
   if (inherits(r_los, "distcrete")) {
     r_los <- r_los$r
   }
   if (!is.function(r_los)) stop("`r_los` must be a function")
-
+  
   if (!is.finite(n_sim)) stop("`n_sim` is not a number")
   if (n_sim[1] < 1) stop("`n_sim` must be >= 1")
-
+  
+  # check if we have no new admissions
+  # if we do, return 0 bed usage
+  if (all(n_admissions < 1)) {
+    empty_proj <- projections::build_projections(x = rep(0, length(dates)), dates = dates)
+    return(projections::merge_projections(lapply(X = 1:n_sim, FUN = function(x){empty_proj})))
+  }
+  
+  
   
   ## Outline:
-
+  
   ## We take a vector of dates and incidence of admissions, and turn this into a
   ## vector of admission dates, whose length is sum(n_admissions). We will
   ## simulate for each date of admission a duration of stay, and a corresponding
@@ -68,14 +75,13 @@ predict_beds <- function(n_admissions, dates, r_los, n_sim = 10) {
   ## counted (summing up all cases) for each day. To account for stochasticity
   ## in duration of stay, this process can be replicated `n_sim` times,
   ## resulting in `n_sim` predictions of bed needs over time.
-
+  
   
   admission_dates <- rep(dates, n_admissions)
   n <- length(admission_dates)
   last_date <- max(dates)
   out <- vector(n_sim, mode = "list")
   
-
   for (j in seq_len(n_sim)) {
     los <- r_los(n)
     list_dates_beds <- lapply(seq_len(n),
@@ -83,18 +89,34 @@ predict_beds <- function(n_admissions, dates, r_los, n_sim = 10) {
                                               length.out = los[i],
                                               by = 1L))
     ## Note: unlist() doesn't work with Date objects
+    
+    # what to do when date has length 0?
+    
     dates_beds <- do.call(c, list_dates_beds)
-    beds_days <- incidence::incidence(dates_beds)
-    if (!is.null(last_date)) {
-      to_keep <- incidence::get_dates(beds_days) <= last_date
-      beds_days <- beds_days[to_keep, ]
+    
+    if (length(dates_beds) == 0){
+      out[[j]] <- projections::build_projections(x = rep(0, length(dates)), dates = dates)
+    } else {
+      
+      beds_days <- incidence::incidence(dates_beds)
+      if (!is.null(last_date)) {
+        to_keep <- incidence::get_dates(beds_days) <= last_date
+        beds_days <- beds_days[to_keep, ]
+      }
+      
+      get_beds_days <- incidence::get_dates(beds_days)
+      
+      if (!is.null(get_beds_days)){
+        
+        out[[j]] <- projections::build_projections(
+          x = beds_days$counts,
+          dates = get_beds_days)
+      } else {
+        out[[j]] <- projections::build_projections(x = rep(0, length(dates)), dates = dates)
+      }
     }
-
-    out[[j]] <- projections::build_projections(
-                                 x = beds_days$counts,
-                                 dates = incidence::get_dates(beds_days))
   }
-
+  
   projections::merge_projections(out)
- 
+  
 }
diff --git a/app/R/q_doubling.R b/app/R/q_doubling.R
index 225af1f..caca40a 100644
--- a/app/R/q_doubling.R
+++ b/app/R/q_doubling.R
@@ -15,7 +15,7 @@ q_doubling <- function(mean, cv, p=c(0.025, 0.975)) {
         q <- invgamma::qinvgamma(p=p, shape = params$shape, rate = params$rate)    
     }
     
-    short_name <- "Inv-&Gamma;"
+    short_name <- "&Gamma;<sup>-1</sup>"
     params_names <- c("&alpha;", "&beta;")
     
     return(list(short_name = short_name,
diff --git a/app/R/q_los.R b/app/R/q_los.R
index d9e42e2..be3fa6e 100644
--- a/app/R/q_los.R
+++ b/app/R/q_los.R
@@ -5,26 +5,22 @@
 #' @author Samuel Clifford
 #'
 
-q_los <- function(distribution, mean, cv, p = c(0.025, 0.975)) {
+q_los <- function(distribution, params, p = c(0.025, 0.975)) {
+    
     
-    if (cv == 0){return(rep(x = mean, times = length(p)))} else{
         
         if (distribution == "gamma"){
-            params <- epitrix::gamma_mucv2shapescale(mean, cv)
             q <- stats::qgamma(p = p, shape =  params$shape, scale = params$scale)
             short_name <- "&Gamma;"
             params_names <- c("k", "&theta;")
         }
         
         if (distribution == "weibull"){
-            params <- weibull_mucv2shapescale(mean, cv)
             q <- stats::qweibull(p = p, shape = params$shape, scale = params$scale)
             short_name <- "W"
             params_names <- c("k", "&lambda;")
         }
         
-    }
-    
     return(list(short_name = short_name,
                 q          = q,
                 params     = params,
diff --git a/app/R/q_r0.R b/app/R/q_r0.R
new file mode 100644
index 0000000..f578ae7
--- /dev/null
+++ b/app/R/q_r0.R
@@ -0,0 +1,23 @@
+#' Generate quantiles of basic reproduction number
+#'
+#' @author Samuel Clifford
+#'
+
+q_r0 <- function(mean, cv, probs = c(0.025, 0.975)) {
+    
+    if (cv == 0){return(rep(x = mean, times = length(p)))} else{
+        
+            params <- epitrix::gamma_mucv2shapescale(mean, cv)
+            params[[2]] <- 1/params[[2]]
+            names(params)[2] <- "rate"
+            q <- stats::qgamma(p = probs, shape =  params$shape, rate = params$rate)
+            short_name <- "&Gamma;"
+            params_names <- c("&alpha;", "&beta;")
+        
+    }
+    
+    return(list(short_name = short_name,
+                q          = q,
+                params     = params,
+                params_names = params_names))
+}
diff --git a/app/R/q_secondary.R b/app/R/q_secondary.R
new file mode 100644
index 0000000..4eb8943
--- /dev/null
+++ b/app/R/q_secondary.R
@@ -0,0 +1,17 @@
+q_secondary <- function(R, dispersion, probs = c(0.025, 0.5, 0.975)){
+    df <- data.frame(R = R$R0)
+    df$y <- rnbinom(n = nrow(df), size = dispersion, mu = df$R)
+    
+    q <- stats::quantile(x = df$y, probs = probs)
+    
+    nz <- mean(df$y > 0)
+    params <- fitdistrplus::fitdist(df$y, "nbinom")$estimate[c(2,1)]
+    short_name <- "NegBin"
+    params_names <- c("&mu;", "k")
+    
+    return(list(short_name = short_name,
+                q          = q,
+                params     = params,
+                params_names = params_names))
+    
+}
\ No newline at end of file
diff --git a/app/R/r_doubling.R b/app/R/r_doubling.R
index 2b6b912..524b171 100644
--- a/app/R/r_doubling.R
+++ b/app/R/r_doubling.R
@@ -6,8 +6,9 @@
 #'
 
 r_doubling <- function(n, mean, cv) {
-    ## parameterised in terms of an inverse gamma to avoid truncating at 0
-    
+  ## parameterised in terms of an inverse gamma to avoid truncating at 0
+  
+  # convert the doubling time to an inverse gamma distribution
   if (cv == 0){
     return(rep(x = mean, times = n))
   } else {
@@ -15,6 +16,5 @@ r_doubling <- function(n, mean, cv) {
       n = n,
       shape = 2 + 1/cv^2,
       rate = mean*(1 + 1/cv^2))    
-    }
-    
+  }
 }
diff --git a/app/R/run_model.R b/app/R/run_model.R
index e672610..7581625 100644
--- a/app/R/run_model.R
+++ b/app/R/run_model.R
@@ -46,12 +46,20 @@
 
 run_model <- function(dates,
                       admissions,
-                      doubling,
+                      doubling = NULL,
+                      R = NULL,
+                      si = NULL, # a distcrete object
+                      dispersion = 1,
                       duration,
                       r_los,
                       reporting = 1,
                       n_sim = 1) {
   ## check input
+  if (all(is.null(c(doubling, si, R)))){
+    msg <- "Must define either doubling times, `doubling`, or basic reproduction number, `R`, and serial interval, `si`"
+    stop(msg)
+  }
+  
   n <- length(dates)
   if (n != length(admissions)) {
     msg <- "`dates` and `admissions` have different length"
@@ -80,18 +88,24 @@ run_model <- function(dates,
     msg <- "some `dates` are missing / invalid"
     stop(msg)
   }
+  
+  out <- list(data = data.frame(date = dates,
+                                n_admissions = admissions))
 
   ord <- order(dates)
   dates <- dates[ord]
   admissions <- admissions[ord]
-  last_date <- dates[n]
-  last_admissions <- admissions[n]
+  last_date <- tail(dates, 1)
+  last_admissions <- tail(admissions,1)
   
   
   ## get projected admissions from the most recent date
-  proj_admissions <- predict_admissions(date_start = last_date,
-                                        n_start = last_admissions,
+  proj_admissions <- predict_admissions(dates = dates,
+                                        n_admissions = admissions,
                                         doubling = doubling,
+                                        R = R,
+                                        si = si,
+                                        dispersion = dispersion,
                                         duration = duration,
                                         reporting = reporting)
 
@@ -116,5 +130,8 @@ run_model <- function(dates,
                                           n_sim = n_sim))
 
   beds <- projections::merge_projections(beds)
-  beds
+  out$beds <- beds
+  out$admissions <- proj_admissions
+  return(out)
+              
 }
diff --git a/app/R/summarise_beds.R b/app/R/summarise_beds.R
index 161eee5..e5d929b 100644
--- a/app/R/summarise_beds.R
+++ b/app/R/summarise_beds.R
@@ -7,6 +7,7 @@ summarise_beds <- function(x) {
   summary_function <- function(x) {
     c(q_025 = round(quantile(x, 0.025)),
       q_25 = round(quantile(x, 0.25)),
+      q_50 = round(quantile(x, 0.5)),
       q_75 = round(quantile(x, 0.75)),
       q_95 = round(quantile(x, 0.975)))
   }
@@ -24,6 +25,7 @@ summarise_beds <- function(x) {
   colnames(summary_table) <- c("Date",
                                "lower 95%",
                                "lower 50%",
+                               "Median",
                                "upper 50%",
                                "upper 95%")
   summary_table
diff --git a/app/R/weibull_q2shapescale.R b/app/R/weibull_q2shapescale.R
new file mode 100644
index 0000000..e5d2d62
--- /dev/null
+++ b/app/R/weibull_q2shapescale.R
@@ -0,0 +1,35 @@
+weibull_q2shapescale <- function(q, p=c(0.025,0.975)){
+# https://www.johndcook.com/quantiles_parameters.pdf    
+    if (length(p) != length(q)){
+        stop("`p` and `q` must have identical lengths")
+    }
+    
+    if (any(q <= 0)){
+        stop("Weibull distribution is strictly non-negative; check quantile probabilities, `q`.")
+    }
+    
+    if (any(duplicated(p))){
+        stop("Quantile probabilities, `p`, must be unique")
+    }
+    
+    if (any(duplicated(q))){
+        stop("Quantile values, `q`, must be unique")
+    }
+    
+    if (any(p <= 0 | p >= 1)){
+        stop("Quantile probabilities, `p` must be in range (0,1)")
+    }
+    
+    if (!all(sign(diff(p[order(p)])) == sign(diff(q[order(p)])))){
+        stop("Quantile pairs do not mutually monotonically increase.")
+    }
+    
+    if (length(p) > 2){
+        warning("Only two quantiles required to specify a two-parameter distribution, using the first two elements of `p` and `q`")
+    }
+    
+    shape <- as.numeric((log(-log(1 - p[2])) - log(-log(1 - p[1]))  )/(log(q[2]) - log(q[1])))
+    scale <- q[1]/(-log(1 - p[1]))^(1/shape)
+    
+    return(list(shape = shape, scale = scale))
+}
diff --git a/app/app.R b/app/app.R
index 5e1b18f..9e30ffe 100644
--- a/app/app.R
+++ b/app/app.R
@@ -23,7 +23,8 @@ library(ggplot2)
 library(invgamma)
 library(markdown)
 library(linelist)
-
+library(shinyWidgets)
+library(fitdistrplus)
 
 ## global variables
 app_title   <- "Hospital bed occupancy projections"
@@ -49,6 +50,7 @@ ui <- navbarPage(
   ## MAIN SIMULATOR PANEL
   tabPanel(
     "Simulator",
+    shinyWidgets::chooseSliderSkin(skin = "Flat", color = "#00AEC7"),
     sidebarLayout(
       position = "left",
       
@@ -120,13 +122,12 @@ ui <- navbarPage(
           condition = sprintf("input.outputPanels == 'Length of Stay'"),
           #"Length of stay in hospital",
           h4("Description", style = sprintf("color:%s", cmmid_color)),
-          p("Parameter inputs specifying the distribution of the length of hospital stay (LoS) for COVID-19 patients. See the 'Inputs' tab for details on these distributions.",
+          p("Specifying the distribution of the length of hospital stay (LoS) for COVID-19 patients by the inter-quartile range (IQR, 25%-75% range). See the 'Inputs' tab for details on these distributions.",
             style = sprintf("color:%s", annot_color)),
           selectInput(
             "los",
             "Length of hospital stay (LoS) distribution",
-            choices = unique(los_parameters$name),
-            selected = "Custom"
+            choices = unique(los_parameters$name)
           ),
           ## Custom LoS distribution
           ## Discretised Gamma param as mean and cv
@@ -136,43 +137,92 @@ ui <- navbarPage(
                        choices = unique(los_parameters$los_dist),
                        selected = "gamma"),
           sliderInput(
-            "mean_los",
-            "Average LoS (in days)",
-            min = 1.1,
+            "los_quantiles",
+            "LoS IQR (in days)",
+            min = 1,
             max = 20,
-            value = 7,
-            step = .1),
-          sliderInput(
-            "cv_los",
-            HTML("Uncertainty as fraction of avg. (<i>c<sub>v</sub></i>)"),
-            min = 0.01,
-            max = 2,
-            value = 0.1,
-            step = .01)),
+            value = c(7, 14),
+            step = .1)),
         
         
         
         ## Epidemic growth inputs
         conditionalPanel(
-          condition = "input.outputPanels == 'Doubling time'",
+          condition = "input.outputPanels == 'Epidemic Parameters'",
           h4("Description", style = sprintf("color:%s", cmmid_color)),
-          p("Parameter inputs specifying the COVID-19 epidemic growth as doubling time and associated uncertainty. See the 'Inputs' tab for details on the doubling time distribution.",
+          p("Parameter inputs specifying the COVID-19 epidemic growth in terms of basic reproduction number and serial interval (and associated uncertainties). See the 'Inputs' tab for details on the serial interval distribution.",
             style = sprintf("color:%s", annot_color)),
-          sliderInput(
-            "doubling_time",
-            "Average doubling time (days):",
-            min = 1,
-            max = 20,
-            value = 7, 
-            step = 0.1
-          ),
-          sliderInput(
-            "uncertainty_doubling_time",
-            HTML("Uncertainty as fraction of avg. (<i>c<sub>v</sub></i>)"),
-            min = 0,
-            max = 0.5,
-            value = 0.1,
-            step = 0.01
+          radioButtons("specifyepi", label = "How do you wish to specify the epidemic growth?",
+                       choices = c("Branching process",
+                                   "Doubling/halving time"),
+                       selected = "Doubling/halving time"),
+          conditionalPanel(
+            condition = "input.specifyepi == 'Branching process'",
+            sliderInput(
+              inputId = "r0",
+              label = HTML("Average basic reproduction number, <i>R</i><sub>0</sub>"),
+              min=0.1, max=5, step=0.1, 
+              value=2.5),  
+            sliderInput(
+              "uncertainty_r0",
+              HTML("Uncertainty as fraction of avg. <i>R</i><sub>0</sub> (<i>c<sub>v,R<sub>0</sub></sub></i>)"),
+              min = 0,
+              max = 0.5,
+              value = 0.26,
+              step = 0.01
+            ),
+            # sliderTextInput(
+            #   "dispersion",
+            #   HTML("Dispersion of <i>R</i><sub>0</sub>"),
+            #   choices = c("0.1", "0.54"),
+            #   selected = "0.54"
+            # ),
+            radioButtons(inputId = "dispersion",
+                         label = HTML("Dispersion of <i>R</i><sub>0</sub>"),
+              choiceNames = c("0.1 (Endo et al.)",
+                              "0.54 (Riou and Althaus)"),
+              choiceValues = c(0.1, 0.54),
+              selected = 0.54
+            ),
+            sliderInput(
+              "serial_interval",
+              "Average serial interval (days):",
+              min = 1,
+              max = 20,
+              value = 7.5, 
+              step = 0.1
+            ),
+            sliderInput(
+              "uncertainty_serial_interval",
+              HTML("Uncertainty as fraction of avg. serial interval (<i>c<sub>v,S</sub></i>)"),
+              min = 0,
+              max = 2,
+              value = 0.45,
+              step = 0.01
+            )),
+          conditionalPanel(
+            condition = "input.specifyepi == 'Doubling/halving time'",
+            radioButtons(inputId = "doublehalf",
+              label = "The number of cases is:",
+              choiceNames = c("Doubling", "Halving"),
+              choiceValues = c(1, -1)
+            ),
+            sliderInput(
+              "doubling_time",
+              "Average time (days):",
+              min = 1,
+              max = 20,
+              value = 7.7, 
+              step = 0.1
+            ),
+            sliderInput(
+              "uncertainty_doubling_time",
+              HTML("Uncertainty as fraction of avg. time (<i>c<sub>v,T</sub></i>)"),
+              min = 0,
+              max = 1,
+              value = 0.33,
+              step = 0.01
+            )
           )
         ),
         
@@ -180,6 +230,8 @@ ui <- navbarPage(
         conditionalPanel(
           condition = "input.outputPanels == 'Main results'",
           h4("Description", style = sprintf("color:%s", cmmid_color)),
+          p("Median and 95% intervals for projected bed occupancy and daily number of new cases.",
+            style = sprintf("color:%s", annot_color)),
           p("Parameter inputs specifying the number and durations of the simulations.",
             style = sprintf("color:%s", annot_color)),
           sliderInput(
@@ -219,24 +271,46 @@ ui <- navbarPage(
             "Length of Stay",
             
             br(),
-            plotOutput("los_plot", width = "30%", height = "300px"),
+            plotOutput("los_plot", width = "100%", height = "300px"),
             br(),
             htmlOutput("los_ci")
           ),
           tabPanel(
-            "Doubling time",
+            "Epidemic Parameters",
             
             br(),
-            plotOutput("doubling_plot", width = "30%", height = "300px"),
-            br(),
-            htmlOutput("doubling_ci")
+            
+            conditionalPanel(condition = "input.specifyepi == 'Doubling/halving time'",
+                             plotOutput("doubling_plot", width = "30%", height = "300px"),
+                             br(),
+                             htmlOutput("doubling_ci")
+            ),
+            # need to amend this panel to have three fluidrows and two columns each,
+            # plot in left col, summary on right
+            conditionalPanel(condition = "input.specifyepi == 'Branching process'",
+                             
+                             fluidRow(
+                               column(6, plotOutput("r0_plot", height = "200px")),
+                               column(6, br(), br(), htmlOutput("r0_ci"))),
+                             
+                             fluidRow(
+                               column(6, plotOutput("secondary_plot", height = "200px")),
+                               column(6, br(), br(), htmlOutput("secondary_ci"))),
+                             br(),
+                             fluidRow(
+                               column(6, plotOutput("serial_plot", height = "200px")),
+                               column(6, br(), br(), htmlOutput("serial_ci")))
+            )
+            
+            
+            
           ),
           tabPanel(
             "Main results",
-          
+            
             br(),
-            plotOutput("main_plot", width = "30%", height = "300px"),
-            checkboxInput("show_table", "Show summary table?", FALSE),
+            plotOutput("main_plot", width = "30%", height = "600px"),
+            checkboxInput("show_table", "Show bed occupancy summary table?", FALSE),
             conditionalPanel(
               condition = sprintf("input.show_table == true"),
               DT::dataTableOutput("main_table", width = "50%"))
@@ -267,6 +341,7 @@ ui <- navbarPage(
 
 
 
+
 #################
 ## SERVER SIDE ##
 #################
@@ -282,12 +357,9 @@ server <- function(input, output, session) {
   observe({
     default=input$los
     
-    updateSliderInput(session, "mean_los",
-                      value = c(los_parameters[los_parameters$name == default, "mean_los"]))
-    
-    updateSliderInput(session, "cv_los",
-                      value = c(los_parameters[los_parameters$name == default, "cv_los"]))
-    
+    updateSliderInput(session, "los_quantiles",
+                      value = as.numeric(los_parameters[los_parameters$name == default, c("los_25", "los_75")]))
+
     updateRadioButtons(session, "los_dist",
                        selected = c(los_parameters[los_parameters$name == default, "los_dist"]))
     
@@ -295,6 +367,9 @@ server <- function(input, output, session) {
   }
   
   )
+  # 
+  # observeEvent(eventExpr = data, 
+  #              {output$main_plot =renderPlot({})})
   
   
   ## GENERAL PROCESSING OF INPUTS: INTERNAL CONSTRUCTS
@@ -319,41 +394,64 @@ server <- function(input, output, session) {
   los <- reactive({
     los_dist(
       distribution          = input$los_dist,
-      mean                  = input$mean_los,
-      cv                    = input$cv_los)
+      q                     = input$los_quantiles)
   })
   
-  ## doubling time (returns a vector or r values)
+  ## doubling time (returns a vector of doubling time values)
   doubling <-  reactive({
-    r_doubling(n = input$number_simulations,
-               mean = input$doubling_time,
-               cv = input$uncertainty_doubling_time)
+    r_doubling(n       = input$number_simulations,
+               mean    = input$doubling_time,
+               cv      = input$uncertainty_doubling_time)
   })
   ## same, but larger sample to plot distribution
   doubling_large <-  reactive({
     r_doubling(n = 1e5,
-               mean = input$doubling_time,
-               cv = input$uncertainty_doubling_time)
+               mean    = input$doubling_time,
+               cv      = input$uncertainty_doubling_time)
   })
   
-  
   output$data_plot <- renderPlot({
     
     ggdata <- data()
     reporting <- input$assumed_reporting
-    simulation_duration <- input$simulation_duration
-    plot_data(data = ggdata, reporting = reporting, simulation_duration)
+    plot_data(data = ggdata, reporting = reporting)
     
   })
   
+  si <- reactive({
+    make_si(mean = input$serial_interval,
+            cv   = input$uncertainty_serial_interval)
+  })
+  
+  R <- reactive({
+    make_r0(n    = input$number_simulations,
+            mean = input$r0,
+            cv   = input$uncertainty_r0)
+  })
+  
+  R_large <-  reactive({
+    make_r0(n    = 1e5,
+            mean = input$r0,
+            cv   = input$uncertainty_r0)
+  })
+  
   ## main results
   results <- eventReactive(
     input$run,
     if (!is.null(data())) {
+      
+      if(input$specifyepi == "Doubling/halving time"){
+        doubling_ <- doubling()} else {
+          doubling_ <- NULL
+        }
+      
       run_model(
         dates = data()$date,
         admissions = data()$n_admissions,
-        doubling = doubling(),
+        doubling = as.numeric(input$doublehalf)*doubling_,
+        R = R()$R0,  
+        si = si(),
+        dispersion = as.numeric(input$dispersion),
         duration = input$simulation_duration,
         r_los = los()$r,
         reporting = input$assumed_reporting / 100,
@@ -364,30 +462,67 @@ server <- function(input, output, session) {
     ignoreNULL = FALSE
   )
   
-  
+  output$doublehalftext <- reactive({
+    if (input$doublehalf == 1){
+      "doubling"
+    } else {
+      "halving"
+    }
+  })
   
   ## PLOTS  
   ## graph for the distribution of length of hospital stay (LoS)
   output$los_plot <- renderPlot(
     plot_los_distribution(
-      los(), "Length of stay in hospital"
-    ), width = 600
+      los(), 
+      title = "Length of stay in hospital",
+      x     = "Days in hospital",
+      y     = "Density"
+    )
+  )
+  
+  output$r0_plot <- renderPlot(
+    plot_r0(
+      R_large()
+    ), width = 300, height = 200
+      #function() {
+      #0.6*session$clientData$output_r0_plot_width
+    #}
+  )
+  
+  output$secondary_plot <- renderPlot({
+
+    plot_secondary(
+      R_large(),
+      dispersion = as.numeric(input$dispersion)
+    )}, width = 300, height = 200
+  )
+  
+  output$serial_plot <- renderPlot({
+    plot_los_distribution(los = si(), 
+                          title = "Serial interval", 
+                          x = "Days to secondary case",
+                          y = "Density")
+  }, width = 300, height = 200
   )
   
   ## graph for the distribution of length of hospital stay (LoS)
-  output$doubling_plot <- renderPlot(
+  output$doubling_plot <- renderPlot({
+    
     plot_doubling_distribution(
-      doubling_large(), "Epidemic doubling time"
-    ), width = 600
+      doubling_large(),
+      title = sprintf("Epidemic %s time", dhlabel(input$doublehalf)), 
+      x = "Days", 
+      y = "Density")
+  }, width = 600
   )
   
   
   ## main plot: predictions of bed occupancy
   output$main_plot <- renderPlot({
-    plot_beds(results(),
-              ribbon_color = slider_color,
-              palette = cmmid_pal,
-              title = "Projected bed occupancy")
+    plot_results(results = results(),
+                 reporting = input$assumed_reporting)
+    
   }, width = 600)
   
   
@@ -396,50 +531,113 @@ server <- function(input, output, session) {
   
   ## summary tables
   output$main_table <- DT::renderDataTable({
-    summarise_beds(results())
+    summarise_beds(results()$beds)
   })
   
   
-  
   ## OTHERS
   
   ## confidence interval for length of stay
+  # rework this
+  los_params_values <- reactive({
+    los_params(distribution = input$los_dist,
+               q = input$los_quantiles)
+  })
+  
   output$los_ci <- reactive({
     q <- q_los(distribution = input$los_dist,
-               mean = input$mean_los, 
-               cv   = input$cv_los,
+               params = los_params_values(),
                p = c(0.025, 0.5, 0.975))
-    #sprintf("<b>LoS distribution:</b> %s(%0.1f, %0.1f)", )
-    sprintf("<b>Median LoS:</b> %0.1f<br>
-            <b>95%% interval</b>: (%0.1f, %0.1f)<br>
+    
+    if (q$q[3] != q$q[1]){
+      sprintf("<b>Median LoS:</b> %0.1f days<br>
+            <b>95%% interval</b>: (%0.1f, %0.1f) days<br>
             <b>Distribution:</b> %s(<i>%s</i>=%0.1f, <i>%s</i>=%0.1f)",
-            q$q[2], q$q[1], q$q[3], q$short_name,
-            q$params_names[1], q$params[1],
-            q$params_names[2], q$params[2])
+              q$q[2], q$q[1], q$q[3], q$short_name,
+              q$params_names[1], q$params[1],
+              q$params_names[2], q$params[2])
+    } else {
+      sprintf("<b>Fixed LoS:</b> %0.1f days<br>",
+              q$q[2])
+    }
+    
+    #sprintf("<b>LoS distribution:</b> %s(%0.1f, %0.1f)", )
+    
   })
   
   ## confidence interval for doubling time 
-  output$doubling_ci <- reactive({
+  output$serial_ci <- reactive({
     
-    if (input$uncertainty_doubling_time > 0){
+    if (input$uncertainty_serial_interval > 0){
+      
+      q <- q_doubling(mean = input$serial_interval, 
+                      cv   = input$uncertainty_serial_interval,
+                      p = c(0.025, 0.5, 0.975))
+      
+      
+      sprintf("<b>Median serial interval:</b> %0.1f days<br>
+            <b>95%% interval:</b> (%0.1f, %0.1f) days<br>
+            <b>Distribution:</b> %s(<i>%s</i>=%0.1f, <i>%s</i>=%0.1f)",
+              q$q[2], q$q[1], q$q[3], q$short_name,
+              q$params_names[1], q$params[1],
+              q$params_names[2], q$params[2])
+    } else {
+      sprintf("<b>Fixed serial interval:</b> %0.1f days<br>",
+              input$serial_interval)
+    }
+  })
+  
+  output$doubling_ci <- reactive({
     
-    q <- q_doubling(mean = input$doubling_time, 
-                    cv   = input$uncertainty_doubling_time,
-                    p = c(0.025, 0.5, 0.975))
+    q <- q_doubling(mean = input$doubling_time,
+                    cv = input$uncertainty_doubling_time, p = c(0.025, 0.5, 0.975))
     
     
-    sprintf("<b>Median doubling time:</b> %0.1f<br>
+    sprintf("<b>Median %s time:</b> %0.1f days<br>
             <b>95%% interval:</b> (%0.1f, %0.1f)<br>
             <b>Distribution:</b> %s(<i>%s</i>=%0.1f, <i>%s</i>=%0.1f)",
+            HTML(dhlabel(input$doublehalf)),
             q$q[2], q$q[1], q$q[3], q$short_name,
             q$params_names[1], q$params[1],
             q$params_names[2], q$params[2])
+  })
+  
+  output$r0_ci <- reactive({
+    
+    if (input$uncertainty_r0 > 0){
+      
+      q <- q_r0(mean = input$r0, 
+                cv   = input$uncertainty_r0,
+                p = c(0.025, 0.5, 0.975))
+      
+      sprintf("<b>Median <i>R</i><sub>0</sub>:</b> %0.1f<br>
+            <b>95%% interval:</b> (%0.1f, %0.1f)<br>
+            <b>Distribution:</b> %s(<i>%s</i>=%0.1f, <i>%s</i>=%0.1f)",
+              q$q[2], q$q[1], q$q[3], q$short_name,
+              q$params_names[1], q$params[1],
+              q$params_names[2], q$params[2])
     } else {
-      sprintf("<b>Fixed doubling time:</b> %0.1f<br>",
-              input$doubling_time)
+      sprintf("<b>Fixed <i>R</i><sub>0</sub>:</b> %0.1f<br>",
+              input$r0)
     }
   })
   
+  output$secondary_ci <- reactive({
+    
+    
+    q <- q_secondary(R_large(),
+                     as.numeric(input$dispersion),
+                     p = c(0.025, 0.5, 0.975))
+    
+    sprintf("<b>Median secondary cases:</b> %0.1f<br>
+            <b>95%% interval:</b> (%0.1f, %0.1f)<br>
+            <b>Distribution:</b> %s(<i>%s</i>=%0.1f, <i>%s</i>=%0.2f)",
+            q$q[2], q$q[1], q$q[3], q$short_name,
+            q$params_names[1], q$params[1],
+            q$params_names[2], q$params[2])
+    
+  })
+  
 }
 
 
diff --git a/app/include/heading_box.md b/app/include/heading_box.md
index 925f3bd..6d7c0cd 100644
--- a/app/include/heading_box.md
+++ b/app/include/heading_box.md
@@ -9,7 +9,7 @@
 <br>
 
 Forecast bed occupancy for COVID-19 patients using recent data on admission and
-an exponential growth model. Click on the <i>simulator</i> tab to use the app. All
+an exponential growth model. Click on the <i>Simulator</i> tab to use the app. All
 inputs are provided using the left panel. Other tabs provide further
 information:
 
diff --git a/app/include/inputs.md b/app/include/inputs.md
index abddc87..deece34 100644
--- a/app/include/inputs.md
+++ b/app/include/inputs.md
@@ -1,5 +1,5 @@
 
-## Data
+## Admissions data
 
 ### Minimum data: admissions on a single day
 
@@ -54,32 +54,58 @@ Currently available options are:
   values from the values of &mu; and &sigma;). Note that the distribution is generated so that LoS
   must be positive.
 
-* **Zhou et al non-critical care**: discretised Weibull (shape: 2, scale: 13) targeting a median of 11
+* **Rees and Nightingale et al. non-critical care**: discretised Weibull (shape: 1.2, scale: 6.9) resulting from a systematic review of length of stay for general (non-ICU) admissions outside of China ([Rees and Nightingale et al., 2020](https://doi.org/10.1101/2020.04.30.20084780)).
+
+* **Rees and Nightingale et al. critical care**: discretised Weibull (shape: 1.5, scale: 8.7) resulting from a systematic review of length of stay for ICU admissions outside of China ([Rees and Nightingale et al., 2020](https://doi.org/10.1101/2020.04.30.20084780)).
+
+* **Zhou et al. non-critical care**: discretised Weibull (shape: 2, scale: 13) targeting a median of 11
     days, IQR 7-14, as reported for general hospitalisation in <a
 	href="https://www.thelancet.com/journals/lancet/article/PIIS0140-6736(20)30566-3/fulltext">Zhou
 	et al., 2020</a>:
 
 
-* **Zhou et al critical care**: discretised Weibull (shape: 2, scale: 10) targeting a median of 8
+* **Zhou et al. critical care**: discretised Weibull (shape: 2, scale: 10) targeting a median of 8
     days, IQR 4-12, as reported for hospitalisation in critical care in <a
 	href="https://www.thelancet.com/journals/lancet/article/PIIS0140-6736(20)30566-3/fulltext">Zhou
 	et al., 2020</a>:
 
-
+The default option is the non-critical setting from non-China studies published by [Rees and Nightingale et al.  (2020)](https://doi.org/10.1101/2020.04.30.20084780)
   
-## Growth parameters
-
-* **Assumed doubling time (days)** This is the estimated (mean) time taken for the epidemic to double in size, and serves as a measure of transmission intensity.
-    + Default: 7
-    + Plausible ranges: 1.8 - 9.3. See [Muniz-Rodriguez et al 2020](https://www.medrxiv.org/content/10.1101/2020.02.05.20020750v4.full.pdf), [Zhao et al 2020](https://www.medrxiv.org/content/medrxiv/early/2020/02/29/2020.02.26.20028449.full.pdf), [Wu et al 2020](https://www.nature.com/articles/s41591-020-0822-7), [Li et al 2020](https://www.nejm.org/doi/full/10.1056/NEJMoa2001316), [Cheng et al 2020](https://link.springer.com/content/pdf/10.1007/s15010-020-01401-y.pdf) and [Granozio 2020](https://arxiv.org/ftp/arxiv/papers/2003/2003.08661.pdf) for references. 
-* **Uncertainty in doubling time (coefficient of variation)**: the sampling
-  distribution for the doubling time is an inverse gamma distribution
-  parameterised in terms of the mean doubling time (defined by the user) and the
-  coefficient of variation (i.e. &sigma;/&mu;) which are then used for moment-matching
+## Epidemic growth parameters
+
+### Doubling/halving time
+
+* **Assumed doubling/halving time (days)** This is the estimated (mean) time taken for the epidemic to double (or halve) in size, and serves as a measure of transmission intensity.
+    + Default: 7.7
+    + Plausible ranges: 1.8 - 9.3. See [Muniz-Rodriguez et al 2020](https://doi.org/10.3201/eid2608.200219), [Zhao et al 2020](https://www.medrxiv.org/content/medrxiv/early/2020/02/29/2020.02.26.20028449.full.pdf), [Wu et al 2020](https://www.nature.com/articles/s41591-020-0822-7), [Li et al 2020](https://www.nejm.org/doi/full/10.1056/NEJMoa2001316), [Cheng et al 2020](https://link.springer.com/content/pdf/10.1007/s15010-020-01401-y.pdf) and [Granozio 2020](https://arxiv.org/ftp/arxiv/papers/2003/2003.08661.pdf) for references. 
+* **Uncertainty in doubling/halving time (coefficient of variation)**: the sampling
+  distribution for the doubling/halving time is an inverse gamma distribution
+  parameterised in terms of the mean doubling/havling time (defined by the user) and the
+  coefficient of variation (i.e. <i>c<sub>v,T</sub> = &sigma;/&mu;</i>) which are then used for moment-matching
   to determine appropriate values of the distribution's shape and rate parameters (&alpha;, &beta;).
-    + Default: &sigma;/&mu; = 0.1
+    + Default: <i>c<sub>v,T</sub></i> = 0.33
+
+The default values here are drawn from the early dynamics reported by <a href="https://www.nejm.org/doi/full/10.1056/NEJMoa2001316">Li et al. (2020)</a>.
 
+When "halving" is selected, the epidemic is considered to be decreasing in intensity by exponential decay in case numbers rather than increasing when characterised by a doubling time.
 
+### Branching process
+
+Alternatively, the epidemic may be parameterised according to a branching process where every generation time (or serial interval) each case generates a number of secondary cases. The number of secondary cases is parameterised with a negative binomial distribution with a mean, <i>&mu;</i>, and dispersion, _k_, parameter.
+
+* **Basic reproduction number**, i.e <i>R</i><sub>0</sub> the average number of secondary infections each case goes on to cause. This is parameterised by an average and coefficient of variation <i>c<sub>v,R<sub>0</sub></sub></i> which are moment matched to the shape, &alpha;, and scale, &beta;, parameters of a gamma distribution. The average behaviour of the epidemic is governed by the size of <i>R</i><sub>0</sub>, namely whether or not sufficient new cases are being generated to replace those who recover.
+    - <i>R</i><sub>0</sub> &gt; 1: the epidemic grows exponentially as cases infect above replacement
+    - <i>R</i><sub>0</sub> &asymp; 1: the epidemic is stable 
+    - <i>R</i><sub>0</sub> &lt; 1: cases are not being replaced and the epidemic growth slows until there are no new cases
+    
+The default value is for <i>R</i><sub>0</sub> to have an average of 2.5 and <i>c<sub>v,R<sub>0</sub></sub></i> of 0.26, corresponding to a 95% interval of (1.4, 3.9) which matches the early dynamics of the outbreak in Wuhan ([Li et al., 2020](https://doi.org/10.1056/NEJMoa2001316)) prior to the introduction of lockdown, travel restrictions or social distancing measures. Check local estimates of reproduction in order to describe how much transmission is occurring in your setting, e.g. https://epiforecasts.io/covid/reports.html.  
+
+* **Dispersion parameter**, <i>k</i>, for a negative binomial distribution with mean <i>&mu; = R</i><sub>0</sub>, describes the number of cases <i>each individual</i> is likely to go on to infect. For smaller values of <i>k</i> and the same <i>R</i><sub>0</sub> it is more probable that an individual will cause 0 new infections. This parameter therefore controls how likely superspreading events are in driving the epidemic. The default value of 0.54 is based on an early estimate of human to human transmission ([Riou and Althaus, 2020](https://doi.org/10.2807/1560-7917.ES.2020.25.4.2000058)) and an alternative is provided based on analysis indicating that 80% of transmission is caused by 10% of cases ([Endo et al., 2020](https://doi.org/10.12688/wellcomeopenres.15842.1)).
+
+* **Serial interval**, the length of time between successive generations of infection. This describes, on average, how long it takes from the infection of a case to a secondary infection. Specified as an average time and coefficient of variation, <i>c<sub>v,S</sub></i>, the specified parameters are moment matched to the shape, &alpha;, and scale, &beta;, of an inverse gamma distribution.
+    * Defaults are based on the early transmission dynamics in [Li et al 2020](https://www.nejm.org/doi/full/10.1056/NEJMoa2001316)
+        - Average serial interval: 7.5 days
+        - <i>c<sub>v,S</sub></i>: 0.45
 
 ## Simulation parameters
 
@@ -93,6 +119,12 @@ Currently available options are:
 
 ## References
 
-Zhou, Fei, et al. "Clinical Course and Risk Factors for Mortality of Adult Inpatients with COVID-19 in Wuhan, China: a Retrospective Cohort Study." _The Lancet_, 2020. <a href="https://doi.org/10.1016/s0140-6736(20)30566-3">doi:10.1016/s0140-6736(20)30566-3</a>.
+Endo, A., et al. "Estimating the overdispersion in COVID-19 transmission using outbreak sizes outside China [version 1; peer review: 1 approved]", _Wellcome Open Research_, 2020. https://doi.org/10.12688/wellcomeopenres.15842.1
+
+Li, Q., et al. "Early Transmission Dynamics in Wuhan, China, of Novel Coronavirus–Infected Pneumonia." _NEJM_, 2020. https://doi.org/10.1056/NEJMoa2001316
+
+Rees, E. M., Nightingale, E. S., et al. "COVID-19 length of hospital stay: a systematic review and data synthesis." _medR&chi;iv_, 2020. https://doi.org/10.1101/2020.04.30.20084780
 
+Riou, J. and Althaus, C. "Pattern of early human-to-human transmission of Wuhan 2019 novel coronavirus (2019-nCoV), December 2019 to January 2020", _Euro Surveillance_, 2020. https://doi.org/10.2807/1560-7917.ES.2020.25.4.2000058
 
+Zhou, Fei, et al. "Clinical Course and Risk Factors for Mortality of Adult Inpatients with COVID-19 in Wuhan, China: a Retrospective Cohort Study." _The Lancet_, 2020. <a href="https://doi.org/10.1016/s0140-6736(20)30566-3">doi:10.1016/s0140-6736(20)30566-3</a>
diff --git a/app/include/model.md b/app/include/model.md
index 2709856..0230726 100644
--- a/app/include/model.md
+++ b/app/include/model.md
@@ -27,22 +27,23 @@ The forecasting approach can be summarised as follows; for each simulation:
    currently done by:
 <div> n<sub>aug</sub> = n<sub>reported</sub> / %<sub>reported</sub> </div>
  
-2. Use a log-linear model, parametrised via the doubling time (drawn from an
-   inverse Gamma distribution, with user-specified mean and coefficient of
-   variation), to simulate future daily admissions trajectories.
+2. Simulate future daily admissions trajectories according to one of the two approaches:
+    - Use a log-linear model, parameterised via the doubling (or halving) time for exponential growth (or decay)
+    - Use a branching process, parameterised via the basic reproduction number, its dispersion, and the serial interval for successive generations. 
 
 3. For each admission, simulate length of stay in hospital from the specified
    distribution.
 
-4. Count beds for each day simulation.
+4. Count beds for each day of the simulation.
 
 
 
 ## Caveats
 
-* The current model assumes exponential growth. This is generally a reasonable
+* The current model assumes exponential growth as its default. This is generally a reasonable
   approximation at the beginning of an epidemic however will become less
-  appropriate as growth slows towards the peak.
+  appropriate as growth slows towards the peak. At this point, it may be more appropriate
+  to switch to the branching process. 
 
 * User inputted data only contributes to bed requirements and any trend in user
   data does not impact on the assumed exponential growth from the final time
@@ -54,7 +55,7 @@ The forecasting approach can be summarised as follows; for each simulation:
   etc.
 
 
-### References
+## References
 
 Jombart et al. 2020 "Forecasting critical care bed requirements for COVID-19 patients
 in England". CMMID post, first online