Going Beyond the Built-in Model Libraries | MonolixSuite Documentation

Introduction

The mlxModelFinder v2.0 package allows you to use your own custom model library instead of the built-in PK library. This is useful when:

You want to search through different types of PD models
You need structural features not available in the built-in library (e.g. saturable absorption, dose-dependent bioavailability)
You want to search over PD or PKPD model structures while keeping the PK part fixed
You have a single validated model and only want to optimize the IIV configuration
You want to leverage any Monolix library (PD, PKPD, TMDD, TGI, ...) through getLibraryModelName()

This vignette demonstrates three real-world use cases, from simple to advanced.

How It Works

To use custom models, you provide:

library = "custom" in findModel()
A model_creation_func(library, filters) that returns a model path (or a lib: reference)
A param_mapping structure if using IIV exploration (settings$findIIV = TRUE)

The model_creation_func receives the current filter combination and must return either:

An absolute path to a .txt Monolix model file, or
A lib: reference from getLibraryModelName()

Example 1: PKPD Search Using the Monolix Library

A common scenario: you have a PK model and want to find the best PD structure. Instead of writing model files manually, you can call getLibraryModelName("pkpd", ...) inside your model_creation_func to pull models directly from the Monolix library.

Running the search for PKPD models sequentially—focusing on the PK model first and then the PD model—is more efficient than conducting searches for both PK and PD simultaneously.

Model Creation Function

The function fixes the PK filters and combines them with the PD filters provided by the search algorithm:

model_func <- function(library, filters) {

  fixed_filters <- list(
    administration = "oral",
    delay = "lagTime",
    absorption = "firstOrder",
    distribution = "1compartment",
    elimination = "linear",
    parametrization = "clearance"
  )

  all_filters <- c(fixed_filters, filters)

  # Inhibition models require specifying full vs partial
  if (filters$drugAction %in% c("productionInhibition", "degradationInhibition")) {
    all_filters$inhibition <- "fullInhibition"
  }

  model_name <- lixoftConnectors::getLibraryModelName("pkpd", all_filters)

  if (length(model_name) != 1) {
    stop(
      "Filters did not resolve to a unique model. Got:\n",
      paste("  -", model_name, collapse = "\n"),
      call. = FALSE
    )
  }

  return(model_name)
}

Running the Search

library(mlxModelFinder)

result <- findModel(
  project = "my_pkpd_project.mlxtran",
  library = "custom",
  filters = list(
    response = "turnover",
    drugAction = c("productionStimulation", "productionInhibition",
                   "degradationStimulation", "degradationInhibition"),
    sigmoidicity = c("true", "false")
  ),
  model_creation_func = model_func,
  algorithm = "exhaustive_search",
  param_mapping = list(
    common_params = c("Tlag", "ka", "V", "Cl", "R0", "kout"),
    filter_params = list(
      drugAction = list(
        productionStimulation  = c("Emax", "EC50"),
        productionInhibition   = c("Imax", "IC50"),
        degradationStimulation = c("Emax", "EC50"),
        degradationInhibition  = c("Imax", "IC50")
      ),
      sigmoidicity = list(
        "true"  = c("gamma"),
        "false" = character(0)
      )
    )
  ),
  settings = list(iiv = TRUE, error = TRUE, plot = FALSE)
)

print(result$results)

This searches 4 drug actions x 2 sigmoidicity options = 8 structural models, each with IIV and error model variations. The lib: prefix returned by getLibraryModelName() is passed directly to Monolix -- no file handling required.

Tip: This pattern works with any Monolix library. Replace "pkpd" with "pd", "tmdd", "tgi", etc. and adjust the fixed/variable filters accordingly.

Example 2: Extended PK Library with Custom Features

The built-in PK library covers standard models. But what if you need saturable absorption or dose-dependent bioavailability? You can write a model_creation_func that generates Monolix model files on the fly.

Model Creation Function

This function builds model code programmatically based on filters:

custom_model_loader <- function(library, filters) {
  stopifnot(!is.null(filters$absorption), !is.null(filters$tlag),
            !is.null(filters$distribution), !is.null(filters$elimination),
            !is.null(filters$bioavailability))

  absorption      <- filters$absorption[[1]]
  tlag_opt        <- filters$tlag[[1]]
  distribution    <- filters$distribution[[1]]
  elimination     <- filters$elimination[[1]]
  bioavailability <- filters$bioavailability[[1]]

  # --- Build parameter list ---
  input_params <- c("CL", "V1", "F")
  if (absorption == "first_order")  input_params <- c(input_params, "KA")
  if (absorption == "zero_order")   input_params <- c(input_params, "D1")
  if (absorption == "saturable")    input_params <- c(input_params, "VMAX_ABS", "KM_ABS")
  if (tlag_opt == "tlag")           input_params <- c(input_params, "TLAG")
  if (distribution == "2_compartment") input_params <- c(input_params, "Q2", "V2")
  if (distribution == "3_compartment") input_params <- c(input_params, "Q2", "V2", "Q3", "V3")
  if (bioavailability == "dose_dependent") input_params <- c(input_params, "D50")

  input_line <- paste0("input = {", paste(input_params, collapse = ", "), "}")

  # --- PK block ---
  pk_lines <- "PK:"
  if (bioavailability == "constant") {
    pk_lines <- c(pk_lines, "Fbio = F")
  } else {
    pk_lines <- c(pk_lines, "Fbio = F * D50/(amtDose + D50)")
  }

  if (tlag_opt == "tlag") {
    pk_lines <- c(pk_lines, "TlagPO = TLAG")
  } else {
    pk_lines <- c(pk_lines, "TlagPO = 0")
  }

  pk_lines <- c(pk_lines,
                "depot(target=Ad, adm=2, p=Fbio, Tlag=TlagPO)",
                "depot(target=A1, adm=1, p=1, Tlag=0)")

  # --- EQUATION block ---
  eq_lines <- c("EQUATION:", "Cc = A1/V1", "k10 = CL/V1")

  if (absorption == "first_order") {
    eq_lines <- c(eq_lines, "rate_abs = KA*Ad")
  } else if (absorption == "zero_order") {
    eq_lines <- c(eq_lines,
                  "if Ad > 0",
                  "    rate_abs = amtDose/D1",
                  "else",
                  "    rate_abs = 0",
                  "end")
  } else if (absorption == "saturable") {
    eq_lines <- c(eq_lines, "rate_abs = VMAX_ABS*Ad/(KM_ABS + Ad)")
  }

  eq_lines <- c(eq_lines, "rate_elim = k10*A1")

  if (distribution == "1_compartment") {
    eq_lines <- c(eq_lines,
                  "ddt_Ad = -rate_abs",
                  "ddt_A1 = rate_abs - rate_elim")
  } else if (distribution == "2_compartment") {
    eq_lines <- c(eq_lines,
                  "k12 = Q2/V1", "k21 = Q2/V2",
                  "ddt_Ad = -rate_abs",
                  "ddt_A1 = rate_abs - rate_elim - k12*A1 + k21*A2",
                  "ddt_A2 = k12*A1 - k21*A2")
  } else if (distribution == "3_compartment") {
    eq_lines <- c(eq_lines,
                  "k12 = Q2/V1", "k21 = Q2/V2",
                  "k13 = Q3/V1", "k31 = Q3/V3",
                  "ddt_Ad = -rate_abs",
                  "ddt_A1 = rate_abs - rate_elim - k12*A1 + k21*A2 - k13*A1 + k31*A3",
                  "ddt_A2 = k12*A1 - k21*A2",
                  "ddt_A3 = k13*A1 - k31*A3")
  }

  out_lines <- c("OUTPUT:", "output = {Cc}")

  # --- Write to temp file ---
  code <- paste(c("[LONGITUDINAL]", input_line, "", pk_lines, "",
                  eq_lines, "", out_lines), collapse = "\n")

  model_path <- tempfile(pattern = "custom_pk_", fileext = ".txt")
  writeLines(code, con = model_path)
  return(model_path)
}

Parameter Mapping

The mapping must reflect which parameters each filter value introduces:

pk_param_mapping <- list(
  common_params = c("CL", "V1", "F"),
  filter_params = list(
    absorption = list(
      first_order = c("KA"),
      zero_order  = c("D1"),
      saturable   = c("VMAX_ABS", "KM_ABS")
    ),
    tlag = list(
      none = character(0),
      tlag = c("TLAG")
    ),
    distribution = list(
      "1_compartment" = character(0),
      "2_compartment" = c("Q2", "V2"),
      "3_compartment" = c("Q2", "V2", "Q3", "V3")
    ),
    elimination = list(
      linear = character(0)
    ),
    bioavailability = list(
      constant       = character(0),
      dose_dependent = c("D50")
    )
  )
)

Running the Search

result <- findModel(
  project = "my_project.mlxtran",
  library = "custom",
  filters = list(
    absorption      = c("first_order", "zero_order", "saturable"),
    tlag            = c("none", "tlag"),
    distribution    = c("1_compartment", "2_compartment", "3_compartment"),
    elimination     = "linear",
    bioavailability = c("constant", "dose_dependent")
  ),
  model_creation_func = custom_model_loader,
  param_mapping = pk_param_mapping,
  algorithm = "ACO",
  settings = list(
    iiv = TRUE,
    error = TRUE
  )
)

This creates a search space of 3 x 2 x 3 x 1 x 2 = 36 structural models, each with IIV and error model exploration.

Example 3: Single Model, IIV-Only Optimization

Sometimes you already know the structural model and just want to find the optimal IIV configuration. This is common with complex models where the structure is established but the random effects need tuning.

Consider a two-compartment PKPD model with target-mediated drug disposition (TMDD). The structural model is fixed, but with 10 parameters, there are 2^10 = 1024 possible IIV configurations to explore.

Model Creation Function

When there's only one structural model, the function simply returns it regardless of filters:

tmdd_model_content <- '[LONGITUDINAL]
input = {CL, V1, Q, V2, kon, koff, kint, ksyn, R0, F1}

PK:
depot(target=Ad, adm=1, p=F1)
compartment(cmt=1, volume=V1, concentration=Cc)
peripheral(k12=Q/V1, k21=Q/V2)

EQUATION:
kdeg = ksyn / R0
ddt_Ad  = -kon * Ad
ddt_L   = Ad*kon/V1 - CL/V1*L - kon*L*R + koff*P
ddt_R   = ksyn - kdeg*R - kon*L*R + koff*P
ddt_P   = kon*L*R - (koff + kint)*P
ddt_A2  = Q/V1*L - Q/V2*A2

Cc = L

OUTPUT:
output = {Cc}
'

tmdd_model_path <- tempfile(pattern = "tmdd_", fileext = ".txt")
writeLines(tmdd_model_content, tmdd_model_path)

# The function always returns the same model
tmdd_model_func <- function(library, filters) {
  return(tmdd_model_path)
}

Parameter Mapping and Fixed IIV

With a single structural model, there are no filter-specific parameters -- everything goes in common_params. Use fixed_iiv to lock parameters that should always (or never) have IIV:

tmdd_param_mapping <- list(
  common_params = c("CL", "V1", "Q", "V2", "kon", "koff", "kint", "ksyn", "R0", "F1"),
  filter_params = list()
)

Running the Search

Use a dummy filter with a single value, since there is only one structural model:

result <- findModel(
  project = "my_tmdd_project.mlxtran",
  library = "custom",
  filters = list(model = "tmdd_qss"),
  model_creation_func = tmdd_model_func,
  param_mapping = tmdd_param_mapping,
  algorithm = "ACO",
  fixed_iiv = list(
    CL = TRUE,   # Always keep IIV on CL
    V1 = TRUE,   # Always keep IIV on V1
    F1 = TRUE    # Always keep IIV on F1
  ),
  settings = list(
    iiv = TRUE,
    error = FALSE,
    initial_iiv_forced_iter = 0
  )
)

With 3 parameters fixed and 7 free, the search space is reduced from 1024 to 2^7 = 128 IIV combinations. The ACO algorithm explores these efficiently without testing all of them.

Tip: You can also use algorithm = "exhaustive_search" if the search space is small enough. With 128 combinations and parallel execution, exhaustive search could be feasible.

Example 4: Searching Over Tumor Growth Inhibition (TGI) Models

Tumor growth inhibition modeling is central to oncology drug development. TGI models describe how tumors grow in the absence of treatment and how drug exposure suppresses that growth. The challenge is that multiple modeling choices interact: the growth law (exponential, logistic, Gompertz, ...), the killing mechanism (log-kill vs. Norton-Simon), the drug effect dynamics (first-order, Michaelis-Menten, ...), and optional features like resistance or immune dynamics.

Manually testing all plausible combinations is impractical -- a typical search space of 4 growth models, 2 killing hypotheses, 2 dynamics, and 3 resistance options already yields 48 structural models before considering IIV or error model variations.

The Monolix TGI library provides a comprehensive collection of these models, accessible through getLibraryModelName("tgi", filters). By wrapping this in a model_creation_func, you can let findModel() systematically explore the TGI model space while keeping certain design choices fixed (e.g., treatment source, initial tumor size handling) and searching over others (e.g., growth law, killing dynamics, resistance).

# Factory function for TGI models using Monolix library
create_tgi_model_func <- function(fixed_filters = list()) {
  function(library, filters) {
    all_filters <- c(fixed_filters, filters)
    model_name <- lixoftConnectors::getLibraryModelName("tgi", all_filters)
    
    if (length(model_name) != 1) {
      stop(
        "Filters did not resolve to a unique TGI model. Got:\n",
        paste("  -", model_name, collapse = "\n"),
        call. = FALSE
      )
    }
    
    return(model_name)
  }
}

# --- Example: compare growth models and killing dynamics ---
# Fix the treatment source, resistance, and delay; search over growth + kill

tgi_model_func <- create_tgi_model_func(
  fixed_filters = list(
    initialTumorSize = "asParameter",
    kinetics = "false",
    treatment = "pkModel",
    resistance = "none",
    delay = "none",
    additionalFeature = "none",
    additionalTreatmentEffect = "none"
  )
)

result <- findModel(
  project = file.path(getDemoPath(), "8.case_studies", "tgi_project.mlxtran"),
  library = "custom",
  filters = list(
    model = c("exponential", "linear", "quadratic", "generalizedExponential"),
    killingHypothesis = c("logKill", "NortonSimon"),
    dynamics = c("firstOrder", "MichaelisMenten")
  ),
  model_creation_func = tgi_model_func,
  algorithm = "ACO",
  param_mapping = list(
    common_params = c("TS0", "V", "k"),
    filter_params = list(
      model = list(
        exponential      = c("kge"),
        linear           = c("kgl"),
        quadratic        = c("kgl", "kg2"),
        generalizedExponential = c("kge", "gamma")
      ),
      killingHypothesis = list(
        logKill     = character(0),
        NortonSimon = character(0)
      ),
      dynamics = list(
        firstOrder       = c("kkill"),
        MichaelisMenten  = c("Kmax", "EC50")
      )
    )
  ),
  settings = list(iiv = TRUE, error = TRUE)
)

Search space: 4 models x 2 killing x 2 dynamics = 16 structural models with IIV and error model variations on top.

The fixed_filters keep the wrapper flexible -- you can change which filters are fixed vs. searched without rewriting the function. For example, to also explore resistance mechanisms:

# Move resistance from fixed into the search space
tgi_model_func_v2 <- create_tgi_model_func(
  fixed_filters = list(
    initialTumorSize = "asParameter",
    kinetics = "false",
    treatment = "pkModel",
    delay = "none",
    additionalFeature = "none",
    additionalTreatmentEffect = "none"
  )
)

result <- findModel(
  project = file.path(getDemoPath(), "8.case_studies", "tgi_project.mlxtran"),
  library = "custom",
  filters = list(
    model = c("exponential", "linear", "quadratic", "generalizedExponential"),
    killingHypothesis = c("logKill", "NortonSimon"),
    dynamics = c("firstOrder", "MichaelisMenten"),
    resistance = c("none", "ClaretExponential", "resistantCells")
  ),
  model_creation_func = tgi_model_func_v2,
  algorithm = "ACO",
  param_mapping = list(
    common_params = c("TS0", "V", "k"),
    filter_params = list(
      model = list(
        exponential      = c("kge"),
        linear           = c("kgl"),
        quadratic        = c("kgl", "kg2"),
        generalizedExponential = c("kge", "gamma")
      ),
      killingHypothesis = list(
        logKill     = character(0),
        NortonSimon = character(0)
      ),
      dynamics = list(
        firstOrder       = c("kkill"),
        MichaelisMenten  = c("Kmax", "EC50")
      ),
      resistance = list(
        none               = character(0),
        ClaretExponential  = c("lambda"),
        resistantCells     = c("f", "kkillr")
      )
    )
  ),
  settings = list(iiv = TRUE, error = TRUE)
)

Additional Features

Custom Scenario Tasks

By default, each model run executes population parameter estimation, conditional distribution sampling (or conditional mode estimation if linearization = TRUE), standard error estimation, and log-likelihood estimation. You can override this with settings$tasks:

result <- findModel(
  project = "my_project.mlxtran",
  library = "custom",
  filters = list(model = "my_model"),
  model_creation_func = my_func,
  param_mapping = my_mapping,
  settings = list(
    iiv = TRUE,
    tasks = c(
      populationParameterEstimation    = TRUE,
      conditionalDistributionSampling  = TRUE,
      conditionalModeEstimation        = FALSE,
      standardErrorEstimation          = FALSE,
      logLikelihoodEstimation          = TRUE
    )
  )
)

You might want to do this if you specify a custom cost function that uses some of the results that are not calculated by default.

Custom Metric Functions

When using custom tasks, make sure your metric_func is compatible. For example, if you disable log-likelihood estimation, you cannot use the default BICc metric:

# Use AIC from importance sampling
my_metric <- function() {
  ll <- lixoftConnectors::getEstimatedLogLikelihood()
  return(ll$importanceSampling[["AIC"]])
}

result <- findModel(
  project = "my_project.mlxtran",
  library = "custom",
  filters = list(model = "my_model"),
  model_creation_func = my_func,
  param_mapping = my_mapping,
  metric_func = my_metric,
  settings = list(iiv = TRUE, error = TRUE)
)

Parameter Distributions

Override the default log-normal distribution for specific parameters:

result <- findModel(
  project = "my_project.mlxtran",
  library = "custom",
  filters = list(absorption = c("first_order", "zero_order")),
  model_creation_func = my_func,
  param_mapping = my_mapping,
  param_distributions = list(
    F1 = list(distribution = "logitNormal", limits = c(0, 1)),
    CL = "logNormal",
    V1 = "logNormal"
  ),
  settings = list(iiv = TRUE)
)

Limitations

Custom models are not compatible with algorithm = "decision_tree" -- use "exhaustive_search" or "ACO"
The requiredFilters parameter is not used with custom libraries
You are responsible for ensuring model files are valid Monolix models
When using getLibraryModelName(), verify that your filter combinations are valid for the chosen Monolix library

Summary

Use Case	model_creation_func returns	Best for
Monolix library (PKPD, PD, ...)	`lib:` reference via `getLibraryModelName()`	Searching PD/PKPD structures with fixed PK
Custom model files	Path to `.txt` file	Non-standard PK features, bespoke models
Single model, IIV search	Same path always	Optimizing random effects on complex models

The key components are:

model_creation_func: Maps filter values to model paths or lib: references
param_mapping: Defines which parameters are available for IIV in each model variant
filters: Defines the search space dimensions
fixed_iiv (optional): Locks IIV on/off for specific parameters to reduce search space