diff --git a/NAMESPACE b/NAMESPACE
index 64940562cce8427561592d4fe2b800a7977ab288..6aeed64b484c17849479b797122650c148fc524b 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,27 +1,27 @@
 # Generated by roxygen2: do not edit by hand
 
 S3method(as.data.frame,seiqhrf)
+S3method(plot,list)
+S3method(plot,seiqhrf)
 S3method(print,control.seiqhrf)
 S3method(print,init.seiqhrf)
 S3method(print,param.seiqhrf)
 S3method(print,seiqhrf)
+S3method(summary,seiqhrf)
 export(arrivals.FUN)
 export(control_seiqhrf)
 export(cum_discr_si)
 export(departures.FUN)
 export(get_prev.FUN)
-export(get_times)
 export(infection.FUN)
 export(init_seiqhrf)
 export(init_status.icm)
 export(initialize.FUN)
 export(param_seiqhrf)
-export(plot_models)
 export(plot_times)
 export(recovery.FUN)
 export(seiqhrf)
 export(simulate_seiqhrf)
-export(summary_seiqhrf)
 export(vary_param)
 import(dplyr)
 import(ggplot2)
diff --git a/R/results-parse.R b/R/FIN_plot.R
similarity index 63%
rename from R/results-parse.R
rename to R/FIN_plot.R
index 83fcbfe2a853d1c6bde81e15f17146a954240976..cd1d0c447bddabc859f62c9a25cae534093672a7 100644
--- a/R/results-parse.R
+++ b/R/FIN_plot.R
@@ -1,10 +1,75 @@
-#' Get data by time
+#' Plot simulation result 
 #'
-#' Function to extract timings and assemble results in a dataframe
+#' Function to extract timings, prevalence etc. from the simulation 
 #'
-#' @param simulate_results results from `simulte()` function.
+#' @param x An seiqhrf object returned from function \code{\link{seiqhrf}}.
+#' @param method If "times", plot Duration frequency distributions.
+#'               If NULL, plot individuals models or mutliple models for comparison.
+#' @param comp_remove Compartments to remove. Suggest c(s.num, r.num)
+#' @param time_lim Number of steps (days) to plot. 
+#' @param trans Y-axis transformation (e.g. log2, log10). Default = none. 
+#' @param known Dataframe with known compartment numbers to plot alongside
+#'        projections
+#' @param start_date Date for day 0. Default: ymd("2020-03-21"),
+#' @param x_axis Title for x-axis. Default: 'Days since beginning of epidemic'
+#' @param plot_title Title for whole plot. Default: 'SEIQHRF plot'
+#' @param ... Additional parameters
+#' 
+#' @return ggplot2 object
+#' 
+#' @importFrom tidyr pivot_longer
+#' @importFrom dplyr mutate
+#' @importFrom dplyr "%>%"
+#' @importFrom dplyr bind_rows
+#' @importFrom dplyr select
+#' @importFrom dplyr filter
+#' @import ggplot2
+#' @export
+plot.seiqhrf <- function(x, method = NULL, comp_remove = "none",
+                         time_lim = 100,
+                         trans = 'na',
+                         known = NULL,
+                         start_date = ymd("2020-03-21"),
+                         x_axis = 'Days since beginning of epidemic',
+                         plot_title = 'SEIQHRF', ...) {
+    
+    if(is.null(method)){
+        
+        plot_model(x,
+                   comp_remove = comp_remove,
+                   time_lim = time_lim,
+                   trans = trans,
+                   known = known,
+                   start_date = start_date,
+                   x_axis = x_axis,
+                   plot_title = plot_title)
+    }else{
+        if(method == "times"){
+            if(!inherits(x, "seiqhrf")) stop("If method == times, x needs to be an seiqhrf object")
+            plot_times(x) 
+        } 
+    }
+    
+}
+
+
+
+#' Plot simulation result 
+#'
+#' Function to extract timings, prevalence etc. from the simulation 
+#'
+#' @param x A list of seiqhrf objects returned from function \code{\link{seiqhrf}}.
+#' @param comp_remove Compartments to remove. Suggest c(s.num, r.num)
+#' @param time_lim Number of steps (days) to plot. 
+#' @param trans Y-axis transformation (e.g. log2, log10). Default = none. 
+#' @param known Dataframe with known compartment numbers to plot alongside
+#'        projections
+#' @param start_date Date for day 0. Default: ymd("2020-03-21"),
+#' @param x_axis Title for x-axis. Default: 'Days since beginning of epidemic'
+#' @param plot_title Title for whole plot. Default: 'SEIQHRF plot'
+#' @param ... Additional parameters
 #' 
-#' @return dataframe with sim.n, period type, and duration
+#' @return ggplot2 object
 #' 
 #' @importFrom tidyr pivot_longer
 #' @importFrom dplyr mutate
@@ -12,10 +77,44 @@
 #' @importFrom dplyr bind_rows
 #' @importFrom dplyr select
 #' @importFrom dplyr filter
+#' @import ggplot2
 #' @export
-get_times <- function(simulate_results) {
+plot.list <- function(x, comp_remove = "none",
+                      time_lim = 100,
+                      trans = 'na',
+                      known = NULL,
+                      start_date = ymd("2020-03-21"),
+                      x_axis = 'Days since beginning of epidemic',
+                      plot_title = 'SEIQHRF', ...){
     
-    sim <- simulate_results$sim
+    plot_model(x,
+               comp_remove = comp_remove,
+               time_lim = time_lim,
+               trans = trans,
+               known = known,
+               start_date = start_date,
+               x_axis = x_axis,
+               plot_title = plot_title)
+}
+
+#' Plot simulation by time
+#'
+#' Function to plot Duration frequency distributions.
+#'
+#' @param sim An seiqhrf object returned from function \code{\link{seiqhrf}}.
+#' 
+#' @return ggplot2 object
+#' 
+#' @import ggplot2
+#' @importFrom tidyr pivot_longer
+#' @importFrom dplyr mutate
+#' @importFrom dplyr "%>%"
+#' @importFrom dplyr bind_rows
+#' @importFrom dplyr select
+#' @importFrom dplyr filter
+#' 
+#' @export
+plot_times <- function(sim) {
     
     for (s in 1:sim$control$nsims) {
         if (s == 1) {
@@ -52,80 +151,40 @@ get_times <- function(simulate_results) {
                                                "Hosp care required duration", 
                                                "Survival time case fatalities"),
                                     ordered = TRUE))
-
-    return(times)
-}
-
-
-#' Plot times
-#'
-#' Function to plot Duration frequency distributions.
-#'
-#' @param times results from `get_times()` function.
-#'
-#' @return ggplot2 object
-#' @import ggplot2
-#' @export
-plot_times <- function(times){
+    
     times %>% filter(duration <= 30) %>% ggplot(aes(x = duration)) + 
         geom_bar() + facet_grid(period_type ~ ., scales = "free_y") + 
         labs(title = "Duration frequency distributions",
              subtitle = "Baseline simulation")
 }
 
-#' Plot models
-#'
-#' Function to plot individuals models or mutliple models for comparison.
-#'
-#' @param sims Single or list of sims to plot
-#' @param sim_id String or list of strings to name each facet
-#' @param comp_remove Compartments to remove. Suggest c(s.num, r.num)
-#' @param time_lim Number of steps (days) to plot. 
-#' @param trans Y-axis transformation (e.g. log2, log10). Default = none. 
-#' @param known Dataframe with known compartment numbers to plot alongside
-#'        projections
-#' @param start_date Date for day 0. Default: ymd("2020-03-21"),
-#' @param x_axis Title for x-axis. Default: 'Days since beginning of epidemic'
-#' @param plot_title Title for whole plot. Default: 'ICM plot'
-#'
-#' @return ggplot2 object
-#' 
-#' @import dplyr
-#' @importFrom tidyr pivot_longer
-#' 
-#' @export
-plot_models <- function(sims = baseline_sim,
-                        sim_id = "baseline",
-                        comp_remove = "none",
-                        time_lim = 100,
-                        trans = 'na',
-                        known = NULL,
-                        start_date = ymd("2020-03-21"),
-                        x_axis = 'Days since beginning of epidemic',
-                        plot_title = 'ICM plot'){
-    
-    # Define a standard set of colours to represent compartments
-    comps <- c("s.num", "e.num", "i.num", "q.num", "h.num", "r.num", "f.num")
-    compcols <- c(s.num = "#4477AA", e.num = "#66CCEE", i.num = "#CCBB44", 
-                  q.num = "#AA3377", h.num = "#EE6677", r.num = "#228833", 
-                  f.num = "#BBBBBB")
-    complabels <- c(s.num = "S: Susceptible", e.num = "E: Asymptomatic", 
-                    i.num = "I: Infected", q.num = "Q: Self-isolated", 
-                    h.num = "H: Hospitalized", r.num = "R: Recovered",
-                    f.num = "F: Case Fatalities")
+
+
+# Internal
+plot_model <- function(sims, comp_remove, time_lim, trans, known, start_date, x_axis, plot_title){
     
     # Merge models to plot together
-    for (i in (1: length(sim_id))) {
-        if (i == 1) {
-            plot_df <- sims[i*2]$df
-            plot_df <- plot_df %>% mutate(experiment = sim_id[i])
-        } else{
-            tmp_df <- sims[i*2]$df
-            tmp_df <- tmp_df %>% mutate(experiment = sim_id[i])
-            plot_df <- plot_df %>% bind_rows(plot_df, tmp_df)
+    if(class(sims) == "seiqhrf"){
+        sim_id <- "Baseline" 
+        plot_df <- as.data.frame(sims, out = "mean")
+        plot_df <- plot_df %>% mutate(experiment = sim_id)
+    }else{
+        
+        sim_id <- names(sims)
+        if(is.null(sim_id)) stop("Please assign a name to each element in sims")
+        
+        plot_df <- as.data.frame(sims[[1]], out = "mean")
+        plot_df <- plot_df %>% mutate(experiment = sim_id[1])
+        if(length(sim_id) > 1){
+            for (i in (2:length(sim_id))) {
+                tmp_df <- as.data.frame(sims[[i]], out = "mean")
+                tmp_df <- tmp_df %>% mutate(experiment = sim_id[i])
+                plot_df <- plot_df %>% bind_rows(plot_df, tmp_df)
+            }
         }
     }
     
+    
     plot_df <- plot_df %>% filter(time <= time_lim) %>% 
         pivot_longer(-c(time, experiment),
                      names_to = "compartment", 
@@ -134,6 +193,17 @@ plot_models <- function(sims = baseline_sim,
     plot_df$Date <- start_date + plot_df$time
     plot_df$time <- NULL
     
+    
+    # Define a standard set of colours to represent compartments
+    comps <- c("s.num", "e.num", "i.num", "q.num", "h.num", "r.num", "f.num")
+    compcols <- c(s.num = "#4477AA", e.num = "#66CCEE", i.num = "#CCBB44", 
+                  q.num = "#AA3377", h.num = "#EE6677", r.num = "#228833", 
+                  f.num = "#BBBBBB")
+    complabels <- c(s.num = "S: Susceptible", e.num = "E: Asymptomatic", 
+                    i.num = "I: Infected", q.num = "Q: Self-isolated", 
+                    h.num = "H: Hospitalized", r.num = "R: Recovered",
+                    f.num = "F: Case Fatalities")
+    
     # Add known data
     if(is.data.frame(known)){
         exps <- unique(plot_df$experiment)
@@ -146,7 +216,7 @@ plot_models <- function(sims = baseline_sim,
     # Filter compartments
     comp_plot <- setdiff(comps, comp_remove)
     plot_df <- plot_df %>% filter(compartment %in% c(comp_plot))
-
+    
     reo_exp <- function(x) { factor(x, levels = sim_id)}
     
     
@@ -197,3 +267,5 @@ plot_models <- function(sims = baseline_sim,
         )
     }
 }
+
+
diff --git a/R/summary.R b/R/FIN_summary.R
similarity index 93%
rename from R/summary.R
rename to R/FIN_summary.R
index dfc0b2bcf5350c6ce7d4490443ef4d7bca2f305f..b8902b9a4973fa6b73764d9c4f5bf8ee4fe9f758 100644
--- a/R/summary.R
+++ b/R/FIN_summary.R
@@ -2,7 +2,8 @@
 #'
 #' Function to extract mean and sdandard deviation of all compartments for SEIQHRF model
 #'
-#' @param object ICM parameters.
+#' @param object seiqhrf object, returned by \code{\link{seiqhrf}}.
+#' @param ... Additional parameters.
 #'
 #' @return A list of compartments, each compartment is a list of three: 
 #' \itemize{
@@ -14,7 +15,7 @@
 #' where t is the total number of time points in the simulation (i.e. object$nstep).
 #'
 #' @export
-summary_seiqhrf <- function(object){
+summary.seiqhrf <- function(object, ...){
     
     nsteps <- object$control$nsteps
     prev_names <- c("s.num", "e.num", "i.num", "q.num", "h.num", "r.num", "f.num")
diff --git a/R/internal_update.R b/R/INTN_update.R
similarity index 100%
rename from R/internal_update.R
rename to R/INTN_update.R
diff --git a/R/PREP_init.R b/R/PREP_init.R
index 7f8147be1816fdede161b2a5f6b72919c2867976..4c311b04442a2eb1efcbeeeb4f2893569e7a349b 100644
--- a/R/PREP_init.R
+++ b/R/PREP_init.R
@@ -37,7 +37,7 @@
 #' @export
 #'
 
-init_seiqhrf <- function(s.num, e.num=0, i.num = 0, q.num=0, h.num=0, r.num = 0, f.num = 0, ...) {
+init_seiqhrf <- function(s.num = 9997, e.num=0, i.num = 3, q.num=0, h.num=0, r.num = 0, f.num = 0, ...) {
 
     p <- list()
     formal.args <- formals(sys.function())
diff --git a/R/original_plot_models.R b/R/original_plot_models.R
new file mode 100644
index 0000000000000000000000000000000000000000..68422188c4219f1fa6ca7d59e697fbdd33a0d1e0
--- /dev/null
+++ b/R/original_plot_models.R
@@ -0,0 +1,123 @@
+#' Plot models
+#'
+#' Function to plot individuals models or mutliple models for comparison.
+#'
+#' @param sims Single or list of sims to plot
+#' @param sim_id String or list of strings to name each facet
+#' @param comp_remove Compartments to remove. Suggest c(s.num, r.num)
+#' @param time_lim Number of steps (days) to plot. 
+#' @param trans Y-axis transformation (e.g. log2, log10). Default = none. 
+#' @param known Dataframe with known compartment numbers to plot alongside
+#'        projections
+#' @param start_date Date for day 0. Default: ymd("2020-03-21"),
+#' @param x_axis Title for x-axis. Default: 'Days since beginning of epidemic'
+#' @param plot_title Title for whole plot. Default: 'SEIQHRF plot'
+#'
+#' @return ggplot2 object
+#' 
+#' @import dplyr
+#' @importFrom tidyr pivot_longer
+
+plot_models <- function(sims = baseline_sim,
+                        sim_id = "baseline",
+                        comp_remove = "none",
+                        time_lim = 100,
+                        trans = 'na',
+                        known = NULL,
+                        start_date = ymd("2020-03-21"),
+                        x_axis = 'Days since beginning of epidemic',
+                        plot_title = 'SEIQHRF'){
+    
+    # Define a standard set of colours to represent compartments
+    comps <- c("s.num", "e.num", "i.num", "q.num", "h.num", "r.num", "f.num")
+    compcols <- c(s.num = "#4477AA", e.num = "#66CCEE", i.num = "#CCBB44", 
+                  q.num = "#AA3377", h.num = "#EE6677", r.num = "#228833", 
+                  f.num = "#BBBBBB")
+    complabels <- c(s.num = "S: Susceptible", e.num = "E: Asymptomatic", 
+                    i.num = "I: Infected", q.num = "Q: Self-isolated", 
+                    h.num = "H: Hospitalized", r.num = "R: Recovered",
+                    f.num = "F: Case Fatalities")
+    
+    # Merge models to plot together
+    for (i in (1: length(sim_id))) {
+        if (i == 1) {
+            plot_df <- sims[i*2]$df
+            plot_df <- plot_df %>% mutate(experiment = sim_id[i])
+        } else{
+            tmp_df <- sims[i*2]$df
+            tmp_df <- tmp_df %>% mutate(experiment = sim_id[i])
+            plot_df <- plot_df %>% bind_rows(plot_df, tmp_df)
+        }
+    }
+    
+    plot_df <- plot_df %>% filter(time <= time_lim) %>% 
+        pivot_longer(-c(time, experiment),
+                     names_to = "compartment", 
+                     values_to = "count")
+    plot_df$sim <- 'sim'
+    plot_df$Date <- start_date + plot_df$time
+    plot_df$time <- NULL
+    
+    # Add known data
+    if(is.data.frame(known)){
+        exps <- unique(plot_df$experiment)
+        for(i in exps){
+            known$experiment <- i
+            plot_df <- rbind(plot_df, known)
+        } 
+    }
+    
+    # Filter compartments
+    comp_plot <- setdiff(comps, comp_remove)
+    plot_df <- plot_df %>% filter(compartment %in% c(comp_plot))
+
+    reo_exp <- function(x) { factor(x, levels = sim_id)}
+    
+    
+    # Plot single model
+    if(trans == "na"){
+        if(length(sim_id) == 1){
+            plot_df %>% ggplot(aes(x = Date, y = count+1, colour = compartment, 
+                                   linetype = sim)) + 
+                geom_line(size = 1.5, alpha = 0.8) + 
+                scale_x_date(date_breaks = "1 week", date_labels = "%m-%d") + 
+                scale_colour_manual(values = compcols, labels = complabels) + 
+                labs(title = plot_title, x = x_axis, y = "Prevalence") +
+                theme_bw() + 
+                theme(axis.text.x = element_text(angle = 90))
+        } else (
+            plot_df %>% ggplot(aes(x = Date, y = count, colour = compartment,
+                                   linetype = sim)) + 
+                facet_grid(reo_exp(experiment) ~ ., scale = 'free') + 
+                scale_x_date(date_breaks = "1 week", date_labels = "%m-%d") + 
+                geom_line(size = 1.5, alpha = 0.8) + 
+                scale_colour_manual(values = compcols, labels = complabels) + 
+                labs(title = plot_title, x = x_axis, y = "Prevalence") +
+                theme_bw() + 
+                theme(axis.text.x = element_text(angle = 90))
+        )
+    }else{
+        if(length(sim_id) == 1){
+            plot_df %>% ggplot(aes(x = Date, y = count+1, colour = compartment, 
+                                   linetype = sim)) + 
+                scale_y_continuous(trans = trans) +
+                scale_x_date(date_breaks = "1 week", date_labels = "%m-%d") + 
+                geom_line(size = 1.5, alpha = 0.8) + 
+                scale_colour_manual(values = compcols, labels = complabels) + 
+                labs(title = plot_title, x = x_axis, y = "Prevalence") +
+                theme_bw() + 
+                theme(axis.text.x = element_text(angle = 90))
+        } else (
+            plot_df %>% ggplot(aes(x = Date, y = count+1, colour = compartment, 
+                                   linetype = sim)) + 
+                facet_grid(reo_exp(experiment) ~ ., scale = 'free') + 
+                scale_y_continuous(trans = trans) +
+                scale_x_date(date_breaks = "1 week", date_labels = "%m-%d") + 
+                geom_line(size = 1.5, alpha = 0.8) + 
+                scale_colour_manual(values = compcols, labels = complabels) + 
+                labs(title = plot_title, x = x_axis, y = "Prevalence") +
+                theme_bw() + 
+                theme(axis.text.x = element_text(angle = 90))
+        )
+    }
+}
diff --git a/R/RUN_simulation.R b/R/original_simulation.R
similarity index 99%
rename from R/RUN_simulation.R
rename to R/original_simulation.R
index 07274ae8d64c12f58a44fdc7864fcdc684c4b3d5..1e6783a1a111536c1e4dcac6b5b164cb05f8b531 100644
--- a/R/RUN_simulation.R
+++ b/R/original_simulation.R
@@ -209,6 +209,7 @@
 #'     to model emigration as well as deaths.
 #' @param out Summary function for the simulation runs. median is
 #'     also available, or percentiles, see the EpiModel documentation.
+#' @param seed Random see for checking consistency.    
 #'
 #' @return list with simulation and simulation dataframe
 #' @export
@@ -276,7 +277,7 @@ simulate_seiqhrf <- function(type = "SEIQHRF",
     dq.rate = (7/365)/1000,
     dh.rate = (20/365)/1000,
     dr.rate = (7/365)/1000,
-    out="mean"
+    out="mean", seed = NULL
 ) {
 
     control <- control_seiqhrf(type = type,
@@ -343,7 +344,7 @@ simulate_seiqhrf <- function(type = "SEIQHRF",
                         dh.rate = dh.rate,
                         dr.rate = dr.rate)
 
-    sim <- seiqhrf(param = param, init = init, control = control)
+    sim <- seiqhrf(param = param, init = init, control = control, seed = seed)
     sim_df <- as.data.frame(sim, out=out)
     class(sim) <- "icm"
 
diff --git a/man/get_times.Rd b/man/get_times.Rd
deleted file mode 100644
index e1b667697f01d809cbae66b030325fa204eca8ae..0000000000000000000000000000000000000000
--- a/man/get_times.Rd
+++ /dev/null
@@ -1,17 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/results-parse.R
-\name{get_times}
-\alias{get_times}
-\title{Get data by time}
-\usage{
-get_times(simulate_results)
-}
-\arguments{
-\item{simulate_results}{results from `simulte()` function.}
-}
-\value{
-dataframe with sim.n, period type, and duration
-}
-\description{
-Function to extract timings and assemble results in a dataframe
-}
diff --git a/man/init_seiqhrf.Rd b/man/init_seiqhrf.Rd
index 435306ad3875570b37a08525dd085470ad237694..1676bedff549e9fb27e12f4d71d160c98f4e7e15 100644
--- a/man/init_seiqhrf.Rd
+++ b/man/init_seiqhrf.Rd
@@ -5,9 +5,9 @@
 \title{Initial Conditions for Stochastic Individual Contact Models}
 \usage{
 init_seiqhrf(
-  s.num,
+  s.num = 9997,
   e.num = 0,
-  i.num = 0,
+  i.num = 3,
   q.num = 0,
   h.num = 0,
   r.num = 0,
diff --git a/man/plot.list.Rd b/man/plot.list.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..820a51f108f71c5b17c2a2cb94410541ba8e81dc
--- /dev/null
+++ b/man/plot.list.Rd
@@ -0,0 +1,44 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/FIN_plot.R
+\name{plot.list}
+\alias{plot.list}
+\title{Plot simulation result}
+\usage{
+\method{plot}{list}(
+  x,
+  comp_remove = "none",
+  time_lim = 100,
+  trans = "na",
+  known = NULL,
+  start_date = ymd("2020-03-21"),
+  x_axis = "Days since beginning of epidemic",
+  plot_title = "SEIQHRF",
+  ...
+)
+}
+\arguments{
+\item{x}{A list of seiqhrf objects returned from function \code{\link{seiqhrf}}.}
+
+\item{comp_remove}{Compartments to remove. Suggest c(s.num, r.num)}
+
+\item{time_lim}{Number of steps (days) to plot.}
+
+\item{trans}{Y-axis transformation (e.g. log2, log10). Default = none.}
+
+\item{known}{Dataframe with known compartment numbers to plot alongside
+projections}
+
+\item{start_date}{Date for day 0. Default: ymd("2020-03-21"),}
+
+\item{x_axis}{Title for x-axis. Default: 'Days since beginning of epidemic'}
+
+\item{plot_title}{Title for whole plot. Default: 'SEIQHRF plot'}
+
+\item{...}{Additional parameters}
+}
+\value{
+ggplot2 object
+}
+\description{
+Function to extract timings, prevalence etc. from the simulation
+}
diff --git a/man/plot.seiqhrf.Rd b/man/plot.seiqhrf.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..86009a2d7ea26ca6c20969e122766709a9b08f7d
--- /dev/null
+++ b/man/plot.seiqhrf.Rd
@@ -0,0 +1,48 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/FIN_plot.R
+\name{plot.seiqhrf}
+\alias{plot.seiqhrf}
+\title{Plot simulation result}
+\usage{
+\method{plot}{seiqhrf}(
+  x,
+  method = NULL,
+  comp_remove = "none",
+  time_lim = 100,
+  trans = "na",
+  known = NULL,
+  start_date = ymd("2020-03-21"),
+  x_axis = "Days since beginning of epidemic",
+  plot_title = "SEIQHRF",
+  ...
+)
+}
+\arguments{
+\item{x}{An seiqhrf object returned from function \code{\link{seiqhrf}}.}
+
+\item{method}{If "times", plot Duration frequency distributions.
+If NULL, plot individuals models or mutliple models for comparison.}
+
+\item{comp_remove}{Compartments to remove. Suggest c(s.num, r.num)}
+
+\item{time_lim}{Number of steps (days) to plot.}
+
+\item{trans}{Y-axis transformation (e.g. log2, log10). Default = none.}
+
+\item{known}{Dataframe with known compartment numbers to plot alongside
+projections}
+
+\item{start_date}{Date for day 0. Default: ymd("2020-03-21"),}
+
+\item{x_axis}{Title for x-axis. Default: 'Days since beginning of epidemic'}
+
+\item{plot_title}{Title for whole plot. Default: 'SEIQHRF plot'}
+
+\item{...}{Additional parameters}
+}
+\value{
+ggplot2 object
+}
+\description{
+Function to extract timings, prevalence etc. from the simulation
+}
diff --git a/man/plot_models.Rd b/man/plot_models.Rd
index a6a19de28f2683d0de54f3ddb9b943482d8bb159..c158b0e2146b997628d4cd950657fa1d2e628ffb 100644
--- a/man/plot_models.Rd
+++ b/man/plot_models.Rd
@@ -1,5 +1,5 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/results-parse.R
+% Please edit documentation in R/original_plot_models.R
 \name{plot_models}
 \alias{plot_models}
 \title{Plot models}
@@ -13,7 +13,7 @@ plot_models(
   known = NULL,
   start_date = ymd("2020-03-21"),
   x_axis = "Days since beginning of epidemic",
-  plot_title = "ICM plot"
+  plot_title = "SEIQHRF"
 )
 }
 \arguments{
@@ -34,7 +34,7 @@ projections}
 
 \item{x_axis}{Title for x-axis. Default: 'Days since beginning of epidemic'}
 
-\item{plot_title}{Title for whole plot. Default: 'ICM plot'}
+\item{plot_title}{Title for whole plot. Default: 'SEIQHRF plot'}
 }
 \value{
 ggplot2 object
diff --git a/man/plot_times.Rd b/man/plot_times.Rd
index cc2a1de0147742a02566fbcee5deda92dca606f0..ca596ed35650bb790f01c11fdfeb6187356c02c6 100644
--- a/man/plot_times.Rd
+++ b/man/plot_times.Rd
@@ -1,13 +1,13 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/results-parse.R
+% Please edit documentation in R/FIN_plot.R
 \name{plot_times}
 \alias{plot_times}
-\title{Plot times}
+\title{Plot simulation by time}
 \usage{
-plot_times(times)
+plot_times(sim)
 }
 \arguments{
-\item{times}{results from `get_times()` function.}
+\item{sim}{An seiqhrf object returned from function \code{\link{seiqhrf}}.}
 }
 \value{
 ggplot2 object
diff --git a/man/simulate_seiqhrf.Rd b/man/simulate_seiqhrf.Rd
index 4e4ae0e602624011bf5608603397daca760c04bc..782e58034c6313fac08397f462426045c8c72af7 100644
--- a/man/simulate_seiqhrf.Rd
+++ b/man/simulate_seiqhrf.Rd
@@ -1,5 +1,5 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/RUN_simulation.R
+% Please edit documentation in R/original_simulation.R
 \name{simulate_seiqhrf}
 \alias{simulate_seiqhrf}
 \title{SEIQHRF Simulation Wrapper}
@@ -67,7 +67,8 @@ simulate_seiqhrf(
   dq.rate = (7/365)/1000,
   dh.rate = (20/365)/1000,
   dr.rate = (7/365)/1000,
-  out = "mean"
+  out = "mean",
+  seed = NULL
 )
 }
 \arguments{
@@ -340,6 +341,8 @@ to model emigration as well as deaths.}
 
 \item{out}{Summary function for the simulation runs. median is
 also available, or percentiles, see the EpiModel documentation.}
+
+\item{seed}{Random see for checking consistency.}
 }
 \value{
 list with simulation and simulation dataframe
diff --git a/man/summary_seiqhrf.Rd b/man/summary.seiqhrf.Rd
similarity index 75%
rename from man/summary_seiqhrf.Rd
rename to man/summary.seiqhrf.Rd
index 27b85cb4b75390422829232e6ae78b51a49571fd..3e68728c0c27cb439931b3edf15caf7813163067 100644
--- a/man/summary_seiqhrf.Rd
+++ b/man/summary.seiqhrf.Rd
@@ -1,13 +1,15 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/summary.R
-\name{summary_seiqhrf}
-\alias{summary_seiqhrf}
+% Please edit documentation in R/FIN_summary.R
+\name{summary.seiqhrf}
+\alias{summary.seiqhrf}
 \title{summary function}
 \usage{
-summary_seiqhrf(object)
+\method{summary}{seiqhrf}(object, ...)
 }
 \arguments{
-\item{object}{ICM parameters.}
+\item{object}{seiqhrf object, returned by \code{\link{seiqhrf}}.}
+
+\item{...}{Additional parameters.}
 }
 \value{
 A list of compartments, each compartment is a list of three: 
diff --git a/tests/testthat/.DS_Store b/tests/testthat/.DS_Store
index 41c1f5f62a267c92f511c48760e81d53b16b7820..9e1fc75f6bee687fada1a430d6e0d1018cdd04a1 100644
Binary files a/tests/testthat/.DS_Store and b/tests/testthat/.DS_Store differ
diff --git a/tests/testthat/test-seiqhrf.R b/tests/testthat/test-all-modules.R
similarity index 91%
rename from tests/testthat/test-seiqhrf.R
rename to tests/testthat/test-all-modules.R
index 615b79d7530c1cf4f1a327a5bc7d90eed8f4d4ad..acda5c05a2195580619a0807fc37aa387e75f577 100644
--- a/tests/testthat/test-seiqhrf.R
+++ b/tests/testthat/test-all-modules.R
@@ -1,4 +1,4 @@
-test_that("Identical output as Churches' original function: seiqhrf", {
+test_that("All modules(MOD) produce the same result as Churches' modules (.seiqhrf.icm)", {
     
     full_params <- set_param(s.num = 1000, nsteps = 10, arec.rate = 0)
     param <- full_params$param
diff --git a/tests/testthat/test-plot-model.R b/tests/testthat/test-plot-model.R
new file mode 100644
index 0000000000000000000000000000000000000000..a42de9683df813ceef53f5a0ee85b06cb1ff0076
--- /dev/null
+++ b/tests/testthat/test-plot-model.R
@@ -0,0 +1,49 @@
+test_that("plot_model is producing the same plot", {
+
+s.num <- 2000  # number susceptible
+i.num <- 15  # number infected 
+q.num <- 5  # number in self-isolation
+h.num <- 1  #     
+vals <- c(10, 7)
+timing <- c(7, 14)
+nsteps <- 90
+
+
+## baseline
+baseline_sim <- simulate_seiqhrf(nsteps = nsteps, s.num = s.num, i.num = i.num,
+                                 q.num = q.num, h.num = h.num)
+
+## Experiment 1
+act_rate <- vary_param(nstep = nsteps, vals = vals, timing = timing)
+
+sim_exp <- simulate_seiqhrf(nsteps = nsteps, s.num = s.num, i.num = i.num,
+                            q.num = q.num, h.num = h.num, act.rate.e = act_rate, 
+                            act.rate.i = act_rate * 0.5)
+
+# Experiment #2
+vals <- c(10, 7, 7, 10)
+timing <- c(7, 14, 21, 28)
+act_rate_relax <- vary_param(nsteps, vals = vals, timing = timing)
+
+sim_exp_relax <- simulate_seiqhrf(nsteps = nsteps, s.num = s.num, i.num = i.num,
+                                  q.num = q.num, h.num = h.num, 
+                                  act.rate.e = act_rate_relax, 
+                                  act.rate.i = 0.5 * act_rate_relax)
+
+
+
+ori_plot <- plot_models(c(baseline_sim, sim_exp, sim_exp_relax),
+                        sim_id = c('Baseline', 'Closures', 'Closures (2 mo)'),
+                        start_date = lubridate::ymd("2020-01-01"),
+                        comp_remove = c('s.num', 'r.num'),
+                        plot_title = 'Closures Experiment')
+
+my_plot <-  plot(list('Baseline' = baseline_sim$sim, 'Closures' = sim_exp$sim, 'Closures (2 mo)' = sim_exp_relax$sim),
+                 start_date = lubridate::ymd("2020-01-01"),
+                 comp_remove = c('s.num', 'r.num'),
+                 plot_title = 'Closures Experiment')
+
+
+expect_equal(all.equal(ori_plot, my_plot), TRUE)
+
+})
diff --git a/tests/testthat/test-seiqhrf-simulate.R b/tests/testthat/test-seiqhrf-simulate.R
new file mode 100644
index 0000000000000000000000000000000000000000..cb9158e662426a2c074e823b4acc962cb4106115
--- /dev/null
+++ b/tests/testthat/test-seiqhrf-simulate.R
@@ -0,0 +1,32 @@
+test_that("seiqhrf and simulate_seiqhrf produce identical output", {
+    
+    s.num = 1000
+    q.num = 10
+    nsteps = 10
+    nsims = 3
+    arec.rate = 0  ### has to be fixed 0 for comparison
+    
+    No_seeds <- 10
+    seed_list <- sample(1:1000, No_seeds)
+    comp <- rep(NA, No_seeds)
+    i <- 1
+    for(seed in seed_list){
+        Churhes_res <- simulate_seiqhrf(nsteps = nsteps, nsims = nsims,
+                                        arec.rate = arec.rate, s.num = s.num, q.num = q.num,
+                                        infection.FUN = infection.seiqhrf.icm, 
+                                        recovery.FUN = progress.seiqhrf.icm, 
+                                        departures.FUN = departures.seiqhrf.icm, 
+                                        arrivals.FUN = arrivals.seiqhrf.icm, seed = seed)$sim
+        class(Churhes_res) <- "seiqhrf"
+        
+        param <- param_seiqhrf(arec.rate = arec.rate)
+        init <- init_seiqhrf(s.num = s.num, q.num = q.num)
+        control <- control_seiqhrf(nsteps = nsteps, nsims = nsims)
+        sirplus_res <- seiqhrf(init, control, param, seed)
+        
+        comp[i] <- identical(Churhes_res, sirplus_res)
+        i <- i + 1
+    }
+    
+    expect_equal(sum(comp), No_seeds)
+})
diff --git a/vignettes/sirplus_intro.Rmd b/vignettes/sirplus_intro.Rmd
index 6dec74ffb89776c55ab18666cba0ce778cfa5183..07d04fe9503f8d9820a286364e9ade737749294b 100644
--- a/vignettes/sirplus_intro.Rmd
+++ b/vignettes/sirplus_intro.Rmd
@@ -75,8 +75,7 @@ functions to examine the distributions of timings for various transitions based
 on these parameters. In the case of a disease with observed data available, these plots can be used to sanity check parameter settings.
 
 ```{r baseline sims, fig.height=8}
-times <- get_times(baseline_sim)
-plot_times(times)
+plot_times(baseline_sim$sim)
 ```
 
 
diff --git a/vignettes/sirplus_intro_v2.Rmd b/vignettes/sirplus_intro_v2.Rmd
index 4a9bdd771ceed584256bb5bdd3f8f6db47a6c09a..2594e00656642a942a951f13c2cc6638e3930e8a 100644
--- a/vignettes/sirplus_intro_v2.Rmd
+++ b/vignettes/sirplus_intro_v2.Rmd
@@ -16,15 +16,54 @@ knitr::opts_chunk$set(echo = TRUE, cache=FALSE, eval=TRUE, collapse = TRUE,
 
 ![](sirplus-logo-small.png)
 
-```{r, set parameters}
+The sirplus package makes it easy to generate stochastic individual compartment 
+models (ICMs) to simulate contagious disease spread using compartments not 
+available in standard SIR packages. This method and most of the code was 
+originally written by Tim Churches (see his [blog post](https://timchurches.github.io/blog/posts/2020-03-18-modelling-the-effects-of-public-health-interventions-on-covid-19-transmission-part-2/)). The sirplus package was developed by 
+the Bioinformatics & Cellular Genomics team at St. Vincent's Institute of 
+Medical Research in order to help St. Vincents' Hospital model the COVID-19
+pandemic.
+
+The compartments available in this package include: 
+
+| Compartment | Functional definition                                                             |
+|-------------|-----------------------------------------------------------------------------------|
+| S           | Susceptible individuals                                                           |
+| E           | Exposed **and** infected, not yet symptomatic but potentially infectious          |
+| I           | Infected, symptomatic **and** infectious                                          |
+| Q           | Infectious, but (self-)isolated                                                   |
+| H           | Requiring hospitalisation (would normally be hospitalised if capacity available)  |
+| R           | Recovered, immune from further infection                                          |
+| F           | Case fatality (death due to infection, not other causes)                          |
+
+Again, for more information about these compartments and about the parameters 
+that are associated with them see Tim Churches's [blog post](https://timchurches.github.io/blog/posts/2020-03-18-modelling-the-effects-of-public-health-interventions-on-covid-19-transmission-part-2/).
+
+```{r, load package}
 library(sirplus)
-s.num <- 1000  # number susceptible
+devtools::load_all()
+```
+
+
+## Simulate and inspect a baseline sirplus model
+
+### Set parameters
+
+Here we will simulate the epidemiological data for a made-up population with 1000 
+susceptible individuals (S), 50 that are infected but not in the hospital or in 
+self-quarantine (I; maybe people that are infected/symptomatic but not tested/
+aware), 10 confirmed cases that have self-isolated (Q), and 1 confirmed case that
+has been hospitalized (H). We call this the baseline model because it uses 
+default parameters for disease spread (i.e. no additional interventions).
+
+```{r, set parameters}
+s.num <- 2000  # number susceptible
 i.num <- 15  # number infected 
 q.num <- 5  # number in self-isolation
 h.num <- 1  # number in the hospital
 nsteps <- 90 # number of steps (e.g. days) to simulate
 control <- control_seiqhrf(nsteps = nsteps)
-param <- param_seiqhrf(arec.rate = 0)
+param <- param_seiqhrf()
 init <- init_seiqhrf(s.num = s.num, i.num = i.num, q.num = q.num, h.num = h.num)
 
 print(init)
@@ -32,7 +71,42 @@ print(control)
 print(param)
 ```
 
+### Simulate baseline
+
+This will produce an seiqhrf object. 
+
 ```{r, run simulation}
 sim <- seiqhrf(init, control, param)
 sim
-```
\ No newline at end of file
+```
+
+### Extract summary of the simulation
+
+```{r, extract summary}
+res <- summary(sim)
+names(res)
+```
+
+### Inspect baseline transition distributions
+
+The sirplus model controls transitions between compartments, i.e. a change in state for an individual (e.g. going from self-isolation to hospital), using a variety of 
+transition parameters. You can use the `plot()` functions, and set parameter `method` as `times` 
+to examine the distributions of timings for various transitions based
+on these parameters. In the case of a disease with observed data available, these plots can be used to sanity check parameter settings.
+
+
+```{r baseline sims, fig.height=8}
+plot(sim, "times")
+```
+
+### Plot baseline sirplus results
+
+To visualise your sirplus model, you can plot the change in prevalence (i.e. people) 
+over time in each compartment by changing parameter `method` to `models`. 
+
+```{r viz prevalence}
+plot(sim, 
+     start_date = lubridate::ymd("2020-01-01"),
+     comp_remove = c('s.num', 'r.num'),
+     plot_title = 'Baseline Model')
+```