diff --git a/DESCRIPTION b/DESCRIPTION
index f50f6b5f8..9c486eb45 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: macpan2
 Title: Fast and Flexible Compartmental Modelling
-Version: 1.14.0
+Version: 2.0.0
 Authors@R: c(
   person("Steve Walker", email="swalk@mcmaster.ca", role=c("cre", "aut")),
   person("Weiguang Guan", role="aut"),
diff --git a/R/formula_list_generators.R b/R/formula_list_generators.R
index 1bdec0a96..9e6340873 100644
--- a/R/formula_list_generators.R
+++ b/R/formula_list_generators.R
@@ -352,8 +352,8 @@ MockChangeModel = function() {
 ##' @examples
 ##' sir = mp_tmb_library("starter_models", "sir", package = "macpan2")
 ##' sir
-##' sir |> mp_euler()             |> mp_expand()
-##' sir |> mp_rk4()               |> mp_expand()
+##' sir |> mp_euler() |> mp_expand()
+##' sir |> mp_rk4() |> mp_expand()
 ##' sir |> mp_euler_multinomial() |> mp_expand()
 ##' 
 ##' @name state_updates
diff --git a/R/mp_tmb_model_spec.R b/R/mp_tmb_model_spec.R
index 31d69d1f3..ef74315e1 100644
--- a/R/mp_tmb_model_spec.R
+++ b/R/mp_tmb_model_spec.R
@@ -299,9 +299,10 @@ must_save_time_args = function(formulas) {
 #' being evaluated. For example, expressions that generate stochasticity should
 #' be listed in \code{sim_exprs} because TMB objective functions must be
 #' continuous.
-#' @param state_update (experimental) Optional character vector for how to 
-#' update the state variables when it is relevant. Options include `"euler"`, 
-#' `"rk4"`, and `"euler_multinomial"`.
+#' @param state_update Optional character vector for how to update the state 
+#' variables when it is relevant. Can be the name of any 
+#' \code{\link{state_updates}} function with the `mp_` prefix removed (e.g., 
+#' `"euler"`, `"rk4"`, and `"euler_multinomial"`).
 #' 
 #' @examples
 #' spec = mp_tmb_model_spec(
diff --git a/R/tmb_model_editors.R b/R/tmb_model_editors.R
index 0cc47caf4..e9eb61307 100644
--- a/R/tmb_model_editors.R
+++ b/R/tmb_model_editors.R
@@ -199,13 +199,34 @@ mp_tmb_delete = function(model
 #' @param incidence_name Name of the incidence variable to be transformed.
 #' @param report_prob Value to use for the reporting probability; the
 #' proportion of cases that get reported.
-#' @param mean_delay Mean of the Gamma distribution of reporting delay times.
+#' @param mean_delay Mean of the Gamma distribution of reporting delay times
+#' measured in time steps. For example, if the mean delay time in days is 14
+#' and each time step in the model corresponds to one week, then `mean_delay`
+#' should be set to `2`.
 #' @param cv_delay Coefficient of variation of the Gamma distribution of
-#' reporting delay times.
+#' reporting delay times, with time measured in time steps (see more on this
+#' in the description of the `mean_delay` argument immediately above).
 #' @param reports_name Name of the new reports variable.
 #' @param report_prob_name Name of the variable containing `report_prob`.
 #' @param mean_delay_name Name of the variable containing `mean_delay`.
 #' @param cv_delay_name Name of the variable containing `cv_delay`.
+#' @examples
+#' sir = mp_tmb_library("starter_models", "sir", package = "macpan2")
+#' sir_delays = mp_tmb_insert_reports(sir
+#'   , incidence_name = "infection"
+#'   , report_prob = 0.1
+#'   , mean_delay = 14
+#'   , cv_delay = 0.25
+#'   , reports_name = "reports"
+#' )
+#' sir_delays |> mp_print_during()
+#' (sir_delays 
+#'   |> mp_simulator(
+#'       time_steps = 50L
+#'     , outputs = c("infection", "reports"), 50L)
+#'   |> mp_trajectory()
+#' )
+#' 
 #' @export
 mp_tmb_insert_reports = function(model
   , incidence_name
@@ -231,7 +252,7 @@ mp_tmb_insert_reports = function(model
   kernel_length = qgamma(0.95, shape, scale = scale) |> ceiling() |> as.integer()
 
   expressions = c(
-      sprintf("%s ~ pgamma(1:%s, 1/(%s), %s * (%s^2))", map$dist, kernel_length + 1L, cv_delay_name, mean_delay_name, cv_delay_name)
+      sprintf("%s ~ pgamma(1:%s, 1/(%s^2), %s * (%s^2))", map$dist, kernel_length + 1L, cv_delay_name, mean_delay_name, cv_delay_name)
     , sprintf("%s ~ %s[1:%s] - %s[0:%s]", map$delta, map$dist, kernel_length, map$dist, kernel_length - 1L)
     , sprintf("%s ~ %s * %s / sum(%s)", map$kernel, report_prob_name, map$delta, map$delta)
     , sprintf("%s ~ convolution(%s, %s)", reports_name, incidence_name, map$kernel)
diff --git a/inst/starter_models/shiver/README.Rmd b/inst/starter_models/shiver/README.Rmd
index 547f5a2ea..fcfa67fb5 100644
--- a/inst/starter_models/shiver/README.Rmd
+++ b/inst/starter_models/shiver/README.Rmd
@@ -338,9 +338,9 @@ Now we update the specification to include these defaults.
 spec = mp_tmb_update(spec
   , default = list(
       N = N
-    , V = V0
-    , I = I0
-    , H = H0
+    , V = 0 # V0
+    , I = 1 # I0
+    , H = 0 # H0
     , a = a
     , b = a
     , alpha = alpha
@@ -369,6 +369,7 @@ shiver_calibrator = mp_tmb_calibrator(
   # we also want to estimate initial E
   , par = c("beta_v","beta_s","E", "sigma", "gamma_h") 
   , outputs = states
+  , time = mp_sim_bounds(-15, 90, "steps")
 )
 # print to check
 shiver_calibrator
@@ -658,8 +659,8 @@ prior_distributions = list(
   , logit_report_prob = mp_normal(qlogis(0.1), sd_par)
   , logit_p = mp_normal(qlogis(1/4), 4)
   , log_E_I_ratio = mp_normal(0, sd_par)
-  , log_I = mp_normal(log(I0), sd_state)
-  , log_H = mp_normal(log(H0), sd_state)
+  , log_I = mp_normal(log(1), sd_state)
+  , log_H = mp_normal(log(1), sd_state)
   , log_R = mp_normal(log(1), sd_state)
 )
 
@@ -669,6 +670,9 @@ dd = rbind(reported_hospitalizations, reported_cases)
 # calibrate
 shiver_calibrator = mp_tmb_calibrator(
     spec = (multi_traj_spec 
+      |> mp_tmb_update(
+        #default = list(log_I = 0, log_H = 0, log_R = 0)
+      )
       |> mp_hazard()
     )
     # row bind both observed data
@@ -684,6 +688,7 @@ shiver_calibrator = mp_tmb_calibrator(
     # a flexible model of time variation.
   , tv = mp_rbf("rbf_beta", 5, sparse_tol = 1e-8)
   , outputs = c(states, "reported_incidence", "beta")
+  , time = mp_sim_bounds(-30, 90, "steps")
 )
 ```
 
@@ -713,7 +718,7 @@ Our prior for `sigma` is similar to the posterior, but `gamma_h` seems to have b
       , fill = "red"
       , alpha = 0.3
     )
-   + geom_point(data = dd, aes(time, value))
+   + geom_point(data = mutate(dd, time = time + 30L), aes(time, value))
    + facet_wrap(vars(matrix), scales = 'free')
    + theme_bw()
 )
diff --git a/inst/starter_models/shiver/README.md b/inst/starter_models/shiver/README.md
index fb1157408..a237f5069 100644
--- a/inst/starter_models/shiver/README.md
+++ b/inst/starter_models/shiver/README.md
@@ -41,7 +41,7 @@ Vaccines are typically subject to resource constraints and distribution
 strategies might prioritize vaccinations for specific subpopulations,
 such as immunocompromised people, to reduce bad outcomes. We model this
 with a flow of susceptibles entering the vaccination class. This flow
-could be a fixed rate, i.e. a constant proportion of the population
+could be a fixed rate, i.e. a constant proportion of the population
 receives a vaccine each time step, but instead we wish to capture a more
 realistic vaccination rate by allowing it to vary (see [Variable
 Vaccination Rate](#variable-vaccination-rate)).
@@ -138,7 +138,7 @@ We can implement vaccine constraints by adding more model complexity.
 Resource limitations create an upper bound on the number of vaccines
 that can be administered to susceptibles per time step. There is also
 the constraint that we can only vaccinate, at most, the current number
-of susceptibles i.e. the vaccination rate can be at most 1. These
+of susceptibles i.e. the vaccination rate can be at most 1. These
 constraints naturally lead us to consider a variable vaccination rate
 $\phi(S(t))$, instead of vaccinating a fixed proportion $\phi > 0$ per
 time step.
diff --git a/inst/starter_models/shiver/figures/diagram-1.png b/inst/starter_models/shiver/figures/diagram-1.png
index 9dfba65f3..083a2eecf 100644
Binary files a/inst/starter_models/shiver/figures/diagram-1.png and b/inst/starter_models/shiver/figures/diagram-1.png differ
diff --git a/inst/starter_models/shiver/figures/mult_traj_fit-1.png b/inst/starter_models/shiver/figures/mult_traj_fit-1.png
index 0798f2e86..e04851e22 100644
Binary files a/inst/starter_models/shiver/figures/mult_traj_fit-1.png and b/inst/starter_models/shiver/figures/mult_traj_fit-1.png differ
diff --git a/inst/starter_models/sir/README.md b/inst/starter_models/sir/README.md
index 3b11aa7dd..7390eb6d1 100644
--- a/inst/starter_models/sir/README.md
+++ b/inst/starter_models/sir/README.md
@@ -150,12 +150,12 @@ print(spec |> mp_expand())
 #> ---------------------
 #> Default values:
 #> ---------------------
-#>  matrix row col value
-#>    beta           0.2
-#>   gamma           0.1
-#>       N         100.0
-#>       I           1.0
-#>       R           0.0
+#>  quantity value
+#>      beta   0.2
+#>     gamma   0.1
+#>         N 100.0
+#>         I   1.0
+#>         R   0.0
 #> 
 #> ---------------------
 #> Before the simulation loop (t = 0):
@@ -165,7 +165,7 @@ print(spec |> mp_expand())
 #> ---------------------
 #> At every iteration of the simulation loop (t = 1 to T):
 #> ---------------------
-#> 1: infection ~ S * (I * beta/N)
+#> 1: infection ~ S * (beta * I/N)
 #> 2: recovery ~ I * (gamma)
 #> 3: S ~ S - infection
 #> 4: I ~ I + infection - recovery
diff --git a/inst/starter_models/sir/figures/diagram-1.png b/inst/starter_models/sir/figures/diagram-1.png
index 1157c89e5..f18eeaf50 100644
Binary files a/inst/starter_models/sir/figures/diagram-1.png and b/inst/starter_models/sir/figures/diagram-1.png differ
diff --git a/man/mp_tmb_insert_reports.Rd b/man/mp_tmb_insert_reports.Rd
index 91966846b..e54d24ff0 100644
--- a/man/mp_tmb_insert_reports.Rd
+++ b/man/mp_tmb_insert_reports.Rd
@@ -24,10 +24,14 @@ mp_tmb_insert_reports(
 \item{report_prob}{Value to use for the reporting probability; the
 proportion of cases that get reported.}
 
-\item{mean_delay}{Mean of the Gamma distribution of reporting delay times.}
+\item{mean_delay}{Mean of the Gamma distribution of reporting delay times
+measured in time steps. For example, if the mean delay time in days is 14
+and each time step in the model corresponds to one week, then \code{mean_delay}
+should be set to \code{2}.}
 
 \item{cv_delay}{Coefficient of variation of the Gamma distribution of
-reporting delay times.}
+reporting delay times, with time measured in time steps (see more on this
+in the description of the \code{mean_delay} argument immediately above).}
 
 \item{reports_name}{Name of the new reports variable.}
 
@@ -45,3 +49,21 @@ a convolution of the simulation history of an incidence variable with
 a kernel that is proportional to a Gamma distribution of reporting
 delay times.
 }
+\examples{
+sir = mp_tmb_library("starter_models", "sir", package = "macpan2")
+sir_delays = mp_tmb_insert_reports(sir
+  , incidence_name = "infection"
+  , report_prob = 0.1
+  , mean_delay = 14
+  , cv_delay = 0.25
+  , reports_name = "reports"
+)
+sir_delays |> mp_print_during()
+(sir_delays 
+  |> mp_simulator(
+      time_steps = 50L
+    , outputs = c("infection", "reports"), 50L)
+  |> mp_trajectory()
+)
+
+}
diff --git a/man/mp_tmb_model_spec.Rd b/man/mp_tmb_model_spec.Rd
index d481e6bd8..e97896143 100644
--- a/man/mp_tmb_model_spec.Rd
+++ b/man/mp_tmb_model_spec.Rd
@@ -65,9 +65,10 @@ being evaluated. For example, expressions that generate stochasticity should
 be listed in \code{sim_exprs} because TMB objective functions must be
 continuous.}
 
-\item{state_update}{(experimental) Optional character vector for how to
-update the state variables when it is relevant. Options include \code{"euler"},
-\code{"rk4"}, and \code{"euler_multinomial"}.}
+\item{state_update}{Optional character vector for how to update the state
+variables when it is relevant. Can be the name of any
+\code{\link{state_updates}} function with the \code{mp_} prefix removed (e.g.,
+\code{"euler"}, \code{"rk4"}, and \code{"euler_multinomial"}).}
 }
 \description{
 Specify a simulation model in the TMB engine.
diff --git a/man/state_updates.Rd b/man/state_updates.Rd
index 47c9eb699..1e9189d05 100644
--- a/man/state_updates.Rd
+++ b/man/state_updates.Rd
@@ -86,8 +86,8 @@ distribution.
 \examples{
 sir = mp_tmb_library("starter_models", "sir", package = "macpan2")
 sir
-sir |> mp_euler()             |> mp_expand()
-sir |> mp_rk4()               |> mp_expand()
+sir |> mp_euler() |> mp_expand()
+sir |> mp_rk4() |> mp_expand()
 sir |> mp_euler_multinomial() |> mp_expand()
 
 }