RXR-2394: lagtime updates

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: PKPDsim
 Type: Package
 Title: Tools for Performing Pharmacokinetic-Pharmacodynamic Simulations
-Version: 1.4.1
-Date: 2025-04-09
+Version: 1.5.0
+Date: 2025-07-08
 Authors@R: c(
         person("Ron", "Keizer", email = "ron@insight-rx.com", role = c("aut", "cre")),
         person("Jasmine", "Hughes", email = "jasmine@insight-rx.com", role = "aut"),

NAMESPACE

@@ -15,7 +15,6 @@ export(advan_create_data)
 export(advan_parse_output)
 export(advan_process_infusion_doses)
 export(apply_duration_scale)
-export(apply_lagtime)
 export(available_default_literature_models)
 export(calc_auc_analytic)
 export(calc_ss_analytic)

R/compile_sim_cpp.R

@@ -270,10 +270,18 @@ compile_sim_cpp <- function(
   if(length(grep("-w", flg)) == 0) {
     Sys.setenv("PKG_CXXFLAGS" = paste(flg, "-w"))
   }
+
   if(compile) {
-    Rcpp::sourceCpp(code=sim_func, rebuild = TRUE, env = globalenv(), verbose = verbose, showOutput = verbose)
+    Rcpp::sourceCpp(
+      code = sim_func,
+      rebuild = TRUE, 
+      env = globalenv(), 
+      verbose = verbose, 
+      showOutput = verbose
+    )
     Sys.setenv("PKG_CXXFLAGS" = flg)
   }
+
   return(list(
     ode = ode_def_cpp,
     cpp = sim_func

R/create_event_table.R

@@ -152,10 +152,9 @@ create_event_table <- function(
   design <- dos[order(dos$t, -dos$dose),]
   if(!is.null(t_obs) && length(t_obs) != 0) { # make sure observation times are in dataset
     t_obs <- round(t_obs, 6)
-    t_diff <- setdiff(t_obs, design$t)
-    if(length(t_diff) > 0) {
-      design[(length(design[,1])+1) : (length(design[,1])+length(t_diff)),] <- cbind(
-         t = t_diff,
+    if(length(t_obs) > 0) {
+      design[(length(design[,1])+1) : (length(design[,1])+length(t_obs)),] <- cbind(
+         t = t_obs,
          dose = 0,
          type = 0,
          dum = 0,

R/parse_lagtime.R

@@ -0,0 +1,27 @@
+#' Parse lagtime specified to main sim() function
+#' 
+#' @inheritParams sim
+#' 
+#' @returns a vector of character or numeric values
+#' 
+parse_lagtime <- function(
+  lagtime,
+  ode,
+  parameters
+) {
+  lagtime_ode <- attr(ode, "lagtime")
+  # override from ode if not specified by user and defined in ode
+  if(is.null(lagtime) && !is.null(lagtime_ode) && lagtime_ode[1] != "NULL" && lagtime_ode[1] != "undefined") {
+    lagtime <- lagtime_ode
+  }
+  if(is.null(lagtime)) {
+    lagtime <- 0
+  }
+  if(inherits(lagtime, "character")) {
+    idx <- grep("[a-zA-Z]", lagtime) # only pick character names, not "0"
+    if(! all(lagtime[idx] %in% names(parameters))) {
+      warning("Lagtime parameter(s) not found. Please check model and parameters.")
+    }
+  } 
+  lagtime
+}

R/sim.R

@@ -147,14 +147,7 @@ sim <- function (ode = NULL,
       regimen <- merge_regimen(list(regimen, regimen_dupl))
     }
   }
-  if(!is.null(attr(ode, "lagtime")) && attr(ode, "lagtime")[1] != "undefined" && attr(ode, "lagtime")[1] != "NULL") {
-    if(is.null(lagtime)) { # only override from metadata if not specified by user
-      lagtime <- attr(ode, "lagtime")
-    }
-  }
-  if(!is.null(lagtime)) {
-    regimen <- apply_lagtime(regimen, lagtime, parameters, attr(ode, "cmt_mapping"))
-  }
+  lagtime <- parse_lagtime(lagtime, ode, parameters)
   if(!is.null(attr(ode, "dose")$duration_scale)) {
     regimen <- apply_duration_scale(
       regimen,
@@ -234,6 +227,9 @@ sim <- function (ode = NULL,
     } else {
       size <- attr(analytical, "size")
     }
+    if(is.null(lagtime)) {
+      lagtime <- rep(0, size) # needs to have at least 1 zero value, cannot be NULL when passed to cpp func
+    }
     if(is.null(ode) && is.null(analytical)) {
       stop("Please specify at least the required arguments 'ode' or 'analytical' for simulations.")
     }
@@ -507,7 +503,7 @@ sim <- function (ode = NULL,
 
     #################### Main call to ODE solver / analytical eq solver #######################
     if(!is.null(ode)) {
-      tmp <- ode(A_init, design_i, p_i, iov_bins, int_step_size)
+      tmp <- ode(A_init, design_i, p_i, iov_bins, lagtime, int_step_size)
     } else {
       tmp <- analytical_eqn_wrapper(analytical, design_i, p_i)
     }
@@ -602,6 +598,12 @@ sim <- function (ode = NULL,
   all_names <- unique(c(par_names, cov_names, var_names))
   all_names <- intersect(all_names, names(comb)) # only cols that appear in data
 
+  ## remove dose-times from regimen
+  ## but leave ones that were also in t_obs
+  ## also leave extra obs when bolus, because this may be needed (controlled below using `extra_t_obs`)
+  dose_times <- design_i[design_i$evid > 0, ]$t
+  comb <- comb[! (comb$t %in% dose_times & duplicated(paste(comb$id, comb$comp, comb$t, comb$obs_type, comb$y, sep="_"))),]
+
   if(!extra_t_obs) {
     ## include the observations at which a bolus dose is added into the output object too
     comb <- comb[!duplicated(paste(comb$id, comb$comp, comb$t, comb$obs_type, sep="_")),]

R/sim_core.R

@@ -1,22 +1,28 @@
 #' Only core function of the simulation function, always just returns observations.
 #' Mostly useful for estimations / optimal design. Has no checks (for speed)!
 #'
+#' @inheritParams sim
 #' @param sim_object list with design and simulation parameters
-#' @param ode ode
-#' @param duplicate_t_obs allow duplicate t_obs in output? E.g. for optimal design calculations when t_obs = c(0,1,2,2,3). Default is FALSE.
-#' @param t_init time of initialization of the ODE system. Usually 0.
+#' 
 #' @export
-#' @return Data frame with simulation results
+#' 
+#' @return data.frame with simulation results
+#' 
 sim_core <- function(
   sim_object = NULL,
   ode,
   duplicate_t_obs = FALSE,
-  t_init = 0) {
-  tmp <- ode(A = sim_object$A_init,
-             design = sim_object$design,
-             par = sim_object$p,
-             iov_bins = sim_object$iov_bins,
-             step_size = sim_object$int_step_size)
+  t_init = 0,
+  lagtime = c(0)
+) {
+  tmp <- ode(
+    A = sim_object$A_init,
+    design = sim_object$design,
+    par = sim_object$p,
+    iov_bins = sim_object$iov_bins,
+    lagtime = lagtime,
+    step_size = sim_object$int_step_size
+  )
   out <- data.frame(t = tmp$time, y = tmp$obs, obs_type = tmp$obs_type)
   if(duplicate_t_obs) {
     # use match to ensure that duplicates in t_obs is possible

inst/cpp/sim.cpp

@@ -63,23 +63,117 @@ void pk_code (int i, std::vector<double> times, std::vector<double> doses, doubl
   // insert custom pk event code
 }
 
+NumericVector lagtime_to_numeric(SEXP lagtime, List parameters) {
+   NumericVector lagtime_numeric;
+  if (TYPEOF(lagtime) == REALSXP) {
+    lagtime_numeric = as<NumericVector>(lagtime);
+  } else if (TYPEOF(lagtime) == STRSXP) {
+    CharacterVector lagtime_char = as<CharacterVector>(lagtime);
+    lagtime_numeric = NumericVector(lagtime_char.size());
+    for (int i = 0; i < lagtime_char.size(); i++) {
+      String param_name = lagtime_char[i];
+      if (parameters.containsElementNamed(param_name.get_cstring())) {
+        lagtime_numeric[i] = as<double>(parameters[param_name]);
+      } else {
+        lagtime_numeric[i] = 0.0; // default value
+      }
+    }
+  } else {
+    stop("lagtime must be either numeric or character vector");
+  }
+  return(lagtime_numeric);
+}
+
+List apply_lagtime(List design, NumericVector tlag, int n_comp) {
+
+  List new_design = clone(design);
+  std::vector<double> times = as<std::vector<double> >(new_design["t"]);
+  std::vector<int> evid = as<std::vector<int> >(new_design["evid"]);
+  std::vector<int> rate = as<std::vector<int> >(new_design["rate"]);
+  std::vector<int> cmt = as<std::vector<int> >(new_design["dose_cmt"]);
+
+  NumericVector lagtime;
+  if(tlag.size() < n_comp) { // fill in with zeroes, if needed
+    int current_size = tlag.size();
+    lagtime = NumericVector(n_comp);
+    for(int i = 0; i < current_size; i++) {
+      lagtime[i] = tlag[i];
+    }
+    for(int i = current_size; i < n_comp; i++) {
+      lagtime[i] = 0.0;
+    }
+  } else {
+    lagtime = tlag;
+  }
+
+  // Apply lagtime to dose events (evid == 1)  
+  for(int i = 0; i < times.size(); i++) {
+    if(evid[i] == 1 | (evid[i] == 2 & rate[i] != 0)) { // dose events or infusion stop events
+      times[i] += lagtime[cmt[i]-1];
+    }
+  }
+
+  // Create sorted index according to "t"
+  std::vector<size_t> indices(times.size());
+  std::iota(indices.begin(), indices.end(), 0); // Fill with 0, 1, 2, ...
+  std::sort(indices.begin(), indices.end(), [&times](size_t i1, size_t i2) {
+    return times[i1] < times[i2];
+  });
+
+  // Reorder all elements in `t` in new_design
+  std::vector<double> sorted_times(times.size());
+  for (size_t i = 0; i < indices.size(); i++) {
+    sorted_times[i] = times[indices[i]];
+  }
+  new_design["t"] = sorted_times;
+
+  // Sort all other vectors in the design object
+  for (const char* key : {"dose", "type", "dum", "dose_cmt", "t_inf", "evid", "bioav", "rate", "obs_type"}) {
+    if (new_design.containsElementNamed(key)) {
+      SEXP vec = new_design[key];
+      if (TYPEOF(vec) == REALSXP) {
+        std::vector<double> old_vec = as<std::vector<double> >(vec);
+        std::vector<double> new_vec(old_vec.size());
+        for (size_t i = 0; i < indices.size(); i++) {
+          new_vec[i] = old_vec[indices[i]];
+        }
+        new_design[key] = new_vec;
+      } else if (TYPEOF(vec) == INTSXP) {
+        std::vector<int> old_vec = as<std::vector<int> >(vec);
+        std::vector<int> new_vec(old_vec.size());
+        for (size_t i = 0; i < indices.size(); i++) {
+          new_vec[i] = old_vec[indices[i]];
+        }
+        new_design[key] = new_vec;
+      }
+    }
+  }
+
+  return new_design;
+}
+
 // [[Rcpp::export]]
-List sim_wrapper_cpp (NumericVector A, List design, List par, NumericVector iov_bins, double step_size) {
+List sim_wrapper_cpp (NumericVector A, List design, List par, NumericVector iov_bins, SEXP lagtime, double step_size) {
   std::vector<double> t;
   std::vector<state_type> y;
   // insert observation variable definition
   double t_start, t_end;
   std::vector<double> times, doses, dummy, rates;
   std::vector<int> dose_cmt, dose_type, evid, obs_type, y_type;
   // insert variable definitions
-  times = as<std::vector<double> >(design["t"]);
-  doses = as<std::vector<double> >(design["dose"]);
-  evid = as<std::vector<int> >(design["evid"]);
-  dummy = as<std::vector<double> >(design["dum"]);
-  rates = as<std::vector<double> >(design["rate"]);
-  dose_cmt = as<std::vector<int> >(design["dose_cmt"]);
-  dose_type = as<std::vector<int> >(design["type"]);
-  obs_type = as<std::vector<int> >(design["obs_type"]);
+
+  // Handle lagtime parameter - can be numeric or character
+  NumericVector lagtime_numeric = lagtime_to_numeric(lagtime, par);
+  List events = apply_lagtime(design, lagtime_numeric, n_comp);
+
+  times = as<std::vector<double> >(events["t"]);
+  doses = as<std::vector<double> >(events["dose"]);
+  evid = as<std::vector<int> >(events["evid"]);
+  dummy = as<std::vector<double> >(events["dum"]);
+  rates = as<std::vector<double> >(events["rate"]);
+  dose_cmt = as<std::vector<int> >(events["dose_cmt"]);
+  dose_type = as<std::vector<int> >(events["type"]);
+  obs_type = as<std::vector<int> >(events["obs_type"]);
   int len = times.size();
   int start;
   memset(rate, 0, sizeof(rate));