diff --git a/Cargo.lock b/Cargo.lock
index 8787dd928..8816ec670 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -1435,9 +1435,9 @@ dependencies = [
 
 [[package]]
 name = "pharmsol"
-version = "0.29.0"
+version = "0.29.2"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "d956a0a3ae74b4a02c16cd3dff3153a7b958f3f8c48b397f469df6151c1fc339"
+checksum = "dfae22a3e9cfc7642c4bc306827090e5341c7cff48be57e47d890c5de9700c8b"
 dependencies = [
  "ahash",
  "argmin",
@@ -1466,9 +1466,9 @@ dependencies = [
 
 [[package]]
 name = "pharmsol-dsl"
-version = "0.29.0"
+version = "0.29.2"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "80ff36016934559d648d3d3dcdf99003c2f314ff7a00e97d8bc471c609e0db0f"
+checksum = "fab86f195b2ded3a73c1e52f9d96ab150162637cc4a54cf95a616b7f856a192a"
 dependencies = [
  "serde",
  "serde_json",
@@ -1477,9 +1477,9 @@ dependencies = [
 
 [[package]]
 name = "pharmsol-macros"
-version = "0.29.0"
+version = "0.29.2"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "44344f94af37484c5f09f75af566b23baf9643f06076aa27ef741f5e657b47ba"
+checksum = "16bded859c1bfaf2ecd34b05cba1074d89c931e86565546c12eb81aea05e9a77"
 dependencies = [
  "pharmsol-dsl",
  "proc-macro-crate",
@@ -1508,9 +1508,9 @@ dependencies = [
 
 [[package]]
 name = "pmcore"
-version = "0.28.0"
+version = "0.29.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "732c974a86e0970239056ed726e0a7269856a957d3095c024461813cd65bd4da"
+checksum = "941b66ea8169eb15284e24ac63b51da18b81c5abb501d396ce28b0750a49b390"
 dependencies = [
  "anyhow",
  "argmin",
diff --git a/R/PM_model.R b/R/PM_model.R
index 4a0ecac25..1faa962d8 100644
--- a/R/PM_model.R
+++ b/R/PM_model.R
@@ -252,17 +252,20 @@ PM_model <- R6::R6Class(
     #'   - `dx[i]` for the change in amount with respect to time (i.e., \eqn{dx/dt}),
     #' where `i` is the compartment number,
     #'   - `x[i]` for the compartment amount, where `i` is the compartment number.
-    #'   - `rateiv[j]` for the infusion rate of input `j`, where `j` is the input number
-    #' in the data corresponding to doses for that input.
-    #'   - Bolus doses are indicated by `DUR = 0` for dose events in the
-    #' data. Currently only one bolus input is allowed, which goes into compartment 1
-    #' and is not modifiable. It does not appear in the differential equations.
+    #'   - `rateiv[j]` or `R[j]` for the infusion rate of input `j`, where `j` is
+    #' the input number in the data. An infusion may be a standalone additive term
+    #' or one exact product with a scale, such as `R[1] * iv_scale`. The scale is
+    #' evaluated continuously in the differential equation.
+    #'   - `B[j]` or `bolus[j]` for a bolus input. A bolus may be a standalone
+    #' additive term or one exact product with a scale, such as `B[1] * v`.
+    #' The scale is applied to the bolus dose and does not affect an infusion with
+    #' the same input number.
     #'
     #'     For example,
     #'     ```
     #'     eqn = function() {
-    #'      dx[1] = -ka * x[1]
-    #'      dx[2] = rateiv[1] + ka * x[1] - ke * x[2]
+    #'      dx[1] = B[1] * fa1 - ka * x[1]
+    #'      dx[2] = rateiv[1] * iv_scale + ka * x[1] - ke * x[2]
     #'     }
     #'     ```
     #' * **Additional equations** in R code can be defined in this block, which are similar to
@@ -286,7 +289,8 @@ PM_model <- R6::R6Class(
     #' As for `eqn`, additional equations in R code can be defined in this block,
     #' but will only be available within the `lag` block.
     #' @param fa A function defining the bioavailability (fraction absorbed) equations,
-    #' similar to `lag`.
+    #' similar to `lag`. The `lag` and `fa` blocks also apply to analytical bolus
+    #' models; they are not valid for analytical infusion-only models.
     #'
     #' Example:
     #' ```
@@ -308,7 +312,8 @@ PM_model <- R6::R6Class(
     #' This sets the initial amount of drug in compartment 2 to the value
     #' of a covariate `init2` multiplied by the volume of the compartment,
     #' `V`, assuming `V` is either a primary parameter or defined in the
-    #' `sec` block.
+    #' `sec` block. Initial conditions are evaluated before the first event of
+    #' every occasion, including occasions started by an `EVID = 4` row.
     #'
     #' As for `eqn`, additional equations in R code can be defined in this block,
     #' but will only be available within the `ini` block.
@@ -471,7 +476,10 @@ PM_model <- R6::R6Class(
             }
           })
           self$arg_list$x <- NULL
-          self$dsl <- x$dsl
+          # `dsl` is rendered from `arg_list`, so it is never inherited. A
+          # PM_model saved before a DSL change would otherwise replay stale
+          # source instead of rendering the model it defines.
+          self$dsl <- NULL
         } else {
           cli::cli_abort(c(
             "x" = "Non supported input for {.arg x}: {typeof(x)}",
diff --git a/R/model_dsl.R b/R/model_dsl.R
index 13dd509ee..af14de010 100644
--- a/R/model_dsl.R
+++ b/R/model_dsl.R
@@ -17,8 +17,8 @@
 # Naming conventions used by the emitter:
 #   * States   `x[i]`            -> `x{i}`          (declared in `states = ...`)
 #   * Outputs  `Y[i]` / `y[i]`   -> `outeq_{i}`     (declared in `outputs = ...`)
-#   * Inputs   `b[j]`/`bolus[j]` -> `bolus(input_{j}) -> x{k}`    (route)
-#              `rateiv[j]`/`r[j]`-> `infusion(input_{j}) -> x{k}` (route)
+#   * ODE inputs `b[j]`/`bolus[j]` -> `bolus(input_{j})`       (RHS term)
+#                `rateiv[j]`/`r[j]`-> `infusion(input_{j})`    (RHS term)
 #   * Params / covariates keep their (lower-cased) names.
 #
 # pharmsol identifies routes and outputs by label and rejects bare numeric
@@ -27,10 +27,9 @@
 # renumbered: the data's `INPUT`/`OUTEQ` values are canonicalised to the same
 # labels before the backend reads them.
 #
-# Infusions and boluses are declared as routes into the compartment in which they
-# appear in the derivative equations; the corresponding `rateiv[j]`/`b[j]` terms
-# are stripped from the derivative because the DSL runtime injects them
-# automatically from the route declaration.
+# ODE inputs and scales remain on the derivative RHS. The DSL infers their routes
+# and lowers them into the existing simulator machinery. Analytical models keep
+# explicit route declarations and route properties.
 # ---------------------------------------------------------------------------
 
 # The set of DSL intrinsic functions and how R names map to them. R names that
@@ -77,10 +76,16 @@ expr_to_dsl <- function(expr, allow_if = TRUE) {
       ))
     }
     idx <- as.integer(idx_raw)
-    if (var %in% c("b", "bolus", "rateiv", "r")) {
+    if (var %in% c("b", "bolus")) {
       cli::cli_abort(c(
-        "x" = "Bolus/infusion inputs may only be used as standalone additive terms in derivative equations.",
-        "i" = "Write, for example, {.code dx[1] = -ke * x[1] + rateiv[1]}, not inside a product."
+        "x" = "Bolus inputs may only be standalone additive terms or one exact product with a scale.",
+        "i" = "Use {.code B[1]}, {.code B[1] * scale}, or {.code scale * B[1]}."
+      ))
+    }
+    if (var %in% c("rateiv", "r")) {
+      cli::cli_abort(c(
+        "x" = "Infusion inputs may only be standalone additive terms or one exact product with a scale.",
+        "i" = "Use {.code R[1]}, {.code R[1] * scale}, or {.code scale * R[1]}."
       ))
     }
     return(sprintf("%s%d", var, idx))
@@ -200,6 +205,36 @@ dsl_route_of <- function(expr) {
   NULL
 }
 
+# Recognize one exact linear route product, with the route on either side.
+dsl_scaled_route_of <- function(expr) {
+  if (!is.call(expr) || !identical(expr[[1]], as.name("*")) || length(expr) != 3) {
+    return(NULL)
+  }
+
+  left <- dsl_route_of(expr[[2]])
+  right <- dsl_route_of(expr[[3]])
+  if (is.null(left) == is.null(right)) {
+    return(NULL)
+  }
+
+  route <- if (!is.null(left)) left else right
+  route$scale <- if (!is.null(left)) expr[[3]] else expr[[2]]
+  if (route$kind == "bolus") {
+    uses_state_or_input <- function(value) {
+      if (!is.call(value)) return(FALSE)
+      if (identical(value[[1]], as.name("["))) {
+        name <- tolower(as.character(value[[2]]))
+        if (name %in% c("x", "dx", "b", "bolus", "r", "rateiv")) return(TRUE)
+      }
+      any(vapply(as.list(value)[-1], uses_state_or_input, logical(1)))
+    }
+    if (uses_state_or_input(route$scale)) {
+      cli::cli_abort("A bolus scale cannot depend directly on state or dose inputs.")
+    }
+  }
+  route
+}
+
 # Return the sorted, unique data input indices referenced by bolus (`b[]` /
 # `bolus[]`) or infusion (`rateiv[]` / `r[]`) terms in a model equation function.
 # Used to validate that the model represents every dose input present in the data.
@@ -235,7 +270,7 @@ dsl_join_terms <- function(terms) {
   pieces <- character(0)
   for (i in seq_along(terms)) {
     t <- terms[[i]]
-    es <- expr_to_dsl(t$expr, allow_if = allow_if)
+    es <- if (!is.null(t$dsl)) t$dsl else expr_to_dsl(t$expr, allow_if = allow_if)
     if (i == 1) {
       pieces <- if (t$sign < 0) sprintf("-(%s)", es) else es
     } else {
@@ -245,7 +280,7 @@ dsl_join_terms <- function(terms) {
   pieces
 }
 
-# Convert a derivative RHS into (routes, stripped DSL expression) for the given
+# Convert a derivative RHS into (routes, explicit-input DSL expression) for the given
 # destination compartment index `comp`.
 dsl_extract_routes <- function(rhs, comp) {
   terms <- dsl_flatten_add(rhs)
@@ -253,6 +288,9 @@ dsl_extract_routes <- function(rhs, comp) {
   kept <- list()
   for (t in terms) {
     route <- dsl_route_of(t$expr)
+    if (is.null(route)) {
+      route <- dsl_scaled_route_of(t$expr)
+    }
     if (!is.null(route)) {
       if (t$sign < 0) {
         cli::cli_abort(c(
@@ -262,6 +300,11 @@ dsl_extract_routes <- function(rhs, comp) {
       }
       route$comp <- comp
       routes[[length(routes) + 1]] <- route
+      input <- sprintf("%s(%s)", route$kind, pm_input_label(route$input))
+      if (!is.null(route$scale)) {
+        input <- sprintf("%s * (%s)", input, expr_to_dsl(route$scale, allow_if = FALSE))
+      }
+      kept[[length(kept) + 1]] <- list(sign = t$sign, dsl = input)
     } else {
       kept[[length(kept) + 1]] <- t
     }
@@ -356,9 +399,11 @@ dsl_sec_block <- function(fun) {
 # `target` is the DSL property name ("lag" or "fa"); the R block assigns to
 # `lag[j]` / `fa[j]` where `j` is the 1-based input index.
 dsl_route_property_block <- function(fun, target) {
-  exprs <- dsl_body_stmts(fun)
+  exprs <- if (is.null(fun)) list() else dsl_body_stmts(fun)
   derived <- character(0)
   lines <- character(0)
+  values <- list()
+  seen_inputs <- character(0)
   for (e in exprs) {
     if (!dsl_is_assign(e)) {
       cli::cli_abort("Only assignments are supported in the {target} block for the DSL backend.")
@@ -371,12 +416,19 @@ dsl_route_property_block <- function(fun, target) {
       if (tgt != target) {
         cli::cli_abort("Unexpected indexed assignment to {.code {tgt}[{idx}]} in {target} block.")
       }
-      lines <- c(lines, sprintf("%s(%s) = %s", target, pm_input_label(idx), expr_to_dsl(rhs)))
+      key <- as.character(idx)
+      if (key %in% seen_inputs) {
+        cli::cli_abort("The {target} block assigns input {idx} more than once.")
+      }
+      seen_inputs <- c(seen_inputs, key)
+      values[[key]] <- expr_to_dsl(rhs)
+      lines <- c(lines, sprintf("%s(%s) = %s", target, pm_input_label(idx), values[[key]]))
     } else {
       derived <- c(derived, sprintf("%s = %s", tolower(as.character(lhs)), expr_to_dsl(rhs)))
     }
   }
-  list(derived = derived, lines = lines)
+
+  list(derived = derived, lines = lines, values = values)
 }
 
 # Emit `init(...)` statements from an initial-conditions block.
@@ -406,12 +458,12 @@ dsl_ini_block <- function(fun) {
 
 # Finalize route usages into concrete DSL routes.
 #
-# Each usage is `list(kind, input, comp)` where `input` is the identifier used in
-# the model equations (`b[j]` / `rateiv[j]`). Route labels are unique *per kind*
+# Each usage is `list(kind, input, comp, scale)` where `scale` is present for an
+# exact route-input product. Route labels are unique *per kind*
 # in the pharmsol DSL, so the same label may declare both a bolus and an
 # infusion; no relabelling of the data is ever required.
 #
-# Returns a list of `{kind, label, comp}` route declarations.
+# Returns a list of `{kind, input, label, comp, scale}` route declarations.
 dsl_finalize_routes <- function(routes) {
   by_input <- list()
   seen_inputs <- integer(0)
@@ -430,10 +482,8 @@ dsl_finalize_routes <- function(routes) {
     grp <- by_input[[as.character(inp)]]
 
     for (kd in c("bolus", "infusion")) {
-      comps <- unique(vapply(
-        Filter(function(r) identical(r$kind, kd), grp),
-        function(r) as.integer(r$comp), integer(1)
-      ))
+      kind_routes <- Filter(function(r) identical(r$kind, kd), grp)
+      comps <- unique(vapply(kind_routes, function(r) as.integer(r$comp), integer(1)))
       if (length(comps) == 0) next
       if (length(comps) > 1) {
         cli::cli_abort(c(
@@ -442,8 +492,17 @@ dsl_finalize_routes <- function(routes) {
         ))
       }
 
+      scaled <- Filter(function(r) !is.null(r$scale), kind_routes)
+      if ((kd == "bolus" || length(scaled) > 0) && length(kind_routes) > 1) {
+        cli::cli_abort(c(
+          "x" = "The {kd} input {inp} is used more than once.",
+          "i" = "Use one input term with a combined scale."
+        ))
+      }
+      scale <- if (length(scaled) == 1) scaled[[1]]$scale else NULL
+
       final_routes[[length(final_routes) + 1L]] <- list(
-        kind = kd, label = pm_input_label(inp), comp = comps[[1]]
+        kind = kd, input = inp, label = pm_input_label(inp), comp = comps[[1]], scale = scale
       )
     }
   }
@@ -584,6 +643,7 @@ model_to_dsl <- function(model) {
   # ---- ODE model ----
   eqn <- dsl_eqn_block(arg_list$eqn)
   derived <- c(derived, eqn$derived)
+  routes <- dsl_finalize_routes(eqn$routes)
 
   out <- dsl_out_block(arg_list$out)
   derived <- c(derived, out$derived)
@@ -602,11 +662,19 @@ model_to_dsl <- function(model) {
     lag_lines <- lag$lines
   }
 
-  fa_lines <- character(0)
+  # Keep the established R fa block, but express its effect on the ODE RHS.
   if (!is.null(arg_list$fa)) {
     fa <- dsl_route_property_block(arg_list$fa, "fa")
     derived <- c(derived, fa$derived)
-    fa_lines <- fa$lines
+    bolus_inputs <- vapply(Filter(function(r) r$kind == "bolus", routes),
+      function(r) as.character(r$input), character(1))
+    for (key in names(fa$values)) {
+      if (!key %in% bolus_inputs) {
+        cli::cli_abort("The fa block references input {key}, which has no bolus term.")
+      }
+      input <- sprintf("bolus(%s)", pm_input_label(key))
+      eqn$dx <- gsub(input, sprintf("%s * (%s)", input, fa$values[[key]]), eqn$dx, fixed = TRUE)
+    }
   }
 
   # Number of states and outputs.
@@ -621,26 +689,17 @@ model_to_dsl <- function(model) {
   states <- paste0("x", seq_len(n_states))
   outputs <- pm_output_label(seq_len(n_out))
 
-  routes <- dsl_finalize_routes(eqn$routes)
-  route_lines <- vapply(routes, function(r) {
-    sprintf("%s(%s) -> x%d", r$kind, r$label, r$comp)
-  }, character(1))
-
-  # Assemble the DSL text in an order that respects definite assignment:
-  # declarations, routes, derived values, route properties, initial conditions,
-  # derivatives, and finally outputs.
+  # Inputs remain explicit in the derivatives; the DSL infers their routes.
+  # Derived values precede lag, initial conditions, and derivatives.
   lines <- c(
     header,
     sprintf("states = %s", paste(states, collapse = ", ")),
     sprintf("outputs = %s", paste(outputs, collapse = ", ")),
     "",
-    route_lines,
-    if (length(route_lines) > 0) "" else NULL,
     derived,
     if (length(derived) > 0) "" else NULL,
     lag_lines,
-    fa_lines,
-    if (length(lag_lines) > 0 || length(fa_lines) > 0) "" else NULL,
+    if (length(lag_lines) > 0) "" else NULL,
     init_lines,
     if (length(init_lines) > 0) "" else NULL,
     eqn$dx,
@@ -670,6 +729,24 @@ dsl_analytical <- function(model, header, derived, parameters) {
   out <- dsl_out_block(arg_list$out)
   derived <- c(derived, out$derived)
 
+  has_absorption <- stringr::str_detect(structure, "absorption")
+  if (!has_absorption && (!is.null(arg_list$lag) || !is.null(arg_list$fa))) {
+    cli::cli_abort("The `lag` and `fa` blocks can only be used with analytical bolus models.")
+  }
+
+  lag_lines <- character(0)
+  fa_lines <- character(0)
+  if (has_absorption && !is.null(arg_list$lag)) {
+    lag <- dsl_route_property_block(arg_list$lag, "lag")
+    derived <- c(derived, lag$derived)
+    lag_lines <- lag$lines
+  }
+  if (has_absorption && !is.null(arg_list$fa)) {
+    fa <- dsl_route_property_block(arg_list$fa, "fa")
+    derived <- c(derived, fa$derived)
+    fa_lines <- fa$lines
+  }
+
   # The DSL analytical structures require specific derived-parameter names (e.g.
   # `kcp`, `kpc`, `vc`). The Pmetrics model-library templates use their own
   # parameter names, so emit derived aliases mapping the library names to the
@@ -689,7 +766,7 @@ dsl_analytical <- function(model, header, derived, parameters) {
   # Declare the dose route. Absorption ("bolus") templates receive a bolus into
   # the depot (x1); IV templates receive an infusion into the central
   # compartment (x1).
-  route_line <- if (stringr::str_detect(structure, "absorption")) {
+  route_line <- if (has_absorption) {
     sprintf("bolus(%s) -> x1", pm_input_label(1))
   } else {
     sprintf("infusion(%s) -> x1", pm_input_label(1))
@@ -705,6 +782,9 @@ dsl_analytical <- function(model, header, derived, parameters) {
     "",
     derived,
     if (length(derived) > 0) "" else NULL,
+    lag_lines,
+    fa_lines,
+    if (length(lag_lines) > 0 || length(fa_lines) > 0) "" else NULL,
     out$out
   )
 
diff --git a/data/NPex.rda b/data/NPex.rda
index be0103957..d3a899e08 100755
Binary files a/data/NPex.rda and b/data/NPex.rda differ
diff --git a/data/badData.rda b/data/badData.rda
index 39d67dbca..79b6f7ff1 100755
Binary files a/data/badData.rda and b/data/badData.rda differ
diff --git a/data/dataEx.rda b/data/dataEx.rda
index 0fe3decaa..6da32b7de 100755
Binary files a/data/dataEx.rda and b/data/dataEx.rda differ
diff --git a/data/modEx.rda b/data/modEx.rda
index 98583c3be..f6960128d 100755
Binary files a/data/modEx.rda and b/data/modEx.rda differ
diff --git a/data/simEx.rda b/data/simEx.rda
index de2f20c36..00866a446 100755
Binary files a/data/simEx.rda and b/data/simEx.rda differ
diff --git a/man/PM_model.Rd b/man/PM_model.Rd
index b944b7ba0..3055f398c 100755
--- a/man/PM_model.Rd
+++ b/man/PM_model.Rd
@@ -301,17 +301,20 @@ Use the following notation in equations:
 \item \code{dx[i]} for the change in amount with respect to time (i.e., \eqn{dx/dt}),
 where \code{i} is the compartment number,
 \item \code{x[i]} for the compartment amount, where \code{i} is the compartment number.
-\item \code{rateiv[j]} for the infusion rate of input \code{j}, where \code{j} is the input number
-in the data corresponding to doses for that input.
-\item Bolus doses are indicated by \code{DUR = 0} for dose events in the
-data. Currently only one bolus input is allowed, which goes into compartment 1
-and is not modifiable. It does not appear in the differential equations.
+\item \code{rateiv[j]} or \code{R[j]} for the infusion rate of input \code{j}, where \code{j} is
+the input number in the data. An infusion may be a standalone additive term
+or one exact product with a scale, such as \code{R[1] * iv_scale}. The scale is
+evaluated continuously in the differential equation.
+\item \code{B[j]} or \code{bolus[j]} for a bolus input. A bolus may be a standalone
+additive term or one exact product with a scale, such as \code{B[1] * v}.
+The scale is applied to the bolus dose and does not affect an infusion with
+the same input number.
 
 For example,
 
 \if{html}{\out{<div class="sourceCode">}}\preformatted{eqn = function() \{
- dx[1] = -ka * x[1]
- dx[2] = rateiv[1] + ka * x[1] - ke * x[2]
+ dx[1] = B[1] * fa1 - ka * x[1]
+ dx[2] = rateiv[1] * iv_scale + ka * x[1] - ke * x[2]
 \}
 }\if{html}{\out{</div>}}
 }
@@ -338,7 +341,8 @@ this code could be used to model delayed absorption if an antacid is present.
 As for \code{eqn}, additional equations in R code can be defined in this block,
 but will only be available within the \code{lag} block.}
       \item{\code{fa}}{A function defining the bioavailability (fraction absorbed) equations,
-similar to \code{lag}.
+similar to \code{lag}. The \code{lag} and \code{fa} blocks also apply to analytical bolus
+models; they are not valid for analytical infusion-only models.
 
 Example:
 
@@ -362,7 +366,8 @@ Example:
 This sets the initial amount of drug in compartment 2 to the value
 of a covariate \code{init2} multiplied by the volume of the compartment,
 \code{V}, assuming \code{V} is either a primary parameter or defined in the
-\code{sec} block.
+\code{sec} block. Initial conditions are evaluated before the first event of
+every occasion, including occasions started by an \code{EVID = 4} row.
 
 As for \code{eqn}, additional equations in R code can be defined in this block,
 but will only be available within the \code{ini} block.}
diff --git a/src/rust/Cargo.toml b/src/rust/Cargo.toml
index 43df0b1bb..fa592aab0 100755
--- a/src/rust/Cargo.toml
+++ b/src/rust/Cargo.toml
@@ -9,7 +9,7 @@ name = 'pm_rs'
 
 [dependencies]
 extendr-api = "=0.9.0"
-pmcore = { version = "=0.28.0", features = ["dsl"] }
+pmcore = { version = "=0.29.0", features = ["dsl"] }
 rayon = "1.12.0"
 anyhow = "1.0.103"
 tracing = "0.1.44"
diff --git a/tests/testthat/test-io-labels.R b/tests/testthat/test-io-labels.R
index 4cfc5b1f2..76e8b5f83 100644
--- a/tests/testthat/test-io-labels.R
+++ b/tests/testthat/test-io-labels.R
@@ -60,8 +60,9 @@ testthat::test_that("An input driving both a bolus and an infusion keeps one lab
     )
     mod$compile(quiet = TRUE)
 
-    testthat::expect_match(mod$dsl, "bolus(input_1) -> x1", fixed = TRUE)
-    testthat::expect_match(mod$dsl, "infusion(input_1) -> x1", fixed = TRUE)
+    testthat::expect_match(mod$dsl, "bolus(input_1)", fixed = TRUE)
+    testthat::expect_match(mod$dsl, "infusion(input_1)", fixed = TRUE)
+    testthat::expect_false(grepl("->", mod$dsl, fixed = TRUE))
     testthat::expect_false(grepl("input_2", mod$dsl, fixed = TRUE))
     # The shared label must be accepted by the pharmsol backend.
     testthat::expect_no_error(model_parameters(mod$dsl))
diff --git a/tests/testthat/test-model-dsl-dose-scaling.R b/tests/testthat/test-model-dsl-dose-scaling.R
new file mode 100644
index 000000000..9b647751f
--- /dev/null
+++ b/tests/testthat/test-model-dsl-dose-scaling.R
@@ -0,0 +1,155 @@
+testthat::skip_on_cran()
+
+model_to_dsl <- getFromNamespace("model_to_dsl", "Pmetrics")
+
+scaled_bolus_model <- function(eqn_block, fa_block = NULL) {
+  PM_model$new(
+    list(
+      pri = list(
+        kes = ab(0, 5),
+        vs = ab(0, 5),
+        ff = ab(0, 1)
+      ),
+      cov = list(wt = interp()),
+      sec = function() {
+        v <- vs * (wt / 70)
+        ke <- kes * (wt / 70)^-0.25
+      },
+      fa = fa_block,
+      eqn = eqn_block,
+      out = function() {
+        y[1] <- x[1] / v
+        y[2] <- x[2]
+      },
+      err = list(
+        additive(1, c(1, 0, 0, 0), outeq = 1),
+        additive(1, c(1, 0, 0, 0), outeq = 2)
+      )
+    ),
+    compile = FALSE
+  )
+}
+
+testthat::test_that("ODE dose inputs and derived scales stay on the RHS", {
+  mod <- scaled_bolus_model(function() {
+    dx[1] <- B[1] * v + R[1] - ke * x[1]
+    dx[2] <- B[2] + x[1] / v * ff
+  })
+
+  dsl <- model_to_dsl(mod)
+  testthat::expect_match(dsl, "bolus(input_1) * (v)", fixed = TRUE)
+  testthat::expect_match(dsl, "infusion(input_1)", fixed = TRUE)
+  testthat::expect_false(grepl("->", dsl, fixed = TRUE))
+  testthat::expect_false(grepl("fa(input_1)", dsl, fixed = TRUE))
+  testthat::expect_match(dsl, "v = vs * (wt / 70.0)", fixed = TRUE)
+})
+
+testthat::test_that("The bolus may be on either side of the scale", {
+  mod <- scaled_bolus_model(function() {
+    dx[1] <- v * B[1] + R[1] - ke * x[1]
+    dx[2] <- B[2] + x[1] / v * ff
+  })
+
+  testthat::expect_match(model_to_dsl(mod), "bolus(input_1) * (v)", fixed = TRUE)
+})
+
+testthat::test_that("Equation scaling and an explicit fa block are multiplied", {
+  mod <- scaled_bolus_model(
+    eqn_block = function() {
+      dx[1] <- B[1] * v + R[1] - ke * x[1]
+      dx[2] <- B[2] + x[1] / v * ff
+    },
+    fa_block = function() {
+      fa[1] <- ff
+    }
+  )
+
+  dsl <- model_to_dsl(mod)
+  testthat::expect_match(dsl, "bolus(input_1) * (ff) * (v)", fixed = TRUE)
+  testthat::expect_false(grepl("fa(input_1)", dsl, fixed = TRUE))
+  matches <- gregexpr("bolus(input_1)", dsl, fixed = TRUE)[[1]]
+  testthat::expect_equal(sum(matches > 0), 1)
+})
+
+testthat::test_that("A linear infusion scale stays explicit in the derivative", {
+  right <- scaled_bolus_model(function() {
+    dx[1] <- R[1] * v - ke * x[1]
+    dx[2] <- B[2] + x[1] / v * ff
+  })
+  left <- scaled_bolus_model(function() {
+    dx[1] <- v * R[1] - ke * x[1]
+    dx[2] <- B[2] + x[1] / v * ff
+  })
+
+  for (mod in list(right, left)) {
+    dsl <- model_to_dsl(mod)
+    testthat::expect_false(grepl("->", dsl, fixed = TRUE))
+    testthat::expect_match(dsl, "infusion(input_1) * (v)", fixed = TRUE)
+    testthat::expect_false(grepl("fa(input_1)", dsl, fixed = TRUE))
+  }
+})
+
+testthat::test_that("Complex route expressions stay unsupported", {
+  nested_infusion <- scaled_bolus_model(function() {
+    dx[1] <- R[1] * v * ff - ke * x[1]
+    dx[2] <- B[2] + x[1] / v * ff
+  })
+  testthat::expect_error(model_to_dsl(nested_infusion), "Infusion inputs may only")
+
+  nested_bolus <- scaled_bolus_model(function() {
+    dx[1] <- B[1] * v * ff + R[1] - ke * x[1]
+    dx[2] <- B[2] + x[1] / v * ff
+  })
+  testthat::expect_error(model_to_dsl(nested_bolus), "Bolus inputs may only")
+
+  state_dependent <- scaled_bolus_model(function() {
+    dx[1] <- B[1] * x[1] + R[1] - ke * x[1]
+    dx[2] <- B[2] + x[1] / v * ff
+  })
+  testthat::expect_error(model_to_dsl(state_dependent), "cannot depend directly on state")
+})
+
+testthat::test_that("Analytical bolus models emit lag and fa", {
+  mod <- PM_model$new(
+    pri = list(Ka = ab(0.5, 6), Ke = ab(0.1, 1.5), V = ab(25, 120)),
+    eqn = function() {
+      one_comp_bolus
+    },
+    lag = function() {
+      lag[1] <- 2
+    },
+    fa = function() {
+      fa[1] <- 0.5
+    },
+    out = function() {
+      y[1] <- x[2] / V
+    },
+    err = list(proportional(5, c(0.1, 0.15, 0, 0))),
+    compile = FALSE
+  )
+
+  dsl <- model_to_dsl(mod)
+  testthat::expect_match(dsl, "bolus(input_1) -> x1", fixed = TRUE)
+  testthat::expect_match(dsl, "lag(input_1) = 2.0", fixed = TRUE)
+  testthat::expect_match(dsl, "fa(input_1) = 0.5", fixed = TRUE)
+  testthat::expect_no_error(model_parameters(dsl))
+})
+
+testthat::test_that("Analytical infusion models reject lag and fa", {
+  mod <- PM_model$new(
+    pri = list(Ke = ab(0.1, 1.5), V = ab(25, 120)),
+    eqn = function() {
+      one_comp_iv
+    },
+    fa = function() {
+      fa[1] <- 0.5
+    },
+    out = function() {
+      y[1] <- x[1] / V
+    },
+    err = list(proportional(5, c(0.1, 0.15, 0, 0))),
+    compile = FALSE
+  )
+
+  testthat::expect_error(model_to_dsl(mod), "analytical bolus models")
+})
diff --git a/tests/testthat/test-pm-model-from-dsl.R b/tests/testthat/test-pm-model-from-dsl.R
index aec9d3428..ae34b5bde 100644
--- a/tests/testthat/test-pm-model-from-dsl.R
+++ b/tests/testthat/test-pm-model-from-dsl.R
@@ -4,8 +4,7 @@ kind = ode
 params = k, v
 states = central
 outputs = cp
-infusion(input_1) -> central
-dx(central) = rate(input_1) - k * central
+dx(central) = infusion(input_1) - k * central
 out(cp) = central / v ~ continuous()
 "
 
@@ -36,9 +35,7 @@ kind = ode
 params = k
 states = central
 outputs = cp
-bolus(oral) -> central
-infusion(iv) -> central
-dx(central) = -k * central
+dx(central) = bolus(oral) + infusion(iv) - k * central
 out(cp) = central ~ continuous()
 "
 
diff --git a/tests/testthat/test-pmcore-pharmsol-migration.R b/tests/testthat/test-pmcore-pharmsol-migration.R
index 407e72331..37db5485d 100644
--- a/tests/testthat/test-pmcore-pharmsol-migration.R
+++ b/tests/testthat/test-pmcore-pharmsol-migration.R
@@ -23,13 +23,14 @@ build_passthrough_ode_model <- function(solver = NULL) {
   )
 }
 
-test_that("ODE generation produces a DSL model with 1-based indices and an infusion route", {
+test_that("ODE generation produces 1-based labels and an explicit infusion term", {
   mod <- build_passthrough_ode_model("TSIT45")
   mod$compile(quiet = TRUE)
   dsl <- mod$dsl
 
   testthat::expect_match(dsl, "kind = ode")
-  testthat::expect_match(dsl, "infusion\\(input_1\\) -> x1")
+  testthat::expect_match(dsl, "infusion(input_1)", fixed = TRUE)
+  testthat::expect_false(grepl("->", dsl, fixed = TRUE))
   testthat::expect_match(dsl, "dx\\(x1\\)")
   testthat::expect_match(dsl, "out\\(outeq_1\\)")
   # Outputs are 1-based to match the Pmetrics data OUTEQ column.