diff --git a/R/change_components.R b/R/change_components.R
index 720515c3..51fc8891 100644
--- a/R/change_components.R
+++ b/R/change_components.R
@@ -21,20 +21,20 @@ print_flow_vec = function(vec) {
 #' @export
 print.PerCapitaFlow = function(x, ...) {
   c(
-      From = x$from
+     `Flow name` = lhs_char(x$rate)
+    , From = x$from
     , To = x$to
     , `Per-capita rate expression` = rhs_char(x$rate)
-    , `Absolute rate symbol` = lhs_char(x$rate)
   ) |> print_flow_vec()
 }
 
 #' @export
 print.AbsoluteFlow = function(x, ...) {
   c(
-      From = x$from
+      `Flow name` = lhs_char(x$rate)
+    , From = x$from
     , To = x$to
     , `Absolute rate expression` = rhs_char(x$rate)
-    , `Absolute rate symbol` = lhs_char(x$rate)
   ) |> print_flow_vec()
 }
 
@@ -42,18 +42,18 @@ print.AbsoluteFlow = function(x, ...) {
 #' @export
 print.AbsoluteInFlow = function(x, ...) {
   c(
-      To = x$to
+      `Flow name` = lhs_char(x$rate)
+    , To = x$to
     , `Absolute rate expression` = rhs_char(x$rate)
-    , `Absolute rate symbol` = lhs_char(x$rate)
   ) |> print_flow_vec()
 }
 
 #' @export
 print.AbsoluteOutFlow = function(x, ...) {
   c(
-      From = x$from
+      `Flow name` = lhs_char(x$rate)
+    , From = x$from
     , `Absolute rate expression` = rhs_char(x$rate)
-    , `Absolute rate symbol` = lhs_char(x$rate)
   ) |> print_flow_vec()
 }
 
diff --git a/R/change_models.R b/R/change_models.R
index 70448e3d..08d2a93f 100644
--- a/R/change_models.R
+++ b/R/change_models.R
@@ -1,7 +1,9 @@
 
-ChangeModelDefaults = function() {
+ChangeModelDefaults = function(delta_t = NULL) {
   self = ChangeModel()
   
+  self$delta_t = handle_delta_t(delta_t)
+  
   self$change_list = list()
   
   self$flow_frame = function() {
diff --git a/R/formula_list_generators.R b/R/formula_list_generators.R
index ed67f0a2..733dafe7 100644
--- a/R/formula_list_generators.R
+++ b/R/formula_list_generators.R
@@ -111,6 +111,8 @@ only_iterable = function(expr_list, states, is_first = FALSE) {
 ChangeModel = function() {
   self = Base()
   
+  self$delta_t = NULL
+  
   # lists of formula expressions to be added to a `before` list
   self$before_loop = function() list()
   
@@ -196,7 +198,8 @@ ChangeComponent = function() {
   ## column - change: unsigned absolute flow rate name.
   ## column - rate: per-capita flow rates (variables or expressions that 
   ## sometimes involve state variables).
-  ## column - abs_rate: absolute flow rate expression
+  ## column - abs_rate: absolute flow rate expression (per-cqpita flow rate 
+  ## multiplied by the size variable)
   ## example:
   ## size, change,    rate,          abs_rate
   ## S,    infection, beta * I / N,  S * (beta * I / N)
@@ -238,8 +241,8 @@ ChangeComponent = function() {
 ##' si$during()
 ##' 
 ##' @noRd
-SimpleChangeModel = function(before = list(), during = list(), after = list()) {
-  self = ChangeModelDefaults()
+SimpleChangeModel = function(before = list(), during = list(), after = list(), delta_t = NULL) {
+  self = ChangeModelDefaults(delta_t)
   
   self$before = before
   self$during = during
@@ -308,7 +311,7 @@ SimpleChangeModel = function(before = list(), during = list(), after = list()) {
 }
 
 AllFormulaChangeModel = function(before = list(), during = list(), after = list()) {
-  self = ChangeModelDefaults()
+  self = ChangeModelDefaults(delta_t = NULL)
   self$before = before
   self$during = during
   self$after = after
@@ -383,7 +386,11 @@ MockChangeModel = function() {
 ##' steps.
 ##' 
 ##' @param model Object with quantities that have been explicitly 
-##' marked as state variables.
+##' marked as state variables. Currently the only valid model object
+##' is a model specification (e.g., \code{\link{mp_tmb_model_spec}}).
+##' @param delta_t Number giving the amount of time that passes during a
+##' single time-step. The default, `NULL`, is equivalent to the value
+##' for `delta_t` that already exists in the `model`.
 ##' 
 ##' @examples
 ##' sir = mp_tmb_library("starter_models", "sir", package = "macpan2")
@@ -401,7 +408,7 @@ NULL
 ##' \code{\link{mp_tmb_model_spec}}, but this default can be changed using
 ##' the functions described below.
 ##' @export
-mp_euler = function(model) UseMethod("mp_euler")
+mp_euler = function(model, delta_t = NULL) UseMethod("mp_euler")
 
 ##' @describeIn state_updates ODE solver using Runge-Kutta 4. Any formulas that
 ##' appear before model flows in the `during` list will only be updated
@@ -424,19 +431,19 @@ mp_euler = function(model) UseMethod("mp_euler")
 ##' be confused. We therefore require that all state variable updates are set
 ##' explicitly (e.g., with \code{\link{mp_per_capita_flow}}).
 ##' @export
-mp_rk4 = function(model) UseMethod("mp_rk4")
+mp_rk4 = function(model, delta_t = NULL) UseMethod("mp_rk4")
 
 ##' @describeIn state_updates Old version of `mp_rk4` that doesn't keep track
 ##' of absolute flows through each time-step. As a result this version is
 ##' more efficient but makes it more difficult to compute things like 
 ##' incidence over a time scale.
 ##' @export
-mp_rk4_old = function(model) UseMethod("mp_rk4_old")
+mp_rk4_old = function(model, delta_t = NULL) UseMethod("mp_rk4_old")
 
 ##' @describeIn state_updates Original and deprecated name for 
 ##' `mp_discrete_stoch`. In all new projects please use `mp_discrete_stoch`.
 ##' @export
-mp_euler_multinomial = function(model) UseMethod("mp_euler_multinomial")
+mp_euler_multinomial = function(model, delta_t = NULL) UseMethod("mp_euler_multinomial")
 
 ##' @describeIn state_updates Update state such that the probability of moving
 ##' from box `i` to box `j` in one time step is given by
@@ -447,33 +454,31 @@ mp_euler_multinomial = function(model) UseMethod("mp_euler_multinomial")
 ##' distribution that determines how many individuals go to each `j` box and 
 ##' how many stay in  `i`.
 ##' @export
-mp_discrete_stoch = function(model) UseMethod("mp_discrete_stoch")
+mp_discrete_stoch = function(model, delta_t = NULL) UseMethod("mp_discrete_stoch")
 
 ##' @describeIn state_updates Update state with hazard steps, which is equivalent
 ##' to taking the step given by the expected value of the Euler-multinomial
 ##' distribution.
 ##' @export
-mp_hazard = function(model) UseMethod("mp_hazard")
+mp_hazard = function(model, delta_t =  NULL) UseMethod("mp_hazard")
 
 ##' @export
-mp_euler.TMBModelSpec = function(model) model$change_update_method("euler")
+mp_euler.TMBModelSpec = function(model, delta_t = NULL) model$change_update_method("euler", delta_t)
 
 ##' @export
-mp_rk4.TMBModelSpec = function(model) model$change_update_method("rk4")
+mp_rk4.TMBModelSpec = function(model, delta_t = NULL) model$change_update_method("rk4", delta_t)
 
 ##' @export
-mp_rk4_old.TMBModelSpec = function(model) model$change_update_method("rk4_old")
-
+mp_rk4_old.TMBModelSpec = function(model, delta_t = NULL) model$change_update_method("rk4_old", delta_t)
 
 ##' @export
-mp_euler_multinomial.TMBModelSpec = function(model) model$change_update_method("euler_multinomial")
+mp_euler_multinomial.TMBModelSpec = function(model, delta_t = NULL) model$change_update_method("euler_multinomial", delta_t)
 
 ##' @export
-mp_discrete_stoch.TMBModelSpec = function(model) model$change_update_method("discrete_stoch")
-
+mp_discrete_stoch.TMBModelSpec = function(model, delta_t = NULL) model$change_update_method("discrete_stoch", delta_t)
 
 ##' @export
-mp_hazard.TMBModelSpec = function(model) model$change_update_method("hazard")
+mp_hazard.TMBModelSpec = function(model, delta_t = NULL) model$change_update_method("hazard", delta_t)
 
 
 #' Expand Model
@@ -501,19 +506,21 @@ mp_expand.TMBModelSpec = function(model) model$expand()
 
 
 ## Utilities
-to_exogenous = function(flow_frame, rand_fn = NULL) {
+to_exogenous = function(flow_frame, rand_fn = NULL, dt = "") {
   frame = flow_frame[flow_frame$size == "", , drop = FALSE]
   if (is.null(rand_fn)) {
-    template = "%s ~ %s"
+    template = sprintf("%%s ~ %%s%s", dt)
   } else {
-    template = sprintf("%%s ~ %s(%%s)", rand_fn)
+    template = sprintf("%%s ~ %s(%%s%s)", rand_fn, dt)
   }
-   
   sprintf(template, frame$change, frame$abs_rate) |> lapply(as.formula)
 }
 to_exogenous_inputs = to_exogenous ## back-compat
-flow_frame_to_absolute_flows = function(flow_frame) {
-  char_vec = with(flow_frame, sprintf("%s ~ %s", change, abs_rate))
+flow_frame_to_absolute_flows = function(flow_frame, delta_t = 1) {
+  delta_t = handle_delta_t(delta_t)
+  delta_t_str = ""
+  if (delta_t != 1) delta_t_str = sprintf("%s * ", delta_t)
+  char_vec = with(flow_frame, sprintf("%s ~ %s%s", change, delta_t_str, abs_rate))
   lapply(char_vec, as.formula)
 }
 to_absolute_flows = function(per_capita_flows) {
@@ -541,13 +548,13 @@ get_state_update_method = function(state_update, change_model) {
   if (state_update == "rk4_old") cls_nm = "RK4OldUpdateMethod"
   get(cls_nm)(change_model)
 }
-get_change_model = function(before, during, after) {
+get_change_model = function(before, during, after, delta_t = NULL) {
   valid_before = all(vapply(before, is_two_sided, logical(1L)))
   if (!valid_before) stop("The before argument must be all two-sided formulas.")
   valid_after = all(vapply(after, is_two_sided, logical(1L)))
   if (!valid_after) stop("The after argument must be all two-sided formulas.")
   any_change_components = any(vapply(during, inherits, logical(1L), "ChangeComponent"))
-  if (any_change_components) return(SimpleChangeModel(before, during, after))
+  if (any_change_components) return(SimpleChangeModel(before, during, after, delta_t))
   AllFormulaChangeModel(before, during, after)
 }
 force_expr_list = function(x) {
@@ -596,7 +603,10 @@ EulerUpdateMethod = function(change_model, existing_global_names = character())
   self$before = function() self$change_model$before_loop()
   self$during = function() {
     before_components = self$change_model$before_flows()
-    components = flow_frame_to_absolute_flows(self$change_model$flow_frame())
+    components = flow_frame_to_absolute_flows(
+        self$change_model$flow_frame()
+      , self$change_model$delta_t
+    )
     before_update = self$change_model$before_state()
     update = self$change_model$update_state()
     
@@ -685,9 +695,10 @@ RK4UpdateMethod = function(change_model) {
   
   self$before = function() self$change_model$before_loop()
   self$during = function() {
-    ## abbr func nms
     unlst = function(x) unlist(x, recursive = FALSE, use.names = FALSE)
     as_forms = function(x) lapply(x, as.formula)
+    dt = self$change_model$delta_t
+    dt = if (dt == 1) "" else sprintf(" * %s", dt)
     
     before_components = self$change_model$before_flows()
     flow_frame = self$change_model$flow_frame()
@@ -730,33 +741,57 @@ RK4UpdateMethod = function(change_model) {
     }
     
     ## rk4 step 1
-    flow_replacements = sprintf("%s ~ %s", flows, flow_step_names$k1) |> as_forms()
+    flow_replacements = sprintf("%s ~ %s"
+      , flows, flow_step_names$k1
+    ) |> as_forms()
     k1_new_before = make_before("k1")
     k1_new_update = update_formulas(rate_formulas, flow_replacements)
     
     ## rk4 step 2
-    state_replacements = sprintf("%s ~ (%s + (%s / 2))", states, states, state_step_names$k1) |> as_forms()
-    flow_replacements = sprintf("%s ~ %s", flows, flow_step_names$k2) |> as_forms()
+    state_replacements = sprintf("%s ~ (%s + (%s%s / 2))"
+      , states, states, state_step_names$k1, dt
+    ) |> as_forms()
+    flow_replacements = sprintf("%s ~ %s"
+      , flows, flow_step_names$k2
+    ) |> as_forms()
     k2_new_before = update_formulas(make_before("k2"), state_replacements)
-    k2_new_update = sprintf("%s ~ %s", state_step_names$k2, rates) |> as_forms() |> update_formulas(flow_replacements)
+    k2_new_update = sprintf("%s ~ %s"
+      , state_step_names$k2, rates
+    ) |> as_forms() |> update_formulas(flow_replacements)
     
     ## rk4 step 3
-    state_replacements = sprintf("%s ~ (%s + (%s / 2))", states, states, state_step_names$k2) |> as_forms()
-    flow_replacements = sprintf("%s ~ %s", flows, flow_step_names$k3) |> as_forms()
+    state_replacements = sprintf("%s ~ (%s + (%s%s / 2))"
+      , states, states, state_step_names$k2, dt
+    ) |> as_forms()
+    flow_replacements = sprintf("%s ~ %s"
+      , flows, flow_step_names$k3
+    ) |> as_forms()
     k3_new_before = update_formulas(make_before("k3"), state_replacements)
-    k3_new_update = sprintf("%s ~ %s", state_step_names$k3, rates) |> as_forms() |> update_formulas(flow_replacements)
+    k3_new_update = sprintf("%s ~ %s"
+      , state_step_names$k3, rates
+    ) |> as_forms() |> update_formulas(flow_replacements)
     
     ## rk4 step 4
-    state_replacements = sprintf("%s ~ (%s + %s)", states, states, state_step_names$k3) |> as_forms()
-    flow_replacements = sprintf("%s ~ %s", flows, flow_step_names$k4) |> as_forms()
+    state_replacements = sprintf("%s ~ (%s + %s%s)"
+      , states, states, state_step_names$k3, dt
+    ) |> as_forms()
+    flow_replacements = sprintf("%s ~ %s"
+      , flows, flow_step_names$k4
+    ) |> as_forms()
     k4_new_before = update_formulas(make_before("k4"), state_replacements)
-    k4_new_update = sprintf("%s ~ %s", state_step_names$k4, rates) |> as_forms() |> update_formulas(flow_replacements)
+    k4_new_update = sprintf("%s ~ %s"
+      , state_step_names$k4, rates
+    ) |> as_forms() |> update_formulas(flow_replacements)
     
     ## final update step
-    final_flow_update = sprintf("%s ~ (%s + 2 * %s + 2 * %s + %s)/6"
-      , flows, flow_step_names$k1, flow_step_names$k2, flow_step_names$k3, flow_step_names$k4
+    final_flow_update = sprintf("%s ~ (%s + 2 * %s + 2 * %s + %s)%s/6"
+      , flows
+      , flow_step_names$k1, flow_step_names$k2
+      , flow_step_names$k3, flow_step_names$k4, dt
     ) |> as_forms()
-    final_state_update = sprintf("%s ~ %s %s", states, states, rates) |> as_forms() |> setNames(states)
+    final_state_update = sprintf("%s ~ %s %s"
+      , states, states, rates
+    ) |> as_forms() |> setNames(states)
     after_components = self$change_model$after_state()
     c(
         k1_new_before, k1_new_update
@@ -779,12 +814,20 @@ DiscreteStochUpdateMethod = function(change_model) {
 EulerMultinomialUpdateMethod = function(change_model) {
   self = Base()
   self$change_model = change_model
+  dt = self$change_model$delta_t
   
   self$vec = function(expr_list, char_fun) {
     vec = vapply(expr_list, char_fun, character(1L))
+
+    # - simple expressions are any non-formula strings 
+    #   (names of variables or state flows)
+    # - expressions that are not simple contain math symbols 
+    #   (e.g., +,-,*,/, etc.)
+
     ## simple expressions are any non-formula strings (names of variables or state flows)
     ## expressions that are not simple contain math symbols (ex. +,-,*,/, etc.)
     ## BMB: allow non-ASCII alpha? "^[[:alpha:]0-9._]+$"
+
     simple_expr = all(grepl("^[a-zA-Z0-9._]+$", vec))
     scalar_expr = length(vec) == 1L
     
@@ -800,17 +843,22 @@ EulerMultinomialUpdateMethod = function(change_model) {
   
   self$before = function() self$change_model$before_loop()
   self$during = function() {
+    dt = self$change_model$delta_t
+    dt_emult = if (dt == 1) "" else sprintf(", %s", dt)
+    dt_pois = if (dt == 1) "" else sprintf(" * %s", dt)
+    
     before_components = c(
         self$change_model$before_flows()
-      , to_exogenous(self$change_model$flow_frame(), rand_fn = "rpois")
+      , to_exogenous(self$change_model$flow_frame(), rand_fn = "rpois", dt = dt_pois)
     )
     flow_list = self$change_model$update_flows()
     components = list()
     for (size_var in names(flow_list)) {
-      components[[size_var]] = sprintf("%s ~ reulermultinom(%s, %s)"
+      components[[size_var]] = sprintf("%s ~ reulermultinom(%s, %s%s)"
         , self$vec(flow_list[[size_var]], lhs_char)
         , size_var
         , self$vec(flow_list[[size_var]], rhs_char)
+        , dt_emult
       )
     }
     ## BMB: add absolute flows here
@@ -833,18 +881,23 @@ EulerMultinomialUpdateMethod = function(change_model) {
 HazardUpdateMethod = function(change_model) {
   self = EulerMultinomialUpdateMethod(change_model)
   self$during = function() {
+    dt = self$change_model$delta_t
+    dt_flow = if (dt == 1) "" else sprintf("%s * ", dt)
+    dt_exog = if (dt == 1) "" else sprintf("/%s", dt)
     before_components = c(
         self$change_model$before_flows()
-      , to_exogenous(self$change_model$flow_frame())
+      , to_exogenous(self$change_model$flow_frame(), dt = dt_exog)
     )
     before_state = self$change_model$before_state()
 
     flow_list = self$change_model$update_flows()
     components = list()
     for (size_var in names(flow_list)) {
-      components[[size_var]] = sprintf("%s ~ %s * (1 - exp(-sum(%s))) * proportions(%s, 0, %s)"
+      components[[size_var]] = sprintf(
+          "%s ~ %s * (1 - exp(-%ssum(%s))) * proportions(%s, 0, %s)"
         , self$vec(flow_list[[size_var]], lhs_char)
         , size_var
+        , dt_flow
         , self$vec(flow_list[[size_var]], rhs_char)
         , self$vec(flow_list[[size_var]], rhs_char)
         , getOption("macpan2_tol_hazard_div")
@@ -883,7 +936,7 @@ HazardUpdateMethod = function(change_model) {
 #' originates.
 #' @param to String giving the name of the compartment to which the flow is
 #' going.
-#' @param rate String giving the expression for the per-capita
+#' @param rate String giving the expression for the per-capita or absolute
 #' flow rate. Alternatively, and for back compatibility, 
 #' a two-sided formula with the left-hand-side giving the name of the absolute 
 #' flow rate per time-step and the right-hand-side giving an expression for 
@@ -931,7 +984,7 @@ HazardUpdateMethod = function(change_model) {
 #' mp_per_capita_flow("S", "V", "(vrate * S/(vrate + S))/S", "vaccination")
 #' 
 #' # importation
-#' # mp_inflow("I", "delta", "importation")
+#' mp_inflow("I", "delta", "importation")
 #' 
 #' @export
 mp_per_capita_flow = function(from, to, rate, flow_name = NULL, abs_rate = NULL) {
@@ -965,6 +1018,18 @@ mp_per_capita_outflow = function(from, rate, flow_name = NULL, abs_rate = NULL)
   PerCapitaOutflow(from, rate, call_string)
 }
 
+#' @describeIn mp_per_capita_flow Alternative to `mp_per_capita_flow` that
+#' allowing specification of flows using absolute rates instead of per-capita
+#' rates.
+#' @param rate_name Deprecated synonym for `flow_name` when using 
+#' `mp_absolute_flow`. Please use `flow_name` in all future work.
+#' @export
+mp_absolute_flow = function(from, to, rate, flow_name = NULL, rate_name = NULL) {
+  call_string = deparse(match.call())
+  rate = handle_abs_rate_args(rate, rate_name, flow_name)
+  AbsoluteFlow(from, to, rate, call_string)
+}
+
 #' @describeIn mp_per_capita_flow Only flow into the `to` compartment.
 #' For adding a birth or immigration process.
 #' @param flow_name String giving the name of the flow
@@ -991,34 +1056,6 @@ mp_outflow = function(from, rate, flow_name = NULL, abs_rate = NULL) {
 }
 
 
-
-
-#' Specify Absolute Flow Between Compartments (Experimental)
-#' 
-#' An experimental alternative to \code{\link{mp_per_capita_flow}} that 
-#' allows users to specify flows using absolute rates instead of 
-#' per-capita rates.
-#' 
-#' @param from String giving the name of the compartment from which the flow
-#' originates.
-#' @param to String giving the name of the compartment to which the flow is
-#' going.
-#' @param rate String giving the expression for the absolute
-#' flow rate per time-step.
-#' @param flow_name String giving the name for the variable that 
-#' will store the `rate`.
-#' @param rate_name Deprecated synonym for `flow_name`. Please use `flow_name`
-#' in all future work.
-#' 
-#' @seealso [mp_per_capita_flow()]
-#' 
-#' @export
-mp_absolute_flow = function(from, to, rate, flow_name = NULL, rate_name = NULL) {
-  call_string = deparse(match.call())
-  rate = handle_abs_rate_args(rate, rate_name, flow_name)
-  AbsoluteFlow(from, to, rate, call_string)
-}
-
 PerCapitaOutflow = function(from, rate, call_string) {
   self = PerCapitaFlow(from, NULL, rate, call_string)
   self$change_frame = function() {
@@ -1080,18 +1117,8 @@ AbsoluteFlow = function(from, to, rate, call_string) {
   self$flow_frame = function() {
     abs_rate = rhs_char(self$rate)
     data.frame(
-        ## BMB: not sure if this is right? there is no 'size'
-        size = "" ## self$from %||% ""
+        size = "" ## absolute flows are coming from 'nowhere' and so there is no 'size'
       , change = lhs_char(self$rate)
-      
-      ## this is the main problem with absolute flows, because it has a 
-      ## `size` variable (e.g., a `from` variable) in the denominator.
-      ## this issue should only arise for update methods that are more 
-      ## naturally expressed for per-capita flows (e.g., Euler-multinomial
-      ## and hazard)
-      ## BMB: why does there need to be a 'from' involved here at all?
-      ##  is this only an issue because we need a Poisson-type stochastic
-      ##  flow to go with the Euler-multinomial?
       , rate = sprintf("%s", abs_rate)
       , abs_rate = abs_rate
     )
diff --git a/R/mp_tmb_calibrator.R b/R/mp_tmb_calibrator.R
index 16177861..12b625ad 100644
--- a/R/mp_tmb_calibrator.R
+++ b/R/mp_tmb_calibrator.R
@@ -85,6 +85,11 @@
 #' if (suppressPackageStartupMessages(require(broom.mixed))) {
 #'   print(mp_tmb_coef(cal))
 #' }
+#' (cal 
+#'   |> mp_optimized_spec() 
+#'   |> mp_default_list() 
+#'   |> getElement("beta")
+#' )
 #' @concept create-model-calibrator
 #' @export
 mp_tmb_calibrator = function(spec
@@ -713,7 +718,8 @@ mp_optimizer_output.TMBCalibrator = function(model, what = c("latest", "all")) {
 #' Optimized Model Specification
 #' 
 #' Create a new model specification using parameter values that have been
-#' optimized.
+#' optimized. See the examples in the help page for 
+#' \code{\link{mp_tmb_calibrator}}.
 #' 
 #' @param model A model object that can be optimized/calibrated. Currently,
 #' only models produced by \code{\link{mp_tmb_calibrator}} are valid.
diff --git a/R/mp_tmb_model_spec.R b/R/mp_tmb_model_spec.R
index 280becdf..45fd98b2 100644
--- a/R/mp_tmb_model_spec.R
+++ b/R/mp_tmb_model_spec.R
@@ -16,6 +16,7 @@ TMBModelSpec = function(
         , "rk4_old"
         , "euler_multinomial"
       )
+    , delta_t = NULL
   ) {
   default = c(default, inits)
   must_not_save = handle_saving_conflicts(must_save, must_not_save)
@@ -24,7 +25,8 @@ TMBModelSpec = function(
   before = force_expr_list(before)
   during = force_expr_list(during)
   after = force_expr_list(after)
-  self$change_model = get_change_model(before, during, after)
+  self$delta_t = handle_delta_t(delta_t)
+  self$change_model = get_change_model(before, during, after, self$delta_t)
   self$state_update = get_state_update_type(match.arg(state_update), self$change_model)
   self$update_method = get_state_update_method(self$state_update, self$change_model)
   self$change_components = function() self$change_model$change_list
@@ -119,10 +121,15 @@ TMBModelSpec = function(
       , must_save = self$must_save, must_not_save = self$must_not_save
       , sim_exprs = self$sim_exprs
       , state_update = self$state_update
+      , delta_t = self$delta_t
     )
   }
   self$change_update_method = function(
-      state_update = c("euler", "rk4", "discrete_stoch", "hazard", "rk4_old", "euler_multinomial")
+      state_update = c(
+          "euler", "rk4", "discrete_stoch", "hazard"
+        , "rk4_old", "euler_multinomial"
+      ),
+      delta_t = NULL
     ) {
     
     if (self$state_update == "no") {
@@ -137,7 +144,9 @@ TMBModelSpec = function(
         before = self$before, during = self$during, after = self$after
       , default = self$default, integers = self$integers
       , must_save = self$must_save, must_not_save = self$must_not_save
-      , sim_exprs = self$sim_exprs, state_update = state_update
+      , sim_exprs = self$sim_exprs
+      , state_update = state_update
+      , delta_t = handle_delta_t(delta_t, self$delta_t)
     )
   }
   self$expand = function() {
@@ -151,6 +160,7 @@ TMBModelSpec = function(
       , must_not_save = self$must_not_save
       , sim_exprs = self$sim_exprs
       , state_update = self$state_update
+      , delta_t = self$delta_t
     )
   }
   self$name_map = function(local_names) {
@@ -205,6 +215,12 @@ TMBModelSpec = function(
   return_object(self, "TMBModelSpec")
 }
 
+handle_delta_t = function(delta_t_new = NULL, delta_t_old = NULL) {
+  if (!is.null(delta_t_new)) return(delta_t_new)
+  if (!is.null(delta_t_old)) return(delta_t_old)
+  return(1)
+}
+
 handle_saving_conflicts = function(must_save, must_not_save) {
   ## must_save takes precedence over must_not_save
   problems = intersect(must_save, must_not_save)
@@ -329,6 +345,8 @@ must_save_time_args = function(formulas) {
 #' @param state_update Optional character vector for how to update the state 
 #' variables when it is relevant. Options include `"euler"` (the default), 
 #' `"rk4"`, and `"discrete_stoch"`.
+#' @param delta_t Number giving the amount of time that passes during a
+#' single time-step.
 #' 
 #' @details
 #' Expressions in the `before`, `during`, and `after` lists can be standard 
@@ -435,6 +453,7 @@ mp_print_before = function(model) {
       model$before
     , c(length(model$before), 0L, 0L)
   )
+  invisible(model$before)
 }
 
 #' @describeIn mp_print_spec Print just the expressions executed during each
@@ -445,6 +464,7 @@ mp_print_during = function(model) {
       model$during
     , c(0L, length(model$during), 0L)
   )
+  invisible(model$during)
 }
 
 #' @describeIn mp_print_spec Print just the expressions executed after the
@@ -455,6 +475,7 @@ mp_print_after = function(model) {
       model$after
     , c(0L, 0L, length(model$after))
   )
+  invisible(model$after)
 }
 
 #' Version Update
diff --git a/R/tmb_model.R b/R/tmb_model.R
index 8cb571cf..8ac00e5f 100644
--- a/R/tmb_model.R
+++ b/R/tmb_model.R
@@ -895,10 +895,11 @@ mp_trajectory_sim.TMBSimulator = function(model
   r = r[, names(r) != "value", drop = FALSE]
   rr = (n
     |> replicate(model$simulate_values()) 
-    |> apply(1, quantile, probs, na.rm = TRUE)
+    |> apply(1, quantile, probs, na.rm = TRUE, simplify = FALSE)
+    |> simplify2array(except = 0L)
     |> t()
   )
-  names(rr) = sprintf("value_%s", names(rr))
+  colnames(rr) = sprintf("value_%s", colnames(rr))
   cbind(r, rr)
 }
 
diff --git a/R/tmb_model_editors.R b/R/tmb_model_editors.R
index 41cf1fb3..465bf45b 100644
--- a/R/tmb_model_editors.R
+++ b/R/tmb_model_editors.R
@@ -112,6 +112,7 @@ mp_tmb_insert = function(model
     , must_not_save = model$must_not_save
     , sim_exprs = model$sim_exprs
     , state_update = model$state_update
+    , delta_t = model$delta_t
   )
 }
 
@@ -155,6 +156,7 @@ mp_tmb_update = function(model
     , must_not_save = model$must_not_save
     , sim_exprs = model$sim_exprs
     , state_update = model$state_update
+    , delta_t = model$delta_t
   )
 }
 
@@ -193,6 +195,7 @@ mp_tmb_delete = function(model
     , must_not_save = model$must_not_save
     , sim_exprs = model$sim_exprs
     , state_update = model$state_update
+    , delta_t = model$delta_t
   )
 }
 
diff --git a/inst/starter_models/one_box/README.Rmd b/inst/starter_models/one_box/README.Rmd
new file mode 100644
index 00000000..039a04b6
--- /dev/null
+++ b/inst/starter_models/one_box/README.Rmd
@@ -0,0 +1,83 @@
+---
+title: "One Box"
+index_entry: "A compartmental model with a single compartment"
+bibliography: ../../references.bib
+link-citations: TRUE
+author: Steve Walker
+output: 
+  github_document:
+    toc: true
+editor_options: 
+  chunk_output_type: console
+---
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>",
+  fig.path = "./figures/"
+)
+```
+
+This is a very simple model of importation with exponential decay due to mortality.
+
+The code in this article uses the following packages.
+
+```{r packages, message=FALSE, warning=FALSE}
+library(macpan2)
+library(ggplot2)
+library(dplyr)
+```
+
+# States
+
+| variable | description           |
+| -------- | --------------------- |
+| N        | Number of individuals |
+
+# Parameters
+
+| variable | description                  |
+| -------- | ---------------------------- |
+| $r$      | absolute rate of importation |
+| $d$      | per capita death rate        |
+
+# Dynamics 
+
+If we interpret this model as an ODE, we have the following dynamics.
+
+$$
+\frac{dN}{dt} = r - (d \times N)  \\
+$$
+
+
+# Model Specification
+
+This model has been specified in the `one_box` directory [here](https://github.com/canmod/macpan2/blob/main/inst/starter_models/one_box/tmb.R) and is accessible from the `macpan2` model library (see [Example Models](https://canmod.github.io/macpan2/articles/example_models.html) for details). 
+
+# Simulation
+
+Simulation of this and other models depends on the kind of [state update](https://canmod.github.io/macpan2/reference/mp_euler) that you use. Here use the discrete-stochastic method.
+
+```{r simulation}
+spec = mp_tmb_library("starter_models"
+  , "one_box"
+  , package = "macpan2"
+)
+set.seed(1)
+(spec
+  |> mp_rk4(delta_t = 0.1)
+  |> mp_simulator(50 / 0.1, "N")
+  |> mp_trajectory(include_initial = TRUE)
+  |> rename(prevalence = value)
+  |> ggplot()
+  + geom_line(aes(time, prevalence))
+  
+  ## asymptotic equilibrium population size
+  + geom_hline(yintercept = with(spec$default, r / d)
+      , colour = "red"
+      , linetype = "dashed"
+    )
+  + theme_bw()
+)
+```
diff --git a/inst/starter_models/one_box/README.md b/inst/starter_models/one_box/README.md
new file mode 100644
index 00000000..161e8c7c
--- /dev/null
+++ b/inst/starter_models/one_box/README.md
@@ -0,0 +1,78 @@
+One Box
+================
+Steve Walker
+
+- [States](#states)
+- [Parameters](#parameters)
+- [Dynamics](#dynamics)
+- [Model Specification](#model-specification)
+- [Simulation](#simulation)
+
+This is a very simple model of importation with exponential decay due to
+mortality.
+
+The code in this article uses the following packages.
+
+``` r
+library(macpan2)
+library(ggplot2)
+library(dplyr)
+```
+
+# States
+
+| variable | description           |
+|----------|-----------------------|
+| N        | Number of individuals |
+
+# Parameters
+
+| variable | description                  |
+|----------|------------------------------|
+| $r$      | absolute rate of importation |
+| $d$      | per capita death rate        |
+
+# Dynamics
+
+If we interpret this model as an ODE, we have the following dynamics.
+
+$$
+\frac{dN}{dt} = r - (d \times N)  \\
+$$
+
+We can also interpret this model as a discrete-time stochastic model, as
+we show below in the [Simulation](#simulation) section.
+
+# Model Specification
+
+This model has been specified in the `one_box` directory
+[here](https://github.com/canmod/macpan2/blob/main/inst/starter_models/one_box/tmb.R)
+and is accessible from the `macpan2` model library (see [Example
+Models](https://canmod.github.io/macpan2/articles/example_models.html)
+for details).
+
+# Simulation
+
+Simulation of this and other models depends on the kind of [state
+update](https://canmod.github.io/macpan2/reference/mp_euler) that you
+use. Here use the discrete-stochastic method.
+
+``` r
+specs = mp_tmb_library("starter_models"
+  , "one_box"
+  , package = "macpan2"
+)
+set.seed(1)
+(specs
+  |> mp_discrete_stoch()
+  |> mp_simulator(75, "N")
+  |> mp_trajectory_replicate(n = 100)
+  |> bind_rows(.id = "replicate")
+  |> rename(prevalence = value)
+  |> ggplot()
+  + geom_line(aes(time, prevalence, group = replicate), alpha = 0.2)
+  + theme_bw()
+)
+```
+
+![](./figures/simulation-1.png)<!-- -->
diff --git a/inst/starter_models/one_box/figures/simulation-1.png b/inst/starter_models/one_box/figures/simulation-1.png
new file mode 100644
index 00000000..96efb898
Binary files /dev/null and b/inst/starter_models/one_box/figures/simulation-1.png differ
diff --git a/inst/starter_models/one_box/tmb.R b/inst/starter_models/one_box/tmb.R
new file mode 100644
index 00000000..d6676c6d
--- /dev/null
+++ b/inst/starter_models/one_box/tmb.R
@@ -0,0 +1,8 @@
+library(macpan2)
+spec = mp_tmb_model_spec(
+    during = list(
+          mp_inflow("N", "r", "importation")
+        , mp_per_capita_outflow("N", "d", "death")
+    )
+  , default = list(N = 0, r = 0.3, d = 0.2)
+)
diff --git a/man/mp_absolute_flow.Rd b/man/mp_absolute_flow.Rd
deleted file mode 100644
index 445c2715..00000000
--- a/man/mp_absolute_flow.Rd
+++ /dev/null
@@ -1,32 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/formula_list_generators.R
-\name{mp_absolute_flow}
-\alias{mp_absolute_flow}
-\title{Specify Absolute Flow Between Compartments (Experimental)}
-\usage{
-mp_absolute_flow(from, to, rate, flow_name = NULL, rate_name = NULL)
-}
-\arguments{
-\item{from}{String giving the name of the compartment from which the flow
-originates.}
-
-\item{to}{String giving the name of the compartment to which the flow is
-going.}
-
-\item{rate}{String giving the expression for the absolute
-flow rate per time-step.}
-
-\item{flow_name}{String giving the name for the variable that
-will store the \code{rate}.}
-
-\item{rate_name}{Deprecated synonym for \code{flow_name}. Please use \code{flow_name}
-in all future work.}
-}
-\description{
-An experimental alternative to \code{\link{mp_per_capita_flow}} that
-allows users to specify flows using absolute rates instead of
-per-capita rates.
-}
-\seealso{
-\code{\link[=mp_per_capita_flow]{mp_per_capita_flow()}}
-}
diff --git a/man/mp_optimized_spec.Rd b/man/mp_optimized_spec.Rd
index 49e64924..c119c8c3 100644
--- a/man/mp_optimized_spec.Rd
+++ b/man/mp_optimized_spec.Rd
@@ -22,7 +22,8 @@ model, with calibrated default values.
 }
 \description{
 Create a new model specification using parameter values that have been
-optimized.
+optimized. See the examples in the help page for
+\code{\link{mp_tmb_calibrator}}.
 }
 \details{
 The following is a list of additional model specification structure that was
diff --git a/man/mp_per_capita_flow.Rd b/man/mp_per_capita_flow.Rd
index e394b820..ec4aef32 100644
--- a/man/mp_per_capita_flow.Rd
+++ b/man/mp_per_capita_flow.Rd
@@ -4,6 +4,7 @@
 \alias{mp_per_capita_flow}
 \alias{mp_per_capita_inflow}
 \alias{mp_per_capita_outflow}
+\alias{mp_absolute_flow}
 \alias{mp_inflow}
 \alias{mp_outflow}
 \title{Specify Flow Into, Out Of, and Between Compartments}
@@ -14,6 +15,8 @@ mp_per_capita_inflow(from, to, rate, flow_name = NULL, abs_rate = NULL)
 
 mp_per_capita_outflow(from, rate, flow_name = NULL, abs_rate = NULL)
 
+mp_absolute_flow(from, to, rate, flow_name = NULL, rate_name = NULL)
+
 mp_inflow(to, rate, flow_name = NULL, abs_rate = NULL)
 
 mp_outflow(from, rate, flow_name = NULL, abs_rate = NULL)
@@ -25,7 +28,7 @@ originates.}
 \item{to}{String giving the name of the compartment to which the flow is
 going.}
 
-\item{rate}{String giving the expression for the per-capita
+\item{rate}{String giving the expression for the per-capita or absolute
 flow rate. Alternatively, and for back compatibility,
 a two-sided formula with the left-hand-side giving the name of the absolute
 flow rate per time-step and the right-hand-side giving an expression for
@@ -35,6 +38,9 @@ the per-capita rate of flow from \code{from} to \code{to}.}
 
 \item{abs_rate}{Deprecated synonym for \code{flow_name}. Please use \code{flow_name}
 in all future work.}
+
+\item{rate_name}{Deprecated synonym for \code{flow_name} when using
+\code{mp_absolute_flow}. Please use \code{flow_name} in all future work.}
 }
 \description{
 Specify different kinds of flows into, out of, and between compartments.
@@ -62,6 +68,10 @@ system (e.g., death). To keep track of the total number of dead individuals
 one can use \code{mp_per_capita_flow} and set \code{to} to be a compartment for
 these individuals (e.g., \code{to = "D"}).
 
+\item \code{mp_absolute_flow()}: Alternative to \code{mp_per_capita_flow} that
+allowing specification of flows using absolute rates instead of per-capita
+rates.
+
 \item \code{mp_inflow()}: Only flow into the \code{to} compartment.
 For adding a birth or immigration process.
 
@@ -105,7 +115,7 @@ mp_per_capita_outflow("R", "mu", "death_R")
 mp_per_capita_flow("S", "V", "(vrate * S/(vrate + S))/S", "vaccination")
 
 # importation
-# mp_inflow("I", "delta", "importation")
+mp_inflow("I", "delta", "importation")
 
 }
 \seealso{
diff --git a/man/mp_tmb_calibrator.Rd b/man/mp_tmb_calibrator.Rd
index 56f495bc..29a3fcbe 100644
--- a/man/mp_tmb_calibrator.Rd
+++ b/man/mp_tmb_calibrator.Rd
@@ -115,5 +115,10 @@ mp_optimize(cal)
 if (suppressPackageStartupMessages(require(broom.mixed))) {
   print(mp_tmb_coef(cal))
 }
+(cal 
+  |> mp_optimized_spec() 
+  |> mp_default_list() 
+  |> getElement("beta")
+)
 }
 \concept{create-model-calibrator}
diff --git a/man/mp_tmb_model_spec.Rd b/man/mp_tmb_model_spec.Rd
index 01b311b1..f99328de 100644
--- a/man/mp_tmb_model_spec.Rd
+++ b/man/mp_tmb_model_spec.Rd
@@ -15,7 +15,8 @@ mp_tmb_model_spec(
   must_not_save = character(),
   sim_exprs = character(),
   state_update = c("euler", "rk4", "discrete_stoch", "hazard", "rk4_old",
-    "euler_multinomial")
+    "euler_multinomial"),
+  delta_t = NULL
 )
 }
 \arguments{
@@ -72,6 +73,9 @@ continuous.}
 \item{state_update}{Optional character vector for how to update the state
 variables when it is relevant. Options include \code{"euler"} (the default),
 \code{"rk4"}, and \code{"discrete_stoch"}.}
+
+\item{delta_t}{Number giving the amount of time that passes during a
+single time-step.}
 }
 \description{
 Specify a simulation model in the TMB engine. A detailed explanation of this
diff --git a/man/state_updates.Rd b/man/state_updates.Rd
index 64097a97..8aba5bd6 100644
--- a/man/state_updates.Rd
+++ b/man/state_updates.Rd
@@ -10,21 +10,26 @@
 \alias{mp_hazard}
 \title{Change How State Variables are Updated}
 \usage{
-mp_euler(model)
+mp_euler(model, delta_t = NULL)
 
-mp_rk4(model)
+mp_rk4(model, delta_t = NULL)
 
-mp_rk4_old(model)
+mp_rk4_old(model, delta_t = NULL)
 
-mp_euler_multinomial(model)
+mp_euler_multinomial(model, delta_t = NULL)
 
-mp_discrete_stoch(model)
+mp_discrete_stoch(model, delta_t = NULL)
 
-mp_hazard(model)
+mp_hazard(model, delta_t = NULL)
 }
 \arguments{
 \item{model}{Object with quantities that have been explicitly
-marked as state variables.}
+marked as state variables. Currently the only valid model object
+is a model specification (e.g., \code{\link{mp_tmb_model_spec}}).}
+
+\item{delta_t}{Number giving the amount of time that passes during a
+single time-step. The default, \code{NULL}, is equivalent to the value
+for \code{delta_t} that already exists in the \code{model}.}
 }
 \description{
 These functions return a modified version of a model specification, such
diff --git a/tests/testthat/test-absolute.R b/tests/testthat/test-absolute.R
index 060aea6c..19ae8690 100644
--- a/tests/testthat/test-absolute.R
+++ b/tests/testthat/test-absolute.R
@@ -1,15 +1,9 @@
 test_that("absolute inflow definition and ode trajectory", {
     ## exponential decay
-    s1 <- mp_tmb_model_spec(
-          during = list(
-                mp_inflow("N", "r", "birth")
-              , mp_per_capita_outflow("N", "d", "death")
-          )
-        , default = list(N = 0, r = 0.3, d = 0.2)
-    )
+    s1 = test_cache_read("SPEC-one_box.rds")
     
     expect_equal(mp_state_vars(s1), "N")
-    expect_equal(mp_flow_vars(s1), c("birth", "death"))
+    expect_equal(mp_flow_vars(s1), c("importation", "death"))
     
     steps = 50
     
@@ -81,5 +75,18 @@ test_that("equivalent absolute and per-capita flows are consistent", {
   }
   expect_equal(sim_fn(sir_pc), sim_fn(sir_ab))
   expect_equal(sim_fn(sir_pc, mp_rk4), sim_fn(sir_ab, mp_rk4))
-  ## expect_equal(sim_fn(sir_pc, mp_hazard), sim_fn(sir_ab2, mp_hazard)) ## -- failing because hazard is behaving like euler with absolute flows
 })
+
+sir = mp_tmb_model_spec(
+    during = list(
+        mp_per_capita_flow("S", "I", "beta * I / N", "infection")
+      , mp_inflow("R", "gamma", "recovery")
+      , mp_absolute_flow("I", "R", "2 * gamma", "recovery2")
+    )
+  , default = list(
+        beta = 0.2, gamma = 0.1
+      , S = 99, I = 1, R = 0, N = 100
+    )
+)
+sir |> mp_expand()
+sir |> mp_discrete_stoch() |> mp_expand()
diff --git a/tests/testthat/test-change-model.R b/tests/testthat/test-change-model.R
index ec11d26c..33bf5d5b 100644
--- a/tests/testthat/test-change-model.R
+++ b/tests/testthat/test-change-model.R
@@ -70,5 +70,3 @@ em = (macpan_base
   |> mp_trajectory()
 )
 el = mp_simulator(macpan_base, 20L, "S.E") |> mp_trajectory()
-
-
diff --git a/tests/testthat/test-delta-t.R b/tests/testthat/test-delta-t.R
new file mode 100644
index 00000000..e484a807
--- /dev/null
+++ b/tests/testthat/test-delta-t.R
@@ -0,0 +1,66 @@
+if (interactive()) {
+  library(macpan2); library(testthat); library(deSolve)
+  source("tests/testthat/setup.R")
+}
+test_that("rk4 time-steps work with the si model", {
+  delta_t = 0.1
+  time = 50
+  si = test_cache_read("SPEC-si.rds") |> mp_rk4(delta_t)
+  solution = si |> mp_simulator(time/delta_t, "I") |> mp_trajectory()
+  state = si$default["I"]
+  params = si$default[c("beta", "N")]
+  reference_solution <- ode(
+      y = unlist(state)
+    , times = seq(from = 0, to = time, by = delta_t)
+    , func = \(t, state, params) with(c(params, state), list((N - I) * (beta * I/N)))
+    , parms = params
+    , method = "rk4"
+  )
+  expect_equal(solution$value, reference_solution[-1, "I"])
+})
+
+test_that("rk4 time-steps work with absolute inflows", {
+  delta_t = 0.1
+  time = 50
+  one_box = test_cache_read("SPEC-one_box.rds") |> mp_rk4(delta_t)
+  solution = one_box |> mp_simulator(time/delta_t, "N") |> mp_trajectory()
+  state = one_box$default["N"]
+  params = one_box$default[c("r", "d")]
+  reference_solution <- ode(
+      y = unlist(state)
+    , times = seq(from = 0, to = time, by = delta_t)
+    , func = \(t, state, params) with(c(params, state), list(r - d * N))
+    , parms = params
+    , method = "rk4"
+  )
+  expect_equal(solution$value, reference_solution[-1, "N"])
+})
+
+test_that("stochastic absolute importation with a time step is plausible", {
+  spec = mp_tmb_model_spec(
+      during = list(
+          N ~ S + I
+        , mp_per_capita_flow("S", "I", "beta * I / N", "infection")
+        , mp_inflow("I", "mu", "importation")
+      )
+    , default = list(S = 100, I = 0, beta = 0.8, mu = 0.1)
+  )
+  time = 10
+  delta_t = 0.1
+  set.seed(1L)
+  traj = (spec
+    |> mp_discrete_stoch(delta_t)
+    |> mp_simulator(time / delta_t, "I") 
+    |> mp_trajectory_sim(10, probs = 0.5)
+    |> mutate(time = delta_t * time)
+    |> filter(time == round(time))
+  )
+  row.names(traj) = 1:10 ## seems this is needed for testthat -- not sure why
+  answer = structure(list(matrix = c("I", "I", "I", "I", "I", "I", "I", 
+  "I", "I", "I"), time = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10), row = c(0, 
+  0, 0, 0, 0, 0, 0, 0, 0, 0), col = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 
+  0), `value_50%` = c(0, 0, 0, 0, 0, 0, 0, 0, 2, 3)), class = "data.frame", 
+  row.names = 1:10)
+  
+  expect_equal(traj, answer)
+})
diff --git a/vignettes/FAQs.Rmd b/vignettes/FAQs.Rmd
index 259aa9b8..740dab74 100644
--- a/vignettes/FAQs.Rmd
+++ b/vignettes/FAQs.Rmd
@@ -113,7 +113,7 @@ You are free to declare any value in the model to be a trajectory or parameter.
 
 ### Can I get a new model specification with fitted parameters as defaults, from a calibrated calibrator object?
 
-Coming soon!
+[Yes](https://canmod.github.io/macpan2/reference/mp_optimized_spec).  [Here](https://canmod.github.io/macpan2/reference/mp_tmb_calibrator.html#ref-examples) is a simple example.
 
 ## Plotting Simulated Trajectories
 
diff --git a/vignettes/state_updaters.Rmd b/vignettes/state_updaters.Rmd
index 910af284..39413b03 100644
--- a/vignettes/state_updaters.Rmd
+++ b/vignettes/state_updaters.Rmd
@@ -1,6 +1,8 @@
 ---
 title: "ODE Solvers, Process Error, and Difference Equations"
-output: rmarkdown::html_vignette
+output: 
+  rmarkdown::html_vignette:
+    toc: true
 vignette: >
   %\VignetteIndexEntry{ODE Solvers, Process Error, and Difference Equations}
   %\VignetteEngine{knitr::rmarkdown}
@@ -17,9 +19,11 @@ set.seed(12L)
 
 [![status](https://img.shields.io/badge/status-working%20draft-red)](https://canmod.github.io/macpan2/articles/vignette-status#working-draft)
 
-The `McMasterPandemic` project has focused on using discrete time difference equations to solve the dynamical system. Here we describe experimental features that allow one to update the data differently. In particular, it is now possible to use a Runge-Kutta 4 ODE solver (to approximate the solution to the continuous time analogue) as well as the Euler Multinomial distribution to generate process error. These features have in principle been available for most of the life of `macpan2`, but only recently have they been given a more convenient interface.
+## Introduction
 
-In order to use this interface, models must be rewritten using explicit declarations of state updates. By explicitly declaring what expressions correspond to state updates, `macpan2` is able to modify the method of updating without further user input. This makes it easier to compare difference equations (i.e. Euler steps), ODEs (i.e. Runge-Kutta), and process error (i.e. Euler Multinomial) runs for the same underlying dynamical model.
+The `McMasterPandemic` project has focused on using discrete time difference equations to solve the dynamical system. Here we describe experimental features that allow one to update the data differently. In particular, it is now possible to use a Runge-Kutta 4 Ordinary Differential Equation (ODE) solver (to approximate the solution to the continuous time analogue) as well as the Euler Multinomial distribution to generate process error. These features have in principle been available for most of the life of `macpan2`, but only recently have they been given a more convenient interface.
+
+In order to use this interface, models must be rewritten using explicit declarations of state updates. By explicitly declaring what expressions correspond to state updates, `macpan2` is able to modify the method of updating without further user input. This makes it easier to compare runs using difference equations (i.e. Euler steps), ODEs (i.e. Runge-Kutta), and stochastic difference equations for the same underlying dynamical model.
 
 To declare a state update one replaces formulas in the `during` argument of a model spec to include calls to the `mp_per_capita_flow()` function.  As an example we will modify the SI model to include explicit state updates.  The original form of the model looks like this.
 ```{r}
@@ -40,7 +44,7 @@ The modified version looks like this.
 ```{r}
 si_explicit = mp_tmb_model_spec(
     before = list(S ~ N - I)
-  , during = list(mp_per_capita_flow("S", "I", infection ~ beta * I / N))
+  , during = list(mp_per_capita_flow("S", "I", "beta * I / N", "infection"))
   , default = list(I = 1, N = 100, beta = 0.25)
 )
 print(si_explicit)
@@ -56,8 +60,8 @@ si_rk4 = (si_explicit
   |> mp_rk4()
   |> mp_simulator(time_steps = 50, outputs = "infection")
 )
-si_euler_multinomial = (si_explicit 
-  |> mp_euler_multinomial()
+si_discrete_stoch = (si_explicit 
+  |> mp_discrete_stoch()
   |> mp_simulator(time_steps = 50, outputs = "infection")
 )
 ```
@@ -68,7 +72,7 @@ library(dplyr)
 trajectory_comparison = list(
     euler = mp_trajectory(si_euler)
   , rk4 = mp_trajectory(si_rk4)
-  , euler_multinomial = mp_trajectory(si_euler_multinomial)
+  , discrete_stoch = mp_trajectory(si_discrete_stoch)
 ) |> bind_rows(.id = "update_method")
 print(head(trajectory_comparison))
 ```
@@ -79,18 +83,20 @@ library(ggplot2)
   |> rename(`Time Step` = time, Incidence = value)
   |> ggplot()
   + geom_line(aes(`Time Step`, Incidence, colour = update_method))
+  + theme_bw()
 )
 ```
 
 The incidence trajectory is different for the three update methods, even though the initial values of the state variables were identical in all three cases.
 ```{r}
 mp_initial(si_euler)
-mp_initial(si_euler_multinomial)
+mp_initial(si_discrete_stoch)
 mp_initial(si_rk4)
 ```
 
 The incidence is different, even during the first time step, because each state update method causes different numbers of susceptible individuals to become infectious at each time step. Further, incidence here is defined for a single time step and so this number of new infectious individuals is exactly the incidence.
 
+
 ## Branching Flows and Process Error
 
 ```{r}
@@ -108,24 +114,75 @@ print(siv)
 
 ```{r}
 (siv 
- |> mp_euler_multinomial()
+ |> mp_discrete_stoch()
  |> mp_simulator(50, "infection")
  |> mp_trajectory()
  |> ggplot()
   + geom_line(aes(time, value))
+  + theme_bw()
+)
+```
+
+
+## Time Steps
+
+```{r sir_prevalence, fig.height = 4, fig.width = 7}
+library(tidyr); library(dplyr)
+delta_t = c(0.1, 0.5, 1)
+time = 50
+sir = mp_tmb_library("starter_models", "sir", package = "macpan2")
+sims = list()
+run = function(model) {
+  (model
+    |> mp_simulator(time/dt, "I")
+    |> mp_trajectory()
+    |> mutate(time = time * dt)
+  )
+}
+for (dt in delta_t) {
+  sims[[sprintf("euler_%s" , dt)]] = sir |> mp_euler(dt)  |> run()
+  sims[[sprintf("hazard_%s", dt)]] = sir |> mp_hazard(dt) |> run()
+  sims[[sprintf("rk4_%s"   , dt)]] = sir |> mp_rk4(dt)    |> run()
+}
+(sims
+  |> bind_rows(.id = "dt")
+  |> separate(dt, c("method", "dt"), "_")
+  |> mutate(method = factor(method, levels = c("rk4", "euler", "hazard")))
+  |> mutate(dt = factor(dt, levels = c("1", "0.5", "0.1")))
+  |> rename(prevalence = value)
+  |> ggplot()
+  + geom_line(aes(time, prevalence, colour = dt, linetype = method))
+  + theme_bw()
 )
 ```
 
 
+$$
+t_n = \delta_t n
+$$
+
+$$
+t_{n + 1} =  t_n \delta_t
+$$
+
+$$
+t_{n + 1} =  \delta_t (n + 1)
+$$
+
+$$
+t_n + \frac{\delta_t}{2} = \delta_t n + \frac{\delta_t}{2} = \delta_t \left(n + \frac{1}{2}\right)
+$$
+
 ## Internal Design
 
 ```{r}
-sir = mp_tmb_model_spec(
+spec = mp_tmb_model_spec(
   during = list(
       N ~ S + I + R
     , mp_per_capita_flow("S", "I", "beta * I / N", "infection")
     , mp_per_capita_flow("I", "R", "(1 - p) * gamma", "recovery")
     , mp_per_capita_flow("I", "D", "p * gamma", "death")
+    , mp_inflow("I", "mu", "importation")
   )
 )
 ```
@@ -141,7 +198,7 @@ The simplest way to define these lists is using two-sided R formulas. But in the
 
 The list of change components, which is just the `during` field but ensuring that all elements are valid `ChangeComponent` objects, can be obtained using the `change_components()` method in a model spec object.
 ```{r}
-cc = sir$change_components()
+cc = spec$change_components()
 print(cc)
 ```
 
@@ -220,12 +277,12 @@ The `update_state()` method makes use of this `flow_frame()` to produce the
 
 
 ```{r}
-si = mp_tmb_library("starter_models", "sir", package = "macpan2")
-si$change_model
-si$change_model$flow_frame()
-si$change_model$change_frame()
-si$change_model$update_state()
-si$change_model$update_flows()
+sir = mp_tmb_library("starter_models", "sir", package = "macpan2")
+sir$change_model
+sir$change_model$flow_frame()
+sir$change_model$change_frame()
+sir$change_model$update_state()
+sir$change_model$update_flows()
 ```
 ```{r}
 spec = mp_tmb_model_spec(
@@ -245,14 +302,34 @@ The constructors of these `UpdateMethod` classes often have a field of class `Ch
 
 
 
-## Relationship to Ordinary Differential Equation Solvers
+## Mathematical Description
+
+We distinguish four kinds of flows that differ in how they are parameterized and how they affect state variables.
+
+### Per-capita Flow Between Two Compartments
+
+This is the most common kind of flow. It is characterized by two compartments and the per-capita rate of flow from one of them to the other. We use the symbol $r$ to indicate per-capita flow, and subscripts to indicate the source and sink compartments. For example, in the standard SIR model there are two flows and both are per-capita flows. The infection flow is from S to I at per-capita rate $r_{SI} = \frac{\beta I}{N}$ and the recovery flow is from `I` to `R` at per-capita rate $r_{IR} = gamma$.
+
+The per-capita flow between two compartemnts is given by 
+
+How S, I, and R are updated in this SIR model depends on the state update method, but the simplest approach is to use differential equations.
+
+$$
+\begin{align}
+\frac{dS}{dt} = \
+\end{align}
+$$
+
 
 It is instructive to view these state update methods as approximate solutions to ordinary differential equations (ODEs). We consider ODEs of the following form.
 
 $$
-\frac{dx_i}{dt} = \underbrace{\sum_j x_j r_{ji}}_{\text{inflow}} - \underbrace{\sum_j x_i r_{ij}}_{\text{outflow}}
+\frac{dx_i}{dt} = 
+  \underbrace{\omega_i}_{\text{exogenous change}} +
+  \underbrace{\sum_j (x_j r_{ji} + \rho_{ji})}_{\text{inflow}} - 
+  \underbrace{\sum_j (x_i r_{ij} + \rho_{ij})}_{\text{outflow}}
 $$
-Where the per-capita rate of flowing from compartment $i$ to compartment $j$ is $r_{ij}$, and $x_i$ is the number of individuals in the $i$th compartment. We allow each $r_{ij}$ to depend on any number of state variables and time-varying parameters. For example, for an SIR model we have  $x_S, x_I, x_R$ (for readability we use $S, I, R$). In this case $r_{SI} = \beta I / N$, $r_{IR} = \gamma$, and all other values of $r_{ij}$ are zero. Here the force of infection, $r_{SI}$, depends on a state variable $I$.
+Where the per-capita rate of flowing from compartment $i$ to compartment $j$ is $r_{ij}$, the associated the rate of exogeneous change of compartment $i$ is $\rho_i$, and $x_i$ is the number of individuals in the $i$th compartment. We allow each $r_{ij}$ to depend on any number of state variables and time-varying parameters. For example, for an SIR model we have  $x_S, x_I, x_R$ (for readability we use $S, I, R$). In this case $r_{SI} = \beta I / N$, $r_{IR} = \gamma$, and all other values of $r_{ij}$ are zero. Here the force of infection, $r_{SI}$, depends on a state variable $I$.
 
 Although each state-update method presented below can be thought of as an approximate solution to this ODE, they can also be thought of as dynamical models in their own right. For example, the Euler-multinomial model is a useful model of process error.