Feature proposition : automatically define $ODE for linear models

[FelicienLL]

Hi @kylebaron,

For a project, I am modelling with the ADVAN5/7 routines in NONMEM, which are convenient for linear models because there are no ODE to write, just to define the rate constants in $PK. But I wanted to transpose my model to mrgsolve, so equations had to be written eventually... So I wrote a little bit of code to automatically achieve that.

I thought it could have been a nice feature to implement into mrgsolve, for instance in the form of a plugin that would automatically define the $ODE when the model is linear... I don't know if this was something you had already considered? In that case, maybe could this code be helpful? Otherwise, please ignore.

Félicien

# Linear model, i.e. in NONMEM, corresponds to ADVAN5 or ADVAN7
# The idea : only write the rate constants in the control stream ($MAIN, $PK)
# The $ODE get automatically constructed as function of rate constants defined and $CMT defined

# Two possible syntaxes for rate constants (see NONMEM documentation)

# 1 = K12, K13...
ksyntax1 <- "^K{1}\\d{2}$"
# 2 = K1T2, K11T12, useful if n compartments > 10, and avoid ambiguous K like K121... 
ksyntax2 <- "^K{1}\\d+T{1}\\d+$"

grepl(ksyntax1, c("K12", "k12", "K121", "K12bis", "CL", "KK12", "K10T11"))
#> [1]  TRUE FALSE FALSE FALSE FALSE FALSE FALSE
grepl(ksyntax2, c("K1T2","K11T2","K1T22","K11T22", "k1T2", "KK10T11", "K1TT1", "K1T1bis", "K12"))
#> [1]  TRUE  TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE

keep_k <- function(p){
  #p, a vector a parameter names
  
  p[c(grep(ksyntax1, p), grep(ksyntax2, p))]
  
  # extracts only the parameters identified as rate constants, to implement in $ODE
}
keep_k(p = c("CL", "AGE", "KA", "K10", "K21", "K15T24", "K1T0"))
#> [1] "K10"    "K21"    "K15T24" "K1T0"

fromk <- function(k){
  #k, a vector of rate constants 
  
  #seems we cannot vectorize through ifelse because regmatches drops missing matches
  which1 <- which(grepl(ksyntax1, k))
  from1 <- substr(k[which1], 2, 2)
  
  which2 <- which(grepl(ksyntax2, k))
  from2 <- regmatches(k[which2], m = regexpr(pattern = "\\d+(?=T\\d+$)", text = k[which2], perl = TRUE))
  
  from <- character(length(k))
  from[which1] <- from1
  from[which2] <- from2
  from
  # returns the compartment number it starts from 
}

tok <- function(k){
  #k, a vector of rate constants 
  
  #seems we cannot vectorize through ifelse because regmatches drops missing matches
  which1 <- which(grepl(ksyntax1, k))
  to1 <- substr(k[which1], 3, 3)
  
  which2 <- which(grepl(ksyntax2, k))
  to2 <- regmatches(k[which2], m = regexpr(pattern = "\\d+$", text = k[which2], perl = TRUE))
  
  to <- character(length(k))
  to[which1] <- to1
  to[which2] <- to2
  to
  
  # returns the compartment number it ends to 
}

fromk(c("K12", "K1T21", "K25", "K12T22"))
#> [1] "1"  "1"  "2"  "12"
tok(c("K12", "K1T21", "K25", "K12T22"))
#> [1] "2"  "21" "5"  "22"

validate_k <- function(k){
  #k, a vector of rate constants 
  from <- fromk(k)
  to <- tok(k)
  
  if(any(from == "0")) stop("Error in K, cannot starts with 0")
  if(any(from == to)) stop("Error in K, `from` is the same than `to`")
  if(any(duplicated(paste0("K", from, "T", to)))) stop("Duplicates of K were found")
  
}
# some tests
validate_k("K01")
#> Error in validate_k("K01"): Error in K, cannot starts with 0
validate_k("K11")
#> Error in validate_k("K11"): Error in K, `from` is the same than `to`
validate_k(c("K10", "K1T0"))
#> Error in validate_k(c("K10", "K1T0")): Duplicates of K were found

validate_k_cmts <- function(k, cmts){
  #k, a vector of rate constants 
  #cmts, a vector of all compartments defined in the model
  from <- fromk(k)
  to <- tok(k)
  cmtsnum <- seq_along(cmts)
  
  if(!all(from %in% c(0, cmtsnum))) stop("One value of K was defined with an invalid `from` compartment number")
  if(!all(to %in% c(0, cmtsnum))) stop("One value of K was defined with an invalid `to` compartment number")
}

# some more tests
validate_k_cmts("K20", "GUT")
#> Error in validate_k_cmts("K20", "GUT"): One value of K was defined with an invalid `from` compartment number
validate_k_cmts("K12", "GUT")
#> Error in validate_k_cmts("K12", "GUT"): One value of K was defined with an invalid `to` compartment number

lhs <- function(i, cmts){
  #i, the number of the compartment for which we want to write the equation
  #cmts, a vector of all compartments defined in the model
  paste0("dxdt_", cmts[i])
  
  #retuns the left-hand side of the equation
}

lhs(1, c("CENTRAL", "PERIPH"))
#> [1] "dxdt_CENTRAL"
lhs(2, c("CENTRAL", "PERIPH"))
#> [1] "dxdt_PERIPH"

removeplus <- function(x) gsub(pattern = "^ ?\\+? ?", replacement =  "", x = x)
# quick helper to avoid equation starting with a sign "+"
removeplus(c(" + K12 * GUT", "- K10 * CENTRAL + K21 * PERIPH"))
#> [1] "K12 * GUT"                      "- K10 * CENTRAL + K21 * PERIPH"

rhs <- function(i, cmts, k){
  #i, the number of the compartment for which we want to write the equation
  #cmts, a vector of all compartments defined in the model
  #k, a vector of rate constants 
  cmtsnum <- seq_along(cmts)
  
  getout <- ""
  from <- fromk(k)
  ifrom <- from == i
  if(any(ifrom)){
    if(length(which(ifrom)) == 1){
      getout <- paste0("- ", k[ifrom], " * ", cmts[i])
    } else {
      getout <- paste0("- (", paste(k[ifrom], collapse = " + "), ") * ", cmts[i])
    }
  }
  
  getin <- ""
  to <- tok(k)
  ito <- to == i
  if(any(ito)){
    getin <- paste("+", (paste(k[ito], cmts[which(cmtsnum %in% from[ito])], sep = " * ")))
  }
  
  if(!any(ito, ifrom)) return("0 ;")
  
  ans <- paste(
    paste(getout, collapse = " "), 
    paste(getin, collapse = " "),
    ";",
    collapse = " "
  )
  
  ans <- removeplus(ans)
  
  ans
  #retuns the right-hand side of the equation
}

rhs(1, c("CENTRAL", "PERIPH"), c("K10", "K21", "K12"))
#> [1] "- (K10 + K12) * CENTRAL + K21 * PERIPH ;"
rhs(2, c("CENTRAL", "PERIPH"), c("K10", "K21", "K12"))
#> [1] "- K21 * PERIPH + K12 * CENTRAL ;"
rhs(1, c("GUT", "CENT"), c("K12"))
#> [1] "- K12 * GUT  ;"
rhs(2, c("GUT", "CENT"), c("K12"))
#> [1] "K12 * GUT ;"
rhs(2, c("CENTRAL", "PERIPH"), c("K10"))
#> [1] "0 ;"



# Final function to do it all ----

linearmodel <- function(cmts, k, .cat = FALSE){
  #cmts, a vector of all compartments defined in the model
  #k, a vector of rate constants 
  #.cat, a logical, should the result be printed to the console with cat
  
  validate_k(k = k)
  validate_k_cmts(k = k, cmts = cmts)
  cmtsnum <- seq_along(cmts)
  
  ans <- paste(
    sapply(cmtsnum, FUN = lhs, cmts = cmts), 
    sapply(cmtsnum, FUN = rhs, cmts = cmts, k = k), 
    sep = " = "
  )
  
  if(.cat) {
    cat(ans, sep = "\n")
  } else {
    return(ans)
  }
}


# Example for an easy bicomp model
linearmodel(cmt = c("CENTRAL", "PERIPH"), 
            k = c("K10", "K12", "K21"), .cat = TRUE)
#> dxdt_CENTRAL = - (K10 + K12) * CENTRAL + K21 * PERIPH ;
#> dxdt_PERIPH = - K21 * PERIPH + K12 * CENTRAL ;

# More complex, two DEPOT compartment...
linearmodel(cmt = c("DEPOT1", "CENTRAL", "DEPOT3", "PERIPH"), 
            k = c("K12", "K32", "K24", "K42", "K20"), .cat = TRUE)
#> dxdt_DEPOT1 = - K12 * DEPOT1  ;
#> dxdt_CENTRAL = - (K24 + K20) * CENTRAL + K12 * DEPOT1 + K32 * DEPOT3 + K42 * PERIPH ;
#> dxdt_DEPOT3 = - K32 * DEPOT3  ;
#> dxdt_PERIPH = - K42 * PERIPH + K24 * CENTRAL ;

# If a compartment has no associated rate constant, then dxdt = 0
linearmodel(cmt = c("DEPOT1", "CENTRAL", "DEPOT3", "PERIPH"), 
            k = c("K12", "K24", "K42", "K20"), .cat = TRUE)
#> dxdt_DEPOT1 = - K12 * DEPOT1  ;
#> dxdt_CENTRAL = - (K24 + K20) * CENTRAL + K12 * DEPOT1 + K42 * PERIPH ;
#> dxdt_DEPOT3 = 0 ;
#> dxdt_PERIPH = - K42 * PERIPH + K24 * CENTRAL ;

# The KiTj syntax works, we can have a lot of compartments...
linearmodel(cmt = paste0("CMT", 1:10), k = c(paste0("K", 1:9, "T", 2:10), "K10T0"), .cat = TRUE)
#> dxdt_CMT1 = - K1T2 * CMT1  ;
#> dxdt_CMT2 = - K2T3 * CMT2 + K1T2 * CMT1 ;
#> dxdt_CMT3 = - K3T4 * CMT3 + K2T3 * CMT2 ;
#> dxdt_CMT4 = - K4T5 * CMT4 + K3T4 * CMT3 ;
#> dxdt_CMT5 = - K5T6 * CMT5 + K4T5 * CMT4 ;
#> dxdt_CMT6 = - K6T7 * CMT6 + K5T6 * CMT5 ;
#> dxdt_CMT7 = - K7T8 * CMT7 + K6T7 * CMT6 ;
#> dxdt_CMT8 = - K8T9 * CMT8 + K7T8 * CMT7 ;
#> dxdt_CMT9 = - K9T10 * CMT9 + K8T9 * CMT8 ;
#> dxdt_CMT10 = - K10T0 * CMT10 + K9T10 * CMT9 ;

^{Created on 2022-06-15 by the reprex package (v2.0.1)}