Pharmacodynamic Models in RDarwin

In addition to supporting PK models, RDarwin also offers tools for defining and customizing pharmacodynamic (PD) models. The create_ModelEmax() function is designed to generate sets of pharmacodynamic parameters for model selection, enabling you to define various options for Emax models.

Emax models arguments

The create_ModelEmax() function is particularly useful for Emax models, supporting various configurations.

Note: create_ModelEmax() currently serves as an alias for create_ModelPD()

Key arguments include:

Baseline: Adds a baseline response to the Emax model.
Fractional: Applies to Emax models with a baseline response, allowing for fractional effects.
Inhibitory: When enabled, the model parameters change from EC50/Emax to IC50/Imax for inhibitory effects.
Sigmoid: Adds a Hill coefficient to model sigmoidal dose-response curves.

Examples of different Emax models

library(Certara.RDarwin)

# Create a PD model set with default options
BaseEmax <- create_ModelEmax()
print(BaseEmax)
#> Emax 
#>  test() {
#>  E = Emax * C / (EC50 + C)
#>  error(EEps = 0.1)
#>  observe(EObs(C) = E * (1 + EEps))
#>  fcovariate(C)
#>  stparm(Emax = tvEmax * exp( nEmax ))
#>  fixef(tvEmax= c(, 1, ))
#>  ranef(diag(nEmax) = c(1))
#>  stparm(EC50 = tvEC50 * exp( nEC50 ))
#>  fixef(tvEC50= c(, 1, ))
#>  ranef(diag(nEC50) = c(1))
#> 
#> }
# Generate a model set with baseline, inhibitory, and sigmoid options
ImaxWithBaseline <- 
  create_ModelEmax(Baseline = TRUE, Inhibitory = TRUE, Sigmoid = TRUE)
print(ImaxWithBaseline)
#> ImaxE0Gam 
#>  test() {
#>  E = E0 - Imax * C^Gam / (IC50^Gam + C^Gam)
#>  error(EEps = 0.1)
#>  observe(EObs(C) = E * (1 + EEps))
#>  fcovariate(C)
#>  stparm(E0 = tvE0 * exp( nE0 ))
#>  fixef(tvE0= c(, 1, ))
#>  ranef(diag(nE0) = c(1))
#>  stparm(Imax = tvImax * exp( nImax ))
#>  fixef(tvImax= c(, 1, ))
#>  ranef(diag(nImax) = c(1))
#>  stparm(Gam = tvGam * exp( nGam ))
#>  fixef(tvGam= c(, 1, ))
#>  ranef(diag(nGam) = c(1))
#>  stparm(IC50 = tvIC50 * exp( nIC50 ))
#>  fixef(tvIC50= c(, 1, ))
#>  ranef(diag(nIC50) = c(1))
#> 
#> }

This function produces a structured list of PD models that can be customized using Certara.RDarwin tools to modify structural parameters, add covariates, and change the observation error type.

Exploring Model Variations with `ByVector` argument

The create_ModelEmax() function also provides a way to explore different Emax model configurations through its ByVector argument.

# ByVector==TRUE will use each argument value sequentially, 
# recycling the arguments with a length==1
Emax_Models_Vec <- create_ModelEmax(
  Baseline = c(FALSE, TRUE),
  Fractional = c(FALSE, TRUE),
  Inhibitory = c(FALSE, TRUE),
  Sigmoid = c(FALSE),
  ByVector = TRUE
)

# it is expected to have 2 models, one simple Emax and another one 
# Imax with a baseline and a Hill coefficient:
length(Emax_Models_Vec)
#> [1] 2

# ByVector==FALSE will expand each argument value over other arguments
Emax_Models_Expanded <- create_ModelEmax(
  Baseline = c(FALSE, TRUE),
  Fractional = c(FALSE, TRUE),
  Inhibitory = c(FALSE, TRUE),
  Sigmoid = c(FALSE, TRUE),
  ByVector = FALSE
)
#> Warning in create_ModelEmax(Baseline = c(FALSE, TRUE), Fractional = c(FALSE, :
#> A model with 'Baseline == FALSE' and 'Fractional == TRUE' is removed from the
#> search since it duplicates the model with 'Baseline == FALSE' and 'Fractional
#> == FALSE'

# note that 4 models are removed from the set since they provide the same models
# 2*2*2*2 - 4 = 12
length(Emax_Models_Expanded)
#> [1] 12

Generating the space of models

This section demonstrates how to generate a space of candidate PD models using the create_ModelEmax() function and further customize it using other RDarwin functions.

WorkingDir <- tempdir()
DataFilePath <- file.path(WorkingDir, "emax.csv")
# This data file is made available in the package examples folder, 
# however, we will copy it to our current working directory for ease of access
file.copy(system.file(package = "Certara.RDarwin", 
                      "examples", 
                      "emax.csv"),
          DataFilePath, 
          overwrite = TRUE)
TemplateFilePath <- 
  file.path(WorkingDir, "template.txt")
TokensFilePath <- 
  file.path(WorkingDir, "tokens.json")

## Create models space
## A warning is expected about models duplication
model <-
  create_ModelEmax(
    Baseline = c(TRUE, FALSE),
    Fractional = c(TRUE, FALSE),
    Inhibitory = TRUE,
    Sigmoid = c(TRUE, FALSE)
  )
#> Warning in create_ModelEmax(Baseline = c(TRUE, FALSE), Fractional = c(TRUE, : A
#> model with 'Baseline == FALSE' and 'Fractional == TRUE' is removed from the
#> search since it duplicates the model with 'Baseline == FALSE' and 'Fractional
#> == FALSE'

## Modify error type for particular models
model <- 
  modify_Observation(
    model,
    ObservationName = "EObs",
    SigmasChosen = Sigmas(
      Proportional = 0,
      MixRatio = list(PropPart = 0.1, AddPart = 0.2)
    ),
    PMLStructures = c("ImaxE0", "ImaxE0Gam")
  )

This generated space of models can be further refined using other RDarwin functions such as modify_Observation() to adjust the observation error model, add_Covariate() to incorporate covariates, and modify_StParm() to customize structural parameters. Here we will let it as is.

# write tokens and the template to the file
generatedOutput <-
  write_ModelTemplateTokens(
    TemplateFilePath = TemplateFilePath,
    TokensFilePath = TokensFilePath,
    Description = "Emax_Sigmoid_Covariate_IIV_REsError",
    Author = "Certara",
    DataFilePath = DataFilePath,
    DataMapping =  c(id = "subject", C = "c", EObs = "e"),
    PMLParametersSets = model
  )

# check tokens file
cat("\ntokens.json:", readLines(TokensFilePath), sep = "\n")
#> 
#> tokens.json:
#> {
#>     "PML": {
#>         "Imax": [
#>             " EObs = e",
#>             "test() {\n\tE = E0 * (1 - C / (IC50 + C))\n\terror(EEps = 0.1)\n\tobserve(EObs(C) = E * (1 + EEps))\n\tfcovariate(C)\n\tstparm(E0 = tvE0 * exp( nE0 ))\n\tfixef(tvE0= c(, 1, ))\n\tranef(diag(nE0) = c(1))\n\tstparm(IC50 = tvIC50 * exp( nIC50 ))\n\tfixef(tvIC50= c(, 1, ))\n\tranef(diag(nIC50) = c(1))\n\n}"
#>         ],
#>         "ImaxGam": [
#>             " EObs = e",
#>             "test() {\n\tE = E0 * (1 - C^Gam / (IC50^Gam + C^Gam))\n\terror(EEps = 0.1)\n\tobserve(EObs(C) = E * (1 + EEps))\n\tfcovariate(C)\n\tstparm(E0 = tvE0 * exp( nE0 ))\n\tfixef(tvE0= c(, 1, ))\n\tranef(diag(nE0) = c(1))\n\tstparm(Gam = tvGam * exp( nGam ))\n\tfixef(tvGam= c(, 1, ))\n\tranef(diag(nGam) = c(1))\n\tstparm(IC50 = tvIC50 * exp( nIC50 ))\n\tfixef(tvIC50= c(, 1, ))\n\tranef(diag(nIC50) = c(1))\n\n}"
#>         ],
#>         "ImaxE0": [
#>             " EObs = e",
#>             "test() {\n\tE = E0 - Imax * C / (IC50 + C)\n\terror(EEps = 0.2)\n\tobserve(EObs(C) = E + EEps*(1 + E*EMixRatio))\n\tstparm(EMixRatio = tvEMixRatio )\n\tfixef(tvEMixRatio= c(, 0.1, ))\n\t\n\tfcovariate(C)\n\tstparm(E0 = tvE0 * exp( nE0 ))\n\tfixef(tvE0= c(, 1, ))\n\tranef(diag(nE0) = c(1))\n\tstparm(Imax = tvImax * exp( nImax ))\n\tfixef(tvImax= c(, 1, ))\n\tranef(diag(nImax) = c(1))\n\tstparm(IC50 = tvIC50 * exp( nIC50 ))\n\tfixef(tvIC50= c(, 1, ))\n\tranef(diag(nIC50) = c(1))\n\n}"
#>         ],
#>         "ImaxE0Gam": [
#>             " EObs = e",
#>             "test() {\n\tE = E0 - Imax * C^Gam / (IC50^Gam + C^Gam)\n\terror(EEps = 0.2)\n\tobserve(EObs(C) = E + EEps*(1 + E*EMixRatio))\n\tstparm(EMixRatio = tvEMixRatio )\n\tfixef(tvEMixRatio= c(, 0.1, ))\n\t\n\tfcovariate(C)\n\tstparm(E0 = tvE0 * exp( nE0 ))\n\tfixef(tvE0= c(, 1, ))\n\tranef(diag(nE0) = c(1))\n\tstparm(Imax = tvImax * exp( nImax ))\n\tfixef(tvImax= c(, 1, ))\n\tranef(diag(nImax) = c(1))\n\tstparm(Gam = tvGam * exp( nGam ))\n\tfixef(tvGam= c(, 1, ))\n\tranef(diag(nGam) = c(1))\n\tstparm(IC50 = tvIC50 * exp( nIC50 ))\n\tfixef(tvIC50= c(, 1, ))\n\tranef(diag(nIC50) = c(1))\n\n}"
#>         ],
#>         "ImaxE01+": [
#>             " EObs = e",
#>             "test() {\n\tE = E0 * (1 - Imax * C / (IC50 + C))\n\terror(EEps = 0.1)\n\tobserve(EObs(C) = E * (1 + EEps))\n\tfcovariate(C)\n\tstparm(E0 = tvE0 * exp( nE0 ))\n\tfixef(tvE0= c(, 1, ))\n\tranef(diag(nE0) = c(1))\n\tstparm(Imax = tvImax * exp( nImax ))\n\tfixef(tvImax= c(, 1, ))\n\tranef(diag(nImax) = c(1))\n\tstparm(IC50 = tvIC50 * exp( nIC50 ))\n\tfixef(tvIC50= c(, 1, ))\n\tranef(diag(nIC50) = c(1))\n\n}"
#>         ],
#>         "ImaxE01+Gam": [
#>             " EObs = e",
#>             "test() {\n\tE = E0 * (1 - Imax * C^Gam / (IC50^Gam + C^Gam))\n\terror(EEps = 0.1)\n\tobserve(EObs(C) = E * (1 + EEps))\n\tfcovariate(C)\n\tstparm(E0 = tvE0 * exp( nE0 ))\n\tfixef(tvE0= c(, 1, ))\n\tranef(diag(nE0) = c(1))\n\tstparm(Imax = tvImax * exp( nImax ))\n\tfixef(tvImax= c(, 1, ))\n\tranef(diag(nImax) = c(1))\n\tstparm(Gam = tvGam * exp( nGam ))\n\tfixef(tvGam= c(, 1, ))\n\tranef(diag(nGam) = c(1))\n\tstparm(IC50 = tvIC50 * exp( nIC50 ))\n\tfixef(tvIC50= c(, 1, ))\n\tranef(diag(nIC50) = c(1))\n\n}"
#>         ]
#>     }
#> }

# check template file
cat("\ntemplate.txt:", readLines(TemplateFilePath), sep = "\n")
#> 
#> template.txt:
#> ##Description: Emax_Sigmoid_Covariate_IIV_REsError
#> ##Author: Certara
#> ##DATA {data_dir}/emax.csv
#> ##MAP  C=c {PML[1]} id = subject
#> ##MODEL {PML[2]}
#> ##ESTARGS
#>  sort=FALSE
#> ##TABLES

# refer to pyDarwin documentation for aliases explanation
workingDir <- "{project_dir}/Results"
outputDir <- "{project_dir}/Results/output"
tempDir <- "{project_dir}/Results/temp"

# Option setup
# RsNLME is expected to be installed
optionSetup <- create_pyDarwinOptions(
  algorithm = "GA",
  engine_adapter = "nlme",
  nlme_dir = Sys.getenv("INSTALLDIR"),
  working_dir = workingDir,
  output_dir = outputDir,
  temp_dir = tempDir,
  gcc_dir = Sys.getenv("NLMEGCCDir64")
)

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

Once your model space is defined, you can utilize pyDarwin to efficiently explore the candidate models and identify the most appropriate model for your data.

# Example of using run_pyDarwin on Linux machine
result <- run_pyDarwin(
  InterpreterPath = "~/darwin/venv/bin/python",
  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: Emax_Sigmoid_Covariate_IIV_REsError
#> ##Author: Certara
#> ##DATA emax.csv
#> ##MAP  C=c  EObs = e id = subject
#> ##MODEL test() {
#>   E = E0 * (1 - Imax * C^Gam / (IC50^Gam + C^Gam))
#>   error(EEps = 0.1)
#>   observe(EObs(C) = E * (1 + EEps))
#>   fcovariate(C)
#>   stparm(E0 = tvE0 * exp( nE0 ))
#>   fixef(tvE0= c(, 1, ))
#>   ranef(diag(nE0) = c(1))
#>   stparm(Imax = tvImax * exp( nImax ))
#>   fixef(tvImax= c(, 1, ))
#>   ranef(diag(nImax) = c(1))
#>   stparm(Gam = tvGam * exp( nGam ))
#>   fixef(tvGam= c(, 1, ))
#>   ranef(diag(nGam) = c(1))
#>   stparm(IC50 = tvIC50 * exp( nIC50 ))
#>   fixef(tvIC50= c(, 1, ))
#>   ranef(diag(nIC50) = c(1))
#> 
#>   }
#> ##ESTARGS
#>   sort=FALSE
#> ##TABLES
#> 
#> 
#> ## Phenotype: ([('PML', 5)])
#> ## Genotype: [1, 1, 1]
#> ## Num non-influential tokens: 0

Summary

The create_ModelEmax() function empowers you to define and customize Emax models within RDarwin. By leveraging its arguments, you can readily generate a diverse range of Emax model options, tailored to your specific PD analysis needs.

This integration of RDarwin and pyDarwin streamlines the process of building and evaluating PD models, ultimately leading to more robust and insightful analyses of drug effects.

Pharmacodynamic Models in RDarwin