Skip to contents

Custom Models in RDarwin

Custom_Models.Rmd

Custom Spaces in RDarwin

RDarwin provides the ability to define and incorporate custom models within its framework, significantly expanding the range of model structures that can be explored and evaluated. This flexibility is achieved through the use of custom spaces, which allow you to specify models using PML (Pharmacometric Modeling Language) code.

Data

The example dataset (OneCpt_1stOrderAbsorpTlagNoRanef.csv) is located in the package’s examples folder. For convenience, we will copy it to the current (temporary) working directory.

WorkingDir <- tempdir()
DataFilePath <- file.path(WorkingDir, "OneCpt_1stOrderAbsorpTlagNoRanef.csv")
file.copy(system.file(package = "Certara.RDarwin", 
                      "examples", 
                      "OneCpt_1stOrderAbsorpTlagNoRanef.csv"),
          DataFilePath, overwrite = TRUE)

Creating Custom Spaces

The create_CustomSpace() function is the cornerstone of custom model definition in RDarwin. It takes a character string containing PML code as input and generates a CustomSpace object, which represents the parsed and structured representation of the custom model.

Example: Two-Transit Compartment Model

Let’s consider an example where we want to define a custom model with two transit compartments:

TwoTransitCpt_CustomCode <- "
  deriv(Aa1 = -Ktr * Aa1)
  deriv(Aa2 = Ktr * (Aa1 - Aa2))
  deriv(A1 = Ktr * Aa2 - Cl * C)
    dosepoint(Aa1)

    ## Drug concentration at the central compartment
    C = A1 / V

    ## Residual error model
    error(CEps = 0.1)
    observe(ConcObs = C * (1 + CEps))

    ## Model parameters
    stparm(V = tvV * exp(nV))
    stparm(Cl = tvCl * exp(nCl))
    stparm(Ktr = tvKtr * exp(nKtr))

    ## Fixed effects
    fixef(tvV = c(, 5, ))
    fixef(tvCl = c(, 1, ))
    fixef(tvKtr = c(, 1, ))

    ## Random effects
    ranef(diag(nV) = c(1))
    ranef(diag(nCl) = c(1))
    ranef(diag(nKtr) = c(1))
"

# Create the custom space
CustomSpace <- 
  create_CustomSpace(CustomCode = TwoTransitCpt_CustomCode)

The create_CustomSpace() function parses the TwoTransitCpt_CustomCode and extracts information about various model components, including derivatives, dosepoints, observations, structural parameters, fixed effects, and random effects. It then constructs a CustomSpace object that encapsulates this information. Later more models could be added to the search with add_CustomSpace() function.

Appending Custom Spaces to Existing Models

You can seamlessly integrate custom spaces with existing model sets generated by RDarwin. For instance, you can append the custom two-transit compartment model to a list of built-in PK models:

## Append the custom models to the built-in model list
models <- create_ModelPK(Absorption = c("First-Order", "Gamma"))

models <- 
  add_CustomSpace(models, TwoTransitCpt_CustomCode)

The add_CustomSpace() function adds the new custom space to the existing list of models, expanding the range of model structures available for exploration.

Modifying Structural Parameters in Custom Spaces

Once a custom space is created, you can further refine its structural parameters using the modify_StParmCustom() function. This function allows you to change the distribution type of a structural parameter, add or modify its associated fixed and random effects (Theta and Omega), and incorporate covariates. Note that this function could fail to remove ranef from the PML code. It could happen when more than one omega is presented in ranef statement. In such case a warning is given. An appropriate way to cope it is to create the PML model without corresponding structural parameter and use add_StParm() function to add it with given arguments.

In the example below, we modify the ‘Ktr’ structural parameter in our custom model.

models <-
  modify_StParmCustom(
    models,
    StParmName = "Ktr",
    Type = "LogNormal2",
    OmegaStParm = Omega("nKtr",
                        State = "Searched")
  )

The modification updates the Ktr parameter’s distribution to LogNormal2 and sets the random effect nKtr to be actively explored (i.e., “Searched”) during the model fitting.

By removing the parameter from the PML code and representing it as a separate structure, we gain greater control over its definition and modifications. Now we can modify it and all structures inside as a usual structural parameter, not custom.

# Ktr is not custom anymore:
list_StParms(models, IncludeCustom = FALSE)
#> [1] "Cl"                 "V"                  "Ka"                
#> [4] "MeanDelayTime"      "ShapeParamMinusOne" "Ktr"

# Modifying theta inside Ktr structural parameter:
models <-
  modify_Theta(
    models,
    Name = "tvKtr",
    InitialEstimates = InitialEstimate(Initial = 0.9),
    Frozen = FALSE
  )

Inspecting and Utilizing Custom Spaces

RDarwin provides functions to inspect the components of the spaces. For example, list_Dosepoints() can be used to examine the dosepoints, list_Observations can be used to examine the observations included in the models. get_ModelTermsToMap() function allows to get the model terms readily for mapping.

# IncludeCustom == TRUE will also include the dosepoints from the PML code (custom models)
list_Dosepoints(models,
                IncludeCustom = TRUE)
#> [1] "Aa"  "A1"  "Aa1"
list_Observations(models,
                  IncludeCustom = TRUE)
#> [1] "CObs"    "ConcObs"

## get the model terms to be mapped:
get_ModelTermsToMap(models)
#> For the mapping purposes, dosing Terms Aa, A1, Aa1 could be substituted to `AMT`.
#> $Required
#> [1] "id"      "time"    "Aa"      "A1"      "Aa1"     "CObs"    "ConcObs"
#> 
#> $Optional
#> [1] "Aa_Duration"  "Aa_Rate"      "A1_Duration"  "A1_Rate"      "Aa1_Duration"
#> [6] "Aa1_Rate"

Generating Template and Token Files

The write_ModelTemplateTokens() function is used to generate template and token files that represent the model space, including the custom spaces. These files then passed to pyDarwin for model selection and evaluation.

TemplateFilePath <- file.path(WorkingDir, "template.txt")
TokensFilePath <- file.path(WorkingDir, "tokens.json")

# please refer to ?specify_EngineParams for more engine parameters arguments available
generatedOutput <-
  write_ModelTemplateTokens(
    TemplateFilePath = TemplateFilePath,
    TokensFilePath = TokensFilePath,
    Description = "SearchAbsorptionType",
    Author = "Certara",
    DataFilePath = DataFilePath,
    DataMapping = c(
      id = "ID",
      time = "Time",
      AMT = "Dose",
      CObs = "CObs",
      ConcObs = "CObs"
    ),
    PMLParametersSets = models,
    EstArgs = specify_EngineParams(
      ODE = "DVERK",
      method = "QRPEM",
      stdErr = "Fisher-Score"
    )
  )
#> For a space named l577 the main dosepoint is set to Aa1
#> information stored in C:\Users\jcraig\AppData\Local\Temp\RtmpmuGajj/template.txt and C:\Users\jcraig\AppData\Local\Temp\RtmpmuGajj/tokens.json

# check tokens file
cat("tokens.json:", readLines(TokensFilePath), sep = "\n")
#> tokens.json:
#> {
#>     "PML": {
#>         "PK1FOC": [
#>             "  Aa = Dose CObs = CObs id = ID time = Time",
#>             "test() {\n\tcfMicro(A1, Cl / V, first = (Aa = Ka))\n\tC = A1 / V\n\tdosepoint(Aa, idosevar = AaDose, infdosevar = AaInfDose, infratevar = AaInfRate)\n\terror(CEps = 0.1)\n\tobserve(CObs = C * (1 + CEps))\n\t\n\tstparm(Cl = tvCl * exp( nCl ))\n\tfixef(tvCl= c(, 1, ))\n\tranef(diag(nCl) = c(1))\n\tstparm(V = tvV * exp( nV ))\n\tfixef(tvV= c(, 1, ))\n\tranef(diag(nV) = c(1))\n\tstparm(Ka = tvKa * exp( nKa ))\n\tfixef(tvKa= c(, 1, ))\n\tranef(diag(nKa) = c(1))\n\n}"
#>         ],
#>         "PK1GC": [
#>             "  A1 = Dose CObs = CObs id = ID time = Time",
#>             "test() {\n\tdelayInfCpt(A1, MeanDelayTime, ShapeParamMinusOne, out = - Cl * C, dist = Gamma)\n\tC = A1 / V\n\tdosepoint(A1, idosevar = A1Dose, infdosevar = A1InfDose, infratevar = A1InfRate)\n\terror(CEps = 0.1)\n\tobserve(CObs = C * (1 + CEps))\n\t\n\tstparm(MeanDelayTime = tvMeanDelayTime * exp( nMeanDelayTime ))\n\tfixef(tvMeanDelayTime= c(, 1, ))\n\tranef(diag(nMeanDelayTime) = c(1))\n\tstparm(ShapeParamMinusOne = tvShapeParamMinusOne * exp( nShapeParamMinusOne ))\n\tfixef(tvShapeParamMinusOne= c(, 1, ))\n\tranef(diag(nShapeParamMinusOne) = c(1))\n\tstparm(Cl = tvCl * exp( nCl ))\n\tfixef(tvCl= c(, 1, ))\n\tranef(diag(nCl) = c(1))\n\tstparm(V = tvV * exp( nV ))\n\tfixef(tvV= c(, 1, ))\n\tranef(diag(nV) = c(1))\n\n}"
#>         ],
#>         "l577": [
#>             "  Aa1 = Dose ConcObs = CObs id = ID time = Time",
#>             "test() {\n\t\n\t  deriv(Aa1 = -Ktr * Aa1)\n\t  deriv(Aa2 = Ktr * (Aa1 - Aa2))\n\t  deriv(A1 = Ktr * Aa2 - Cl * C)\n\t\tdosepoint(Aa1)\n\t\n\t\t## Drug concentration at the central compartment\n\t\tC = A1 / V\n\t\n\t\t## Residual error model\n\t\terror(CEps = 0.1)\n\t\tobserve(ConcObs = C * (1 + CEps))\n\t\n\t\t## Model parameters\n\t\tstparm(V = tvV * exp(nV))\n\t\tstparm(Cl = tvCl * exp(nCl))\n\t\t\n\t\n\t\t## Fixed effects\n\t\tfixef(tvV = c(, 5, ))\n\t\tfixef(tvCl = c(, 1, ))\n\t\t\n\t\n\t\t## Random effects\n\t\tranef(diag(nV) = c(1))\n\t\tranef(diag(nCl) = c(1))\n\t\t\n\t\n\t\n\tstparm(Ktr = exp( tvKtr {_nKtr[2]} ))\n\tfixef(tvKtr= c(, 0.9, ))\n\t{_nKtr[1]}\n\n}"
#>         ]
#>     },
#>     "_nKtr": [
#>         [
#>             "",
#>             ""
#>         ],
#>         [
#>             "ranef(diag(nKtr) = c(1))",
#>             "+ nKtr"
#>         ]
#>     ]
#> }

# check template file
cat("template.txt:", readLines(TemplateFilePath), sep = "\n")
#> template.txt:
#> ##Description: SearchAbsorptionType
#> ##Author: Certara
#> ##DATA {data_dir}/OneCpt_1stOrderAbsorpTlagNoRanef.csv
#> ##MAP {PML[1]}
#> ##MODEL {PML[2]}
#> ##ESTARGS
#>  sort=FALSE ODE=DVERK method=QRPEM stdErr=Fisher-Score
#> ##TABLES

Next section sets up the options for running the pyDarwin model search using Certara.RsNLME. Make sure you have Certara.RsNLME installed and its environment variables properly configured.

# Option setup
# RsNLME is expected to be installed
optionSetup <- create_pyDarwinOptions(
  algorithm = "EX",
  engine_adapter = "nlme",
  nlme_dir = Sys.getenv("INSTALLDIR"),
  working_dir = "{project_dir}/Results",
  output_dir = "{project_dir}/Results/output",
  temp_dir = "{project_dir}/Results/temp",
  gcc_dir = Sys.getenv("NLMEGCCDir64")
)

# Generate option file
write_pyDarwinOptions(pyDarwinOptions = optionSetup,
                      file = file.path(WorkingDir, "options.json"))

Running the Model Search with pyDarwin

To kick off the model search process using pyDarwin, we’ll use a function run_pyDarwin(). This function needs to know where your Python interpreter is located and where to find the template, tokens, and options files we created earlier.

# Example of using run_pyDarwin on Windows machine
result <- run_pyDarwin(
  InterpreterPath = "d:/Python/Projects/venv/Scripts/python.exe",
  DirectoryPath = WorkingDir,
  TemplatePath = "template.txt",
  TokensPath = "tokens.json",
  OptionsPath = "options.json",
  Wait = TRUE
)

# the best found model file could be found in
cat(result$FinalControlFile, sep = "\n")
#> ##Description: SearchAbsorptionType
#> ##Author: Certara
#> ##DATA OneCpt_1stOrderAbsorpTlagNoRanef.csv
#> ##MAP   Aa1 = Dose ConcObs = CObs id = ID time = Time
#> ##MODEL test() {
#> 
#>   deriv(Aa1 = -Ktr * Aa1)
#>   deriv(Aa2 = Ktr * (Aa1 - Aa2))
#>   deriv(A1 = Ktr * Aa2 - Cl * C)
#>   dosepoint(Aa1)
#> 
#> ## Drug concentration at the central compartment
#>   C = A1 / V
#> 
#> ## Residual error model
#>   error(CEps = 0.1)
#>   observe(ConcObs = C * (1 + CEps))
#> 
#> ## Model parameters
#>   stparm(V = tvV * exp(nV))
#>   stparm(Cl = tvCl * exp(nCl))
#> 
#> 
#> ## Fixed effects
#>   fixef(tvV = c(, 5, ))
#>   fixef(tvCl = c(, 1, ))
#> 
#> 
#> ## Random effects
#>   ranef(diag(nV) = c(1))
#>   ranef(diag(nCl) = c(1))
#> 
#> 
#> 
#>   stparm(Ktr = exp( tvKtr  ))
#>   fixef(tvKtr= c(, 0.9, ))
#> 
#> 
#>   }
#> ##ESTARGS
#>   sort=FALSE ODE=DVERK method=QRPEM stdErr=Fisher-Score
#> ##TABLES
#> 
#> 
#> ## Phenotype: ([('PML', 2), ('_nKtr', 0)])
#> ## Genotype: [2, 0]
#> ## Num non-influential tokens: 0

Summary

Custom spaces in RDarwin provide a powerful mechanism for extending the capabilities of the framework by incorporating user-defined models. The create_CustomSpace() and add_CustomSpace() functions allow you to seamlessly integrate custom PML code into your model building workflows, enabling the exploration of a wider range of model structures and enhancing the flexibility of your pharmacodynamic analyses. By combining custom spaces with RDarwin’s existing functionalities and pyDarwin’s model selection capabilities, you can achieve more comprehensive and tailored modeling approaches.