diff --git a/tests/testthat/setup.R b/tests/testthat/setup.R
index 906863f..ed1138f 100644
--- a/tests/testthat/setup.R
+++ b/tests/testthat/setup.R
@@ -1,6 +1,13 @@
+# =============================================================================
+# Shared models compiled once at test session startup
+# =============================================================================
+
+# --- Library models (standard PK models) ---
 mod_1cmt_iv <- new_ode_model("pk_1cmt_iv")
 mod_2cmt_iv <- new_ode_model("pk_2cmt_iv")
 mod_1cmt_oral <- new_ode_model("pk_1cmt_oral")
+
+# --- Custom 1cmt oral with lagtime ---
 mod_1cmt_oral_lagtime <- new_ode_model(
   code = "
     dAdt[0] = -KA * A[0]
@@ -11,7 +18,9 @@ mod_1cmt_oral_lagtime <- new_ode_model(
   dose = list(cmt = 1, bioav = 1),
   parameters = list(CL = 5, V = 50, KA = 0.5, TLAG = 0.83)
 )
-oral_1cmt_allometric <- new_ode_model( # also timevarying and dose-dependence factor
+
+# --- 1cmt oral with allometric scaling, time-varying CL, dose-dependent V ---
+oral_1cmt_allometric <- new_ode_model(
   code = "
       if(t<168.0) {
         CLi = CL * pow(WT/70, 0.75)
@@ -32,3 +41,174 @@ oral_1cmt_allometric <- new_ode_model( # also timevarying and dose-dependence fa
   declare_variables = c("CLi", "Vi"),
   parameters = list(KA = 0.5, CL = 5, V = 50)
 )
+
+# --- From test_iov.R: IOV models ---
+pars_iov <- list(
+  "kappa_CL_1" = 0,
+  "kappa_CL_2" = 0,
+  "kappa_CL_3" = 0,
+  "eta_CL" = 0,
+  "CL" = 5,
+  "V" = 50,
+  "KA" = 1
+)
+pars_iov_no_iov <- list(
+  "CL" = 5,
+  "V" = 50,
+  "KA" = 1
+)
+pk_iov_none <- new_ode_model(
+  code = "
+      dAdt[1] = -KA * A[1]
+      dAdt[2] = +KA * A[1] -(CL/V) * A[2]
+    ",
+  obs = list(cmt = 2, scale = "V"),
+  dose = list(cmt = 1, bioav = 1),
+  parameters = names(pars_iov_no_iov),
+  cpp_show_code = FALSE
+)
+pk_iov <- new_ode_model(
+  code = "
+      CL_iov = CL * exp(kappa_CL + eta_CL);
+      dAdt[1] = -KA * A[1]
+      dAdt[2] = +KA * A[1] -(CL_iov/V) * A[2]
+    ",
+  iov = list(
+    cv = list(CL = 0.2),
+    n_bins = 3
+  ),
+  obs = list(cmt = 2, scale = "V"),
+  dose = list(cmt = 1, bioav = 1),
+  declare_variables = c("kappa_CL", "CL_iov"),
+  parameters = names(pars_iov),
+  cpp_show_code = FALSE
+)
+
+# --- From test_compare_results.R ---
+dose_in_cmt_2 <- new_ode_model(
+  code = "
+      dAdt[1] = -KA * A[1];
+      dAdt[2] = KA*A[1] -(CL/V) * A[2]
+      dAdt[3] = S2*(A[2]-A[3])
+    ",
+  obs = list(cmt = 2, scale = "V"),
+  dose = list(cmt = 2),
+  cpp_show_code = FALSE
+)
+
+# --- From test_cmt_mapping.R: 1cmt oral with cmt_mapping ---
+pk1cmt_oral_cmt_mapping <- new_ode_model(
+  code = "dAdt[1] = -KA*A[1]; dAdt[2] = KA*A[1] - (CL/V)*A[2];",
+  obs = list(cmt = 2, scale = "V"),
+  cmt_mapping = list(oral = 1, infusion = 2, bolus = 2)
+)
+
+# --- From test_multi_obs.R: Multi-observation model ---
+vars_multi_obs <- c("CONC", "METAB", "METAB2", "ACT")
+pk_multi_obs <- new_ode_model(
+  code = "dAdt[1] = -(CL/V)*A[1]; CONC = 1000*A[1]/V; METAB = CONC/2; METAB2 = CONC * t; ACT = 15",
+  obs = list(variable = vars_multi_obs),
+  declare_variables = vars_multi_obs,
+  cpp_show_code = FALSE
+)
+
+# --- From test_mixture_model.R ---
+covs_mixture <- list(WT = new_covariate(70))
+mod_mixture <- new_ode_model(
+  code = "
+      dAdt[0] = -(CL*(WT/70.0)/V)*A[0];
+    ",
+  pk_code = " ",
+  obs = list(cmt = 1, scale = "V"),
+  mixture = list(CL = list(values = c(5, 15), probability = 0.3)),
+  covariates = covs_mixture
+)
+
+# --- Conditional models (skip on CRAN due to compilation time) ---
+if (identical(Sys.getenv("NOT_CRAN"), "true")) {
+
+  # Library models
+  pk_3cmt_iv <- new_ode_model("pk_3cmt_iv")
+  pk_2cmt_oral <- new_ode_model("pk_2cmt_oral")
+  pk_3cmt_oral <- new_ode_model("pk_3cmt_oral")
+  mod_1cmt_iv_mm <- new_ode_model("pk_1cmt_iv_mm")
+
+  # --- From test_advan_with_auc.R: Models with AUC compartments ---
+  parameters_advan_auc <- list(
+    CL = 10, V = 50, KA = 0.5, Q = 5, V2 = 100, Q2 = 3, V3 = 150, F1 = 1
+  )
+  mod_1cmt_auc <- new_ode_model(
+    code = "dAdt[1] = -(CL/V)*A[1]; dAdt[2] = A[1]/V;",
+    parameters = parameters_advan_auc
+  )
+  mod_2cmt_auc <- new_ode_model(
+    code = "
+    dAdt[1] = -(CL/V)*A[1] - (Q/V)*A[1] + (Q/V2)*A[2];
+    dAdt[2] = +(Q/V)*A[1] - (Q/V2)*A[2];
+    dAdt[3] = A[1]/V;
+  ",
+    parameters = parameters_advan_auc
+  )
+  mod_3cmt_auc <- new_ode_model(
+    code = "
+    dAdt[1] = -(CL/V)*A[1] - (Q/V)*A[1] + (Q/V2)*A[2] - (Q2/V)*A[1] + (Q2/V3)*A[3];
+    dAdt[2] =                (Q/V)*A[1]  -(Q/V2)*A[2]                             ;
+    dAdt[3] =                                           (Q2/V)*A[1] - (Q2/V3)*A[3];
+    dAdt[4] = A[1]/V;
+  ",
+    parameters = parameters_advan_auc
+  )
+
+  # --- From test_timevar_cov.R: 2cmt with time-varying covariates ---
+  mod_2cmt_timevar <- new_ode_model(
+    code = "
+      dAdt[1] = -(Q/V)*A[1] + (Q/V2)*A[2] -(CLi/V)*A[1];
+      dAdt[2] = -(Q/V2)*A[2] + (Q/V)*A[1];
+    ",
+    pk_code = "CLi = CL + CRCL",
+    obs = list(cmt = 2, scale = "V"),
+    covariates = list(CRCL = new_covariate(5)),
+    declare_variables = "CLi",
+    cpp = FALSE
+  )
+
+  # --- From test_t_init.R: Model with state initialization ---
+  mod_t_init <- new_ode_model(
+    code = "CLi = CL; Vi = V; dAdt[1] = -(CLi/Vi)*A[1]; CONC = A[1]/Vi",
+    state_init = "A[1] = TDM_INIT * Vi",
+    parameters = list(CL = 7.67, V = 97.7, TDM_INIT = 500),
+    obs = list(cmt = 1, scale = "Vi"),
+    declare_variables = c("CONC", "CLi", "Vi"),
+    cpp_show_code = FALSE
+  )
+
+  # --- From test_reparametrization.R: Carreno 2cmt model ---
+  covs_carreno <- list(
+    CRCL = new_covariate(5),
+    CL_HEMO = new_covariate(0)
+  )
+  model_carreno <- new_ode_model(
+    code = "
+    CLi = SCLInter + SCLSlope * (CRCL*16.667) + CL_HEMO \
+    Vi = V \
+    Qi = K12 * Vi \
+    V2i = Qi / K21 \
+    dAdt[0] = -(CLi/V)*A[0] - K12*A[0] + K21*A[1] \
+    dAdt[1] = K12*A[0] - K21*A[1] \
+    dAdt[2] = A[0]/V
+  ",
+    parameters = list(
+      V = 25.76, SCLSlope = 0.036, K12 = 2.29, K21 = 1.44,
+      SCLInter = 0.18, TDM_INIT = 0
+    ),
+    declare_variables = c("CLi", "Qi", "Vi", "V2i"),
+    covariates = covs_carreno,
+    obs = list(cmt = 1, scale = "V"),
+    reparametrization = list(
+      "CL" = "SCLInter + SCLSlope * (CRCL*16.667) + CL_HEMO",
+      "V" = "V",
+      "Q" = "K12 * V",
+      "V2" = "(K12 * V) / K21"
+    )
+  )
+}
diff --git a/tests/testthat/test_advan_with_auc.R b/tests/testthat/test_advan_with_auc.R
index ced7ef9..0cc9e8a 100644
--- a/tests/testthat/test_advan_with_auc.R
+++ b/tests/testthat/test_advan_with_auc.R
@@ -1,43 +1,14 @@
+# Uses models and parameters defined in setup.R (conditional, NOT_CRAN only):
+# - mod_1cmt_auc, mod_2cmt_auc, mod_3cmt_auc
+# - parameters_advan_auc
+
 if (identical(Sys.getenv("NOT_CRAN"), "true")) {
   dose <- 100
   interval <- 12
   n_days <- 5
   t_inf <- 1.5
-  parameters <- list(
-    CL = 10,
-    V = 50,
-    KA = 0.5,
-    Q = 5,
-    V2 = 100,
-    Q2 = 3,
-    V3 = 150,
-    F1 = 1
-  )
   t_obs <- c(3, 6, 8, 23, 47)
 
-  ## ODE models for testing
-  mod_1cmt <- new_ode_model(
-    code="dAdt[1] = -(CL/V)*A[1]; dAdt[2] = A[1]/V;",
-    parameters = parameters
-  )
-  mod_2cmt <- new_ode_model(
-    code="
-    dAdt[1] = -(CL/V)*A[1] - (Q/V)*A[1] + (Q/V2)*A[2];
-    dAdt[2] = +(Q/V)*A[1] - (Q/V2)*A[2];
-    dAdt[3] = A[1]/V;
-  ",
-  parameters = parameters
-  )
-  mod_3cmt <- new_ode_model(
-    code="
-    dAdt[1] = -(CL/V)*A[1] - (Q/V)*A[1] + (Q/V2)*A[2] - (Q2/V)*A[1] + (Q2/V3)*A[3];
-    dAdt[2] =                (Q/V)*A[1]  -(Q/V2)*A[2]                             ;
-    dAdt[3] =                                           (Q2/V)*A[1] - (Q2/V3)*A[3];
-    dAdt[4] = A[1]/V;
-  ",
-  parameters = parameters
-  )
-
   ## bolus dataset
   reg_bolus <- new_regimen(
     amt = dose,
@@ -47,7 +18,7 @@ if (identical(Sys.getenv("NOT_CRAN"), "true")) {
   )
   data_bolus <- advan_create_data(
     reg_bolus,
-    parameters = parameters,
+    parameters = parameters_advan_auc,
     cmts = 5,
     t_obs = t_obs
   )
@@ -61,7 +32,7 @@ if (identical(Sys.getenv("NOT_CRAN"), "true")) {
   )
   data_infusion <- advan_create_data(
     reg_infusion,
-    parameters = parameters,
+    parameters = parameters_advan_auc,
     cmts = 6,
     t_obs = t_obs
   )
@@ -71,7 +42,7 @@ test_that("One compartment bolus ADVAN runs", {
   skip_on_cran()
   res1_iv_r <- advan("1cmt_iv_bolus", cpp=FALSE)(data_bolus)
   res1_iv_c <- advan("1cmt_iv_bolus", cpp=TRUE)(data_bolus)
-  res1_iv_ode <- sim(ode = mod_1cmt, regimen = reg_bolus, parameters = parameters, t_obs = t_obs)
+  res1_iv_ode <- sim(ode = mod_1cmt_auc, regimen = reg_bolus, parameters = parameters_advan_auc, t_obs = t_obs)
 
   # AUC R
   expect_equal(round(res1_iv_r[res1_iv_r$TIME %in% t_obs,]$AUC, 5), round(res1_iv_ode[res1_iv_ode$comp == 2,]$y, 5))
@@ -84,7 +55,7 @@ test_that("Two compartment bolus ADVAN runs", {
   skip_on_cran()
   res2_iv_r <- advan("2cmt_iv_bolus", cpp=FALSE)(data_bolus)
   res2_iv_c <- advan("2cmt_iv_bolus", cpp=TRUE)(data_bolus)
-  res2_iv_ode <- sim(ode = mod_2cmt, regimen = reg_bolus, parameters = parameters, t_obs = t_obs)
+  res2_iv_ode <- sim(ode = mod_2cmt_auc, regimen = reg_bolus, parameters = parameters_advan_auc, t_obs = t_obs)
   expect_equal(
     round(res2_iv_r[res2_iv_r$TIME %in% t_obs,]$AUC, 5),
     round(res2_iv_c[res2_iv_c$TIME %in% t_obs,]$AUC, 5)
@@ -106,7 +77,7 @@ test_that("Two compartment infusion ADVAN runs", {
   skip_on_cran()
   res2_inf_r <- advan("2cmt_iv_infusion", cpp=FALSE)(data_infusion)
   res2_inf_c <- advan("2cmt_iv_infusion", cpp=TRUE)(data_infusion)
-  res2_inf_ode <- sim(ode = mod_2cmt, regimen = reg_infusion, parameters = parameters, t_obs = t_obs)
+  res2_inf_ode <- sim(ode = mod_2cmt_auc, regimen = reg_infusion, parameters = parameters_advan_auc, t_obs = t_obs)
 
   expect_equal(
     round(res2_inf_r[res2_inf_r$TIME %in% t_obs,]$AUC, 5),
@@ -131,7 +102,7 @@ test_that("Three compartment bolus ADVAN runs", {
   skip_on_cran()
   res3_iv_r <- advan("3cmt_iv_bolus", cpp=FALSE)(data_bolus)
   res3_iv_c <- advan("3cmt_iv_bolus", cpp=TRUE)(data_bolus)
-  res3_iv_ode <- sim(ode = mod_3cmt, regimen = reg_bolus, parameters = parameters, t_obs = t_obs)
+  res3_iv_ode <- sim(ode = mod_3cmt_auc, regimen = reg_bolus, parameters = parameters_advan_auc, t_obs = t_obs)
   expect_equal(
     round(res3_iv_r[res3_iv_r$TIME %in% t_obs,]$AUC, 5),
     round(res3_iv_c[res3_iv_c$TIME %in% t_obs,]$AUC, 5)
@@ -153,7 +124,7 @@ test_that("Three compartment iv ADVAN runs", {
   skip_on_cran()
   res3_iv_r <- advan("3cmt_iv_infusion", cpp=FALSE)(data_infusion)
   res3_iv_c <- advan("3cmt_iv_infusion", cpp=TRUE)(data_infusion)
-  res3_iv_ode <- sim(ode = mod_3cmt, regimen = reg_infusion, parameters = parameters, t_obs = t_obs)
+  res3_iv_ode <- sim(ode = mod_3cmt_auc, regimen = reg_infusion, parameters = parameters_advan_auc, t_obs = t_obs)
   expect_equal(
     round(res3_iv_r[res3_iv_r$TIME %in% t_obs,]$AUC, 5),
     round(res3_iv_c[res3_iv_c$TIME %in% t_obs,]$AUC, 5)
@@ -170,4 +141,3 @@ test_that("Three compartment iv ADVAN runs", {
     round(res3_iv_ode[res3_iv_ode$comp == 4,]$y, 5)
   )
 })
-
diff --git a/tests/testthat/test_calc_ss_analytic.R b/tests/testthat/test_calc_ss_analytic.R
index 4363472..54ca63c 100644
--- a/tests/testthat/test_calc_ss_analytic.R
+++ b/tests/testthat/test_calc_ss_analytic.R
@@ -1,3 +1,7 @@
+# Uses models defined in setup.R:
+# - mod_1cmt_iv, mod_2cmt_iv, mod_1cmt_oral
+# - pk_3cmt_iv, pk_2cmt_oral, pk_3cmt_oral (conditional, NOT_CRAN only)
+
 # shared parameters
 dose <- 100
 interval <- 12
@@ -6,13 +10,6 @@ reg_oral <- new_regimen(amt = dose, interval = interval, n = n_ss, type = "oral"
 reg_bolus <- new_regimen(amt = dose, interval = interval, n = n_ss, type = "bolus")
 reg_inf <- new_regimen(amt = dose, interval = interval, n = n_ss, type = "infusion")
 t_obs <- max(reg_oral$dose_times) + interval
-# Uses models defined in setup.R
-
-if (identical(Sys.getenv("NOT_CRAN"), "true")) {
-  pk_3cmt_iv <- new_ode_model("pk_3cmt_iv")
-}
-
-#delta <- function(x, ref) { abs(x-ref)/ref }
 
 test_that("1 cmt oral", {
   par <- list(CL = 5, V = 100, KA = 1)
@@ -72,7 +69,6 @@ test_that("1-cmt iv infusion", {
 
 test_that("2-cmt oral", {
   skip_on_cran()
-  pk_2cmt_oral <- new_ode_model("pk_2cmt_oral")
   par <- list(CL = 5, V = 100, Q = 3, V2 = 150, KA = 1)
   res_ana <- calc_ss_analytic(f = "2cmt_oral", dose = dose, interval = interval, parameters = par)
   res_ode <- sim(pk_2cmt_oral, parameters = par, regimen = reg_oral, t_obs = t_obs, only_obs = F)$y
@@ -98,7 +94,6 @@ test_that("2-cmt infusion", {
 
 test_that("3-cmt oral", {
   skip_on_cran()
-  pk_3cmt_oral <- new_ode_model("pk_3cmt_oral")
   par <- list(CL = 5, V = 100, Q = 3, V2 = 150, Q2 = 6, V3 = 250, KA = 1)
   res_ana <- calc_ss_analytic(f = "3cmt_oral", dose = dose, interval = interval, parameters = par)
   res_ode <- sim(pk_3cmt_oral, parameters = par, regimen = reg_oral, t_obs = t_obs, only_obs = F)$y
diff --git a/tests/testthat/test_cmt_mapping.R b/tests/testthat/test_cmt_mapping.R
index 48af64e..dfeea73 100644
--- a/tests/testthat/test_cmt_mapping.R
+++ b/tests/testthat/test_cmt_mapping.R
@@ -1,12 +1,9 @@
-pk1cmt_oral_code <- new_ode_model(
-  code = "dAdt[1] = -KA*A[1]; dAdt[2] = KA*A[1] - (CL/V)*A[2];",
-  obs = list(cmt = 2, scale="V"),
-  cmt_mapping = list(oral = 1, infusion = 2, bolus = 2)
-)
+# Uses model defined in setup.R:
+# - pk1cmt_oral_cmt_mapping
 
 test_that("Compartment mapping is added to attributes", {
-  expect_equal(attr(pk1cmt_oral_code, "cmt_mapping")[["oral"]], 1)
-  expect_equal(attr(pk1cmt_oral_code, "cmt_mapping")[["infusion"]], 2)
+  expect_equal(attr(pk1cmt_oral_cmt_mapping, "cmt_mapping")[["oral"]], 1)
+  expect_equal(attr(pk1cmt_oral_cmt_mapping, "cmt_mapping")[["infusion"]], 2)
 })
 
 test_that("Admin route is interpreted and simulated correctly", {
@@ -18,7 +15,7 @@ test_that("Admin route is interpreted and simulated correctly", {
   )
   p <- list(KA = 1, CL = 5, V = 50)
   res <-  sim_ode(
-    ode = pk1cmt_oral_code,
+    ode = pk1cmt_oral_cmt_mapping,
     parameters = p,
     regimen = regimen,
     only_obs = FALSE
diff --git a/tests/testthat/test_compare_results.R b/tests/testthat/test_compare_results.R
index 73652e9..61120b0 100644
--- a/tests/testthat/test_compare_results.R
+++ b/tests/testthat/test_compare_results.R
@@ -1,31 +1,18 @@
-## models: shared between tests and take a while to compile
-#  - oral models
-## Uses model defined in setup.R
+# Uses models defined in setup.R:
+# - mod_1cmt_oral
+# - mod_1cmt_iv
+# - dose_in_cmt_2
+
 pk1cmt_oral_anal = function(t, dose, KA, V, CL) {
   dose*KA/(V*(KA-CL/V))*(exp(-(CL/V) * t)-exp(-KA * t))
 }
+
 pk1cmt_oral_code <- new_ode_model(
   code = "dAdt[1] = -KA*A[1]; dAdt[2] = KA*A[1] - (CL/V)*A[2];",
-  obs=list(cmt = 2, scale="V")
-)
-
-#   - iv models
-## Uses model defined in setup.R
-
-#   - model with dose cmt specified
-dose_in_cmt_2  <- new_ode_model(
-  code = "
-      dAdt[1] = -KA * A[1];
-      dAdt[2] = KA*A[1] -(CL/V) * A[2]
-      dAdt[3] = S2*(A[2]-A[3])
-    ",
-  obs = list(cmt=2, scale="V"),
-  dose = list(cmt = 2),
-  cpp_show_code = FALSE
+  obs = list(cmt = 2, scale = "V")
 )
 
-
-test_that("Library and custom C++ and code matches analytic soln", {
+test_that("Model from library and custom C++ and code matches analytic solution", {
   p <- list(KA = 1, CL = 5, V = 50)
   t_obs <- c(0:72)
   t_obs2 <- t_obs + 0.1234 # also needs to be producing results with non-integer times
@@ -43,7 +30,7 @@ test_that("Library and custom C++ and code matches analytic soln", {
     only_obs=TRUE
   )
 
-  pk1cmt_oral_code <- sim_ode(
+  pk1cmt_oral_code_res <- sim_ode(
     ode = pk1cmt_oral_code,
     parameters = p,
     duplicate_t_obs = TRUE,
@@ -55,7 +42,7 @@ test_that("Library and custom C++ and code matches analytic soln", {
 
   pk1cmt_oral_anal_res <- pk1cmt_oral_anal(t_obs, dose, p$KA, p$V, p$CL)
   expect_equal(round(pk1cmt_oral_lib$y, 3), round(pk1cmt_oral_anal_res, 3))
-  expect_equal(round(pk1cmt_oral_code$y, 3), round(pk1cmt_oral_anal_res, 3))
+  expect_equal(round(pk1cmt_oral_code_res$y, 3), round(pk1cmt_oral_anal_res, 3))
 })
 
 
@@ -87,6 +74,7 @@ test_that("test bug EmCo 20150925", {
 })
 
 test_that("model size is appropriate (bug: JeHi 20151204)", {
+  skip_on_cran()
   pk3cmt <- new_ode_model(
     code = "
       dAdt[1] = -KA*A[1];
diff --git a/tests/testthat/test_get_var_y.R b/tests/testthat/test_get_var_y.R
index 264843f..fb13027 100644
--- a/tests/testthat/test_get_var_y.R
+++ b/tests/testthat/test_get_var_y.R
@@ -1,6 +1,9 @@
+# Uses models defined in setup.R:
+# - mod_1cmt_iv
+# - mod_2cmt_iv
+# - mod_1cmt_iv_mm (conditional, NOT_CRAN only)
 
 ## Set up simulations to test variance:
-## Uses model defined in setup.R
 reg <- new_regimen(
   amt = 100,
   n = 3,
@@ -130,13 +133,12 @@ test_that("Two compartment model", {
 
 test_that("One compartment with MM kinetics", {
   skip_on_cran()
-  mod3 <- new_ode_model("pk_1cmt_iv_mm")
   par3 <- list(VMAX = 5, KM = 5, V = 10)
   omega3 <- c(0.1,
               0.05, 0.1,
               0.01, 0.01, 0.1)
   res <- sim_ode(
-    mod3,
+    mod_1cmt_iv_mm,
     parameters = par3,
     t_obs = t_obs,
     regimen = reg,
@@ -144,14 +146,14 @@ test_that("One compartment with MM kinetics", {
   )
 
   v1 <- get_var_y(
-    model = mod3,
+    model = mod_1cmt_iv_mm,
     parameters = par3,
     t_obs = t_obs,
     regimen = reg,
     omega = omega3
   )
   v2 <- get_var_y(
-    model = mod3,
+    model = mod_1cmt_iv_mm,
     parameters = par3,
     t_obs = t_obs,
     regimen = reg,
diff --git a/tests/testthat/test_iov.R b/tests/testthat/test_iov.R
index 650d138..e224100 100644
--- a/tests/testthat/test_iov.R
+++ b/tests/testthat/test_iov.R
@@ -1,44 +1,10 @@
-pars <- list(
-  "kappa_CL_1" = 0,
-  "kappa_CL_2" = 0,
-  "kappa_CL_3" = 0,
-  "eta_CL" = 0,
-  "CL" = 5,
-  "V" = 50,
-  "KA" = 1
-)
-pars0 <- list(
-  "CL" = 5,
-  "V" = 50,
-  "KA" = 1
-)
-pk0 <- new_ode_model( # no IOV
-  code = "
-      dAdt[1] = -KA * A[1]
-      dAdt[2] = +KA * A[1] -(CL/V) * A[2]
-    ",
-  obs = list(cmt = 2, scale = "V"),
-  dose = list(cmt = 1, bioav = 1),
-  parameters = names(pars0),
-  cpp_show_code = F
-)
-pk1 <- new_ode_model(
-  code = "
-      CL_iov = CL * exp(kappa_CL + eta_CL);
-      dAdt[1] = -KA * A[1]
-      dAdt[2] = +KA * A[1] -(CL_iov/V) * A[2]
-    ",
-  iov = list(
-    cv = list(CL = 0.2),
-    n_bins = 3
-  ),
-  obs = list(cmt = 2, scale = "V"),
-  dose = list(cmt = 1, bioav = 1),
-  declare_variables = c("kappa_CL", "CL_iov"),
-  parameters = names(pars),
-  cpp_show_code = F
-)
-reg1 <- new_regimen(
+# Uses models and parameters defined in setup.R:
+# - pk_iov_none (pk0)
+# - pk_iov (pk1)
+# - pars_iov
+# - pars_iov_no_iov
+
+reg_iov <- new_regimen(
   amt = 100,
   interval = 24,
   n = 5,
@@ -49,9 +15,9 @@ iov_var <- 0.3 ^ 2 # 30% IOV
 test_that("Throws error when `iov_bins` supplied but not present in model", {
   expect_error({
     sim(
-      ode = pk0,
-      parameters = pars0,
-      regimen = reg1,
+      ode = pk_iov_none,
+      parameters = pars_iov_no_iov,
+      regimen = reg_iov,
       omega = c(
         0.3 # IIV in CL
       ),
@@ -67,9 +33,9 @@ test_that("Throws error when `iov_bins` supplied but not present in model", {
 test_that("Throws error when number of `iov_bins` is higher than allowed for model", {
   expect_error({
     sim(
-      ode = pk1,
-      parameters = pars,
-      regimen = reg1,
+      ode = pk_iov,
+      parameters = pars_iov,
+      regimen = reg_iov,
       omega = c(
         iov_var, # IOV in CL
         0, iov_var,
@@ -89,9 +55,9 @@ test_that("Throws error when number of `iov_bins` is higher than allowed for mod
 test_that("Throws warning when number of `iov_bins` is lower than allowed for model", {
   expect_warning({
     sim(
-      ode = pk1,
-      parameters = pars,
-      regimen = reg1,
+      ode = pk_iov,
+      parameters = pars_iov,
+      regimen = reg_iov,
       omega = c(
         iov_var, # IOV in CL
         0, iov_var,
@@ -113,9 +79,9 @@ test_that("IOV is added to parameters", {
   set.seed(32)
 
   dat <- sim(
-    ode = pk1,
-    parameters = pars,
-    regimen = reg1,
+    ode = pk_iov,
+    parameters = pars_iov,
+    regimen = reg_iov,
     omega = c(
       iov_var, # IOV in CL
       0, iov_var,
@@ -149,17 +115,18 @@ test_that("IOV is added to parameters", {
 })
 
 test_that("Change in F in 2nd bin is applied in 2nd bin and not later.", {
+  skip_on_cran()
   # Previously this was an issue because F, when defined in pk_code(), was not updated before the
   # dose was applied to the state vector, so the bioavailability was not applied at the right time.
   # This was fixed by rearranging the order of execution in sim.cpp in the main loop.
 
-  pars <- list(
+  pars_iov_f <- list(
     "CL" = 5,
     "V" = 50,
     "KA" = 1,
     "F" = 1
   )
-  pk1 <- new_ode_model(
+  pk_iov_f <- new_ode_model(
     code = "
         dAdt[1] = -KA * A[1]
         dAdt[2] = +KA * A[1] -(CLi/V) * A[2]
@@ -175,18 +142,18 @@ test_that("Change in F in 2nd bin is applied in 2nd bin and not later.", {
     obs = list(cmt = 2, scale = "V"),
     dose = list(cmt = 1, bioav = "Fi"),
     declare_variables = c("kappa_F", "Fi", "CLi"),
-    parameters = names(pars),
-    cpp_show_code = F
+    parameters = names(pars_iov_f),
+    cpp_show_code = FALSE
   )
   reg <- new_regimen(amt = 800, interval = 24, n = 10, type = "oral")
 
   # For a first simulation, we're simulating with no variability across the IOV bins:
-  pars$kappa_F_1 <- 0
-  pars$kappa_F_2 <- 0
-  pars$kappa_F_3 <- 0
+  pars_iov_f$kappa_F_1 <- 0
+  pars_iov_f$kappa_F_2 <- 0
+  pars_iov_f$kappa_F_3 <- 0
   args_sim1 <- args <- list(
-    ode = pk1,
-    parameters = pars,
+    ode = pk_iov_f,
+    parameters = pars_iov_f,
     regimen = reg,
     only_obs = TRUE,
     t_obs = seq(0, 50, .25),
@@ -194,10 +161,10 @@ test_that("Change in F in 2nd bin is applied in 2nd bin and not later.", {
   )
   # For a second simulation, we're applying a change in parameter for the 2nd bin (24-48 hrs).
   # This should affect predictions from 24 onward.
-  pars$kappa_F_2 <- 1 # 2nd bin
+  pars_iov_f$kappa_F_2 <- 1 # 2nd bin
   args_sim2 <- args <- list(
-    ode = pk1,
-    parameters = pars,
+    ode = pk_iov_f,
+    parameters = pars_iov_f,
     regimen = reg,
     only_obs = TRUE,
     t_obs = seq(0, 50, .25),
@@ -217,9 +184,9 @@ test_that("error is not invoked when using parameters_table", {
   # specifying both parameters_table but for a model with IOV should not fail!
   expect_silent(
     dat <- sim(
-      ode = pk1,
+      ode = pk_iov,
       parameters_table = parameters_table,
-      regimen = reg1,
+      regimen = reg_iov,
       omega = c(
         iov_var, # IOV in CL
         0, iov_var,
@@ -238,10 +205,10 @@ test_that("error is not invoked when using parameters_table", {
   # specifying both parameters and parameters_table should fail
   expect_error(
     dat <- sim(
-      ode = pk1,
-      parameters = pars,
+      ode = pk_iov,
+      parameters = pars_iov,
       parameters_table = parameters_table,
-      regimen = reg1,
+      regimen = reg_iov,
       omega = c(
         iov_var, # IOV in CL
         0, iov_var,
diff --git a/tests/testthat/test_mixture_model.R b/tests/testthat/test_mixture_model.R
index 4289592..447239b 100644
--- a/tests/testthat/test_mixture_model.R
+++ b/tests/testthat/test_mixture_model.R
@@ -1,22 +1,15 @@
-## Setup model + params
-covs <- list(WT = PKPDsim::new_covariate(70))
-mod <- new_ode_model(
-  code = "
-      dAdt[0] = -(CL*(WT/70.0)/V)*A[0];
-    ",
-  pk_code = " ",
-  obs = list(cmt = 1, scale = "V"),
-  mixture = list(CL = list(values = c(5, 15), probability = 0.3)),
-  covariates = covs,
-)
-par <- list(CL = 3, V = 50)
-reg <- new_regimen(amt = 250, n = 5, interval = 6, type = 'infusion', t_inf = 1)
-t_obs <- seq(0, 36, 4)
+# Uses model and covariates defined in setup.R:
+# - mod_mixture
+# - covs_mixture
+
+par_mixture <- list(CL = 3, V = 50)
+reg_mixture <- new_regimen(amt = 250, n = 5, interval = 6, type = 'infusion', t_inf = 1)
+t_obs_mixture <- seq(0, 36, 4)
 
 test_that("mixture model works properly for single patient", {
-  res0 <- sim_ode(mod, parameters = par, regimen = reg, covariates = covs, t_obs = t_obs, only_obs=T) # mixture_group not supplied
-  res1 <- sim(mod, parameters = par, regimen = reg, t_obs = t_obs, covariates = covs, mixture_group = 1, only_obs=T)
-  res2 <- sim(mod, parameters = par, regimen = reg, t_obs = t_obs, covariates = covs, mixture_group = 2, only_obs=T)
+  res0 <- sim_ode(mod_mixture, parameters = par_mixture, regimen = reg_mixture, covariates = covs_mixture, t_obs = t_obs_mixture, only_obs=T) # mixture_group not supplied
+  res1 <- sim(mod_mixture, parameters = par_mixture, regimen = reg_mixture, t_obs = t_obs_mixture, covariates = covs_mixture, mixture_group = 1, only_obs=T)
+  res2 <- sim(mod_mixture, parameters = par_mixture, regimen = reg_mixture, t_obs = t_obs_mixture, covariates = covs_mixture, mixture_group = 2, only_obs=T)
   expect_equal(round(res0[res0$t == 24,]$y, 2), 9.07) # should use whatever is in `parameters`
   expect_equal(round(res1[res1$t == 24,]$y, 2), 5.82)
   expect_equal(round(res2[res2$t == 24,]$y, 2), 1.15)
@@ -25,9 +18,9 @@ test_that("mixture model works properly for single patient", {
 test_that("mixture model works properly when vectorized (using parameters_table)", {
   partab <- data.frame(CL = rep(0, 6), V = rep(50, 6))
   suppressMessages({
-    expect_error(sim_ode(mod, parameters_table = partab, regimen = reg, t_obs = t_obs, covariates = covs, mixture_group = 1, only_obs=T))
-    res1 <- sim(mod, parameters_table = partab, regimen = reg, t_obs = t_obs, covariates = covs, mixture_group = rep(1, 6), only_obs=T)
-    res2 <- sim(mod, parameters_table = partab, regimen = reg, t_obs = t_obs, covariates = covs, mixture_group = rep(c(1,2), 3), only_obs=T, output_include = list(parameters = TRUE))
+    expect_error(sim_ode(mod_mixture, parameters_table = partab, regimen = reg_mixture, t_obs = t_obs_mixture, covariates = covs_mixture, mixture_group = 1, only_obs=T))
+    res1 <- sim(mod_mixture, parameters_table = partab, regimen = reg_mixture, t_obs = t_obs_mixture, covariates = covs_mixture, mixture_group = rep(1, 6), only_obs=T)
+    res2 <- sim(mod_mixture, parameters_table = partab, regimen = reg_mixture, t_obs = t_obs_mixture, covariates = covs_mixture, mixture_group = rep(c(1,2), 3), only_obs=T, output_include = list(parameters = TRUE))
   })
   expect_equal(round(res1[res1$t == 24,]$y, 2), rep(5.82, 6))
   expect_equal(round(res2[res2$t == 24,]$y, 2), rep(c(5.82, 1.15), 3))
@@ -39,8 +32,8 @@ test_that("mixture model works properly when vectorized (using parameters_table)
 test_that("mixture model works properly when vectorized (using covariates_table)", {
   covtab <- data.frame(ID = 1:8, WT = rep(seq(40, 130, 30), 2))
   suppressMessages({
-    expect_error(sim(mod, parameters = par, covariates_table = covtab, regimen = reg, t_obs = t_obs, mixture_group = 1, only_obs=T))
-    res <- sim(mod, parameters = par, covariates_table = covtab, regimen = reg, t_obs = t_obs, mixture_group = rep(c(1, 2), each=4), only_obs=T)
+    expect_error(sim(mod_mixture, parameters = par_mixture, covariates_table = covtab, regimen = reg_mixture, t_obs = t_obs_mixture, mixture_group = 1, only_obs=T))
+    res <- sim(mod_mixture, parameters = par_mixture, covariates_table = covtab, regimen = reg_mixture, t_obs = t_obs_mixture, mixture_group = rep(c(1, 2), each=4), only_obs=T)
   })
   expect_equal(round(res[res$t == 24,]$y, 2), c(9.39, 5.82, 3.83, 2.65, 2.99, 1.15, 0.52, 0.25))
 })
diff --git a/tests/testthat/test_multi_obs.R b/tests/testthat/test_multi_obs.R
index f64511f..15be190 100644
--- a/tests/testthat/test_multi_obs.R
+++ b/tests/testthat/test_multi_obs.R
@@ -1,24 +1,17 @@
-# Example multiple observation types (e.g. different compartments! not just different residual errors)
+# Uses model defined in setup.R:
+# - pk_multi_obs
+# - vars_multi_obs
 
-## define parameters
-vars <- c("CONC", "METAB", "METAB2", "ACT")
-pk1 <- new_ode_model(
-  code = "dAdt[1] = -(CL/V)*A[1]; CONC = 1000*A[1]/V; METAB = CONC/2; METAB2 = CONC * t; ACT = 15",
-  obs = list(variable = vars),
-  declare_variables = vars,
-  cpp_show_code = F
-)
-
-regimen  <- new_regimen(amt = 100, interval = 12, n = 5, type="infusion", t_inf = 1)
-parameters   <- list("CL" = 15, "V" = 150)
-omega <- PKPDsim::cv_to_omega(list("CL" = 0.2, "V" = 0.2), parameters[1:2])
+regimen_multi_obs <- new_regimen(amt = 100, interval = 12, n = 5, type="infusion", t_inf = 1)
+parameters_multi_obs <- list("CL" = 15, "V" = 150)
+omega_multi_obs <- PKPDsim::cv_to_omega(list("CL" = 0.2, "V" = 0.2), parameters_multi_obs[1:2])
 
 test_that("obs types are output by `sim`", {
   obs_type <- c(1,2,1,3,1)
   data <- sim(
-    ode = pk1,
+    ode = pk_multi_obs,
     parameters = list(CL = 20, V = 200),
-    regimen = regimen,
+    regimen = regimen_multi_obs,
     int_step_size = 0.1,
     only_obs = TRUE,
     obs_type = obs_type,
@@ -33,9 +26,9 @@ test_that("obs types are output by `sim`", {
 test_that("check obs at same timepoint but with different obs_type", {
   obs_type <- c(1,2,1,3,1)
   t_same <- sim(
-    ode = pk1,
+    ode = pk_multi_obs,
     parameters = list(CL = 20, V = 200),
-    regimen = regimen,
+    regimen = regimen_multi_obs,
     int_step_size = 0.1,
     only_obs = TRUE,
     obs_type = obs_type,
@@ -58,9 +51,9 @@ test_that("check that residual error correctly applied to right var", {
   ruv_term1 <- list(prop = c(.1, 0, 0), add = c(1, 0, 0))
 
   data2 <- sim(
-    ode = pk1,
+    ode = pk_multi_obs,
     parameters = list(CL = 20, V = 200),
-    regimen = regimen,
+    regimen = regimen_multi_obs,
     int_step_size = 0.1,
     only_obs = TRUE,
     obs_type = obs_type,
@@ -75,9 +68,9 @@ test_that("check that residual error correctly applied to right var", {
   expect_true(data2$y[-4][4] != y[4])
 
   data3 <- sim(
-    ode = pk1,
+    ode = pk_multi_obs,
     parameters = list(CL = 20, V = 200),
-    regimen = regimen,
+    regimen = regimen_multi_obs,
     int_step_size = 0.1,
     only_obs = TRUE,
     obs_type = obs_type,
@@ -99,10 +92,10 @@ test_that("specifying ruv as multi-type when only 1 obs_type", {
   ruv_multi <- list(prop = c(0.1, 1), add = c(1, 20))
 
   s_single <- sim(
-    ode = pk1,
-    parameters = parameters,
+    ode = pk_multi_obs,
+    parameters = parameters_multi_obs,
     n_ind = 1,
-    regimen = regimen,
+    regimen = regimen_multi_obs,
     only_obs = TRUE,
     t_obs = c(2,4,6,8),
     seed = 123456,
@@ -112,10 +105,10 @@ test_that("specifying ruv as multi-type when only 1 obs_type", {
   ## specified as multi, but obs_type is all 1, so should
   ## give same results as s_single
   s_multi1 <- sim(
-    ode = pk1,
-    parameters = parameters,
+    ode = pk_multi_obs,
+    parameters = parameters_multi_obs,
     n_ind = 1,
-    regimen = regimen,
+    regimen = regimen_multi_obs,
     only_obs = TRUE,
     t_obs = c(2,4,6,8),
     obs_type = c(1, 1, 1, 1),
@@ -128,10 +121,10 @@ test_that("specifying ruv as multi-type when only 1 obs_type", {
   ## specifying ruv as multi-type when multiple obs_type
   ## should now give different results
   s_multi1 <- sim(
-    ode = pk1,
-    parameters = parameters,
+    ode = pk_multi_obs,
+    parameters = parameters_multi_obs,
     n_ind = 1,
-    regimen = regimen,
+    regimen = regimen_multi_obs,
     only_obs = TRUE,
     t_obs = c(2,4,6,8),
     obs_type = c(1, 2, 1, 2),
@@ -145,9 +138,9 @@ test_that("specifying ruv as multi-type when only 1 obs_type", {
 
 test_that("multi-obs with baseline and obs_time = dose_time works correctly", {
   tmp <- sim(
-    pk1,
-    parameters = parameters,
-    regimen = regimen,
+    pk_multi_obs,
+    parameters = parameters_multi_obs,
+    regimen = regimen_multi_obs,
     t_obs = c(0, 0, 6, 6),
     obs_type = c(1, 4, 1, 4),
     only_obs = TRUE,
diff --git a/tests/testthat/test_parameters_table.R b/tests/testthat/test_parameters_table.R
index 1f9af03..e274d7c 100644
--- a/tests/testthat/test_parameters_table.R
+++ b/tests/testthat/test_parameters_table.R
@@ -5,19 +5,11 @@ test_that("Simulating with table of parameters works", {
     V = rnorm(10, 5, 0.5)
   )
 
-  mod <- new_ode_model(
-    code = "
-      dAdt[1] = -(CL/V) * A[1];
-    ",
-    obs = list(cmt = 1, scale = "V"),
-    covariates = list(WT = new_covariate(70)),
-    cpp_show_code=FALSE
-  )
   par <- list(CL = 5, V = 50)
   reg <- new_regimen(amt = 2000, interval = 24, type = "infusion")
   covariates = list(WT = new_covariate(1))
   res <- sim_ode(
-    ode = mod,
+    ode = oral_1cmt_allometric,
     parameters_table = parameters_table,
     covariates = covariates,
     regimen = reg,
diff --git a/tests/testthat/test_reparametrization.R b/tests/testthat/test_reparametrization.R
index 747d373..f3352de 100644
--- a/tests/testthat/test_reparametrization.R
+++ b/tests/testthat/test_reparametrization.R
@@ -1,3 +1,9 @@
+# Uses model and covariates defined in setup.R (conditional, NOT_CRAN only):
+# - model_carreno
+# - covs_carreno
+
+skip_on_cran()
+
 par_orig <- list(
   V = 25.76,
   SCLSlope = 0.036,
@@ -6,33 +12,8 @@ par_orig <- list(
   SCLInter = 0.18,
   TDM_INIT = 0)
 
-covs <- list(
-  CRCL = PKPDsim::new_covariate(5),
-  CL_HEMO = PKPDsim::new_covariate(0)
-)
-model <- new_ode_model( # Carreno et al
-  code = "
-    CLi = SCLInter + SCLSlope * (CRCL*16.667) + CL_HEMO \
-    Vi = V \
-    Qi = K12 * Vi \
-    V2i = Qi / K21 \
-    dAdt[0] = -(CLi/V)*A[0] - K12*A[0] + K21*A[1] \
-    dAdt[1] = K12*A[0] - K21*A[1] \
-    dAdt[2] = A[0]/V
-  ",
-  parameters = par_orig,
-  declare_variables = c("CLi", "Qi", "Vi", "V2i"),
-  covariates= covs,
-  obs = list(cmt = 1, scale = "V"),
-  reparametrization = list(
-    "CL" = "SCLInter + SCLSlope * (CRCL*16.667) + CL_HEMO",
-    "V" = "V",
-    "Q" = "K12 * V",
-    "V2" = "(K12 * V) / K21"
-  )
-)
-pars_covs_comb <- join_cov_and_par(covs, par_orig)
-pars <- reparametrize(model, pars_covs_comb, covariates = covs)
+pars_covs_comb <- join_cov_and_par(covs_carreno, par_orig)
+pars <- reparametrize(model_carreno, pars_covs_comb, covariates = covs_carreno)
 
 reg <- new_regimen(
   amt = 1000,
@@ -44,10 +25,10 @@ reg <- new_regimen(
   n_ss = 50
 )
 s <- sim(
-  ode = model,
+  ode = model_carreno,
   parameters = par_orig,
   regimen = reg,
-  covariates = covs,
+  covariates = covs_carreno,
   t_obs = c(0,24)
 )
 
@@ -96,4 +77,3 @@ test_that("Reparametrized model and analytics equations match", {
   cmin_ss_ode <- s[s$comp == "obs",]$y[1]
   expect_equal(cmin_ss_ode, cmin_ss_lin)
 })
-
diff --git a/tests/testthat/test_t_init.R b/tests/testthat/test_t_init.R
index f0e1a63..abd4e2e 100644
--- a/tests/testthat/test_t_init.R
+++ b/tests/testthat/test_t_init.R
@@ -1,23 +1,19 @@
-# Test t_init functionality
+# Uses model defined in setup.R (conditional, NOT_CRAN only):
+# - mod_t_init
+
+skip_on_cran()
 
+# Test t_init functionality
 ## e.g. TDM before first dose:
 ## at t=-8, conc=10000
 ## Use this as true value in the simulations
-par   <- list(CL = 7.67, V = 97.7, TDM_INIT = 500)
-mod <- new_ode_model(
-  code = "CLi = CL; Vi = V; dAdt[1] = -(CLi/Vi)*A[1]; CONC = A[1]/Vi",
-  state_init = "A[1] = TDM_INIT * Vi",
-  parameters = par,
-  obs = list(cmt = 1, scale = "Vi"),
-  declare_variables = c("CONC", "CLi", "Vi"),
-  cpp_show_code = F
-)
-reg <- new_regimen(amt = 100000, times=c(0, 24), type="bolus")
+par_t_init <- list(CL = 7.67, V = 97.7, TDM_INIT = 500)
+reg_t_init <- new_regimen(amt = 100000, times=c(0, 24), type="bolus")
 s <- sim(
-  ode = mod,
-  parameters = par,
+  ode = mod_t_init,
+  parameters = par_t_init,
   n_ind = 1,
-  regimen = reg,
+  regimen = reg_t_init,
   only_obs = TRUE,
   output_include = list(variables = TRUE),
   t_init = 10
@@ -33,5 +29,3 @@ test_that("TDM before first dose is considered a true initial value", {
 test_that("Variables are set (also in first row) when TDM before first dose", {
   expect_equal(round(s$CONC[1:5], 1), c(500, 462.2, 427.3, 395.1, 365.3))
 })
-
-
diff --git a/tests/testthat/test_timevar_cov.R b/tests/testthat/test_timevar_cov.R
index 60866db..b9b9248 100644
--- a/tests/testthat/test_timevar_cov.R
+++ b/tests/testthat/test_timevar_cov.R
@@ -1,23 +1,16 @@
+# Uses model defined in setup.R (conditional, NOT_CRAN only):
+# - mod_2cmt_timevar
+
 skip_on_cran()
 
-par <- list(CL = 3, V = 50, Q = 2.5, V2 = 70)
-reg <- new_regimen(
+par_timevar <- list(CL = 3, V = 50, Q = 2.5, V2 = 70)
+reg_timevar <- new_regimen(
   amt = 250,
   n = 60,
   interval = 6,
   type = 'infusion',
   t_inf = 1
 )
-mod <- new_ode_model(
-  code = "
-      dAdt[1] = -(Q/V)*A[1] + (Q/V2)*A[2] -(CLi/V)*A[1];
-      dAdt[2] = -(Q/V2)*A[2] + (Q/V)*A[1];
-    ",
-  pk_code = "CLi = CL + CRCL",
-  obs = list(cmt = 2, scale = "V"),
-  covariates = list(CRCL = new_covariate(5)), declare_variables = "CLi",
-  cpp = FALSE
-)
 
 test_that("timevarying covariates handled", {
   # CLi changes by several orders of magnitude after
@@ -32,9 +25,9 @@ test_that("timevarying covariates handled", {
   )
   t_obs <- seq(0, 360, 0.1)
   sim1 <- sim_ode(
-    mod,
-    parameters = par,
-    regimen = reg,
+    mod_2cmt_timevar,
+    parameters = par_timevar,
+    regimen = reg_timevar,
     covariates = covs,
     only_obs = TRUE,
     t_obs = t_obs,
@@ -62,18 +55,18 @@ test_that("timevarying covariates are interpolated and affect PK", {
   )
   t_obs <- seq(0, 120, 2)
   sim2_inter <- sim_ode(
-    mod,
-    parameters = par,
-    regimen = reg,
+    mod_2cmt_timevar,
+    parameters = par_timevar,
+    regimen = reg_timevar,
     covariates = covs_inter,
     only_obs = TRUE,
     t_obs = t_obs,
     output_include = list(parameters = TRUE, covariates = TRUE, variables = TRUE)
   )
   sim2_locf <- sim_ode(
-    mod,
-    parameters = par,
-    regimen = reg,
+    mod_2cmt_timevar,
+    parameters = par_timevar,
+    regimen = reg_timevar,
     covariates = covs_locf,
     only_obs = TRUE,
     t_obs = t_obs,