Create an ODE-based model specification

Create a dynamic ODE-based model object suitably for translation into fast C code

Usage

rxode2(
  model,
  modName = basename(wd),
  wd = getwd(),
  filename = NULL,
  extraC = NULL,
  debug = FALSE,
  calcJac = NULL,
  calcSens = NULL,
  calcSens2 = NULL,
  calcSens3 = NULL,
  collapseModel = FALSE,
  package = NULL,
  ...,
  linCmtSens = c("linCmtA", "linCmtB"),
  indLin = FALSE,
  eventSens = NULL,
  verbose = FALSE,
  fullPrint = getOption("rxode2.fullPrint", FALSE),
  envir = parent.frame()
)

RxODE(
  model,
  modName = basename(wd),
  wd = getwd(),
  filename = NULL,
  extraC = NULL,
  debug = FALSE,
  calcJac = NULL,
  calcSens = NULL,
  calcSens2 = NULL,
  calcSens3 = NULL,
  collapseModel = FALSE,
  package = NULL,
  ...,
  linCmtSens = c("linCmtA", "linCmtB"),
  indLin = FALSE,
  eventSens = NULL,
  verbose = FALSE,
  fullPrint = getOption("rxode2.fullPrint", FALSE),
  envir = parent.frame()
)

rxode(
  model,
  modName = basename(wd),
  wd = getwd(),
  filename = NULL,
  extraC = NULL,
  debug = FALSE,
  calcJac = NULL,
  calcSens = NULL,
  calcSens2 = NULL,
  calcSens3 = NULL,
  collapseModel = FALSE,
  package = NULL,
  ...,
  linCmtSens = c("linCmtA", "linCmtB"),
  indLin = FALSE,
  eventSens = NULL,
  verbose = FALSE,
  fullPrint = getOption("rxode2.fullPrint", FALSE),
  envir = parent.frame()
)

Arguments

model

This is the ODE model specification. It can be:

a string containing the set of ordinary differential equations (ODE) and other expressions defining the changes in the dynamic system.
a file name where the ODE system equation is contained

An ODE expression enclosed in \{\}

(see also the filename argument). For details, see the sections “Details” and rxode2 Syntax below.

modName

a string to be used as the model name. This string is used for naming various aspects of the computations, including generating C symbol names, dynamic libraries, etc. Therefore, it is necessary that modName consists of simple ASCII alphanumeric characters starting with a letter.

wd

character string with a working directory where to create a subdirectory according to modName. When specified, a subdirectory named after the “modName.d” will be created and populated with a C file, a dynamic loading library, plus various other working files. If missing, the files are created (and removed) in the temporary directory, and the rxode2 DLL for the model is created in the current directory named rx_????_platform, for example rx_129f8f97fb94a87ca49ca8dafe691e1e_i386.dll

filename

A file name or connection object where the ODE-based model specification resides. Only one of model or filename may be specified.

extraC

Extra c code to include in the model. This can be useful to specify functions in the model. These C functions should usually take double precision arguments, and return double precision values.

debug

is a boolean indicating if the executable should be compiled with verbose debugging information turned on.

calcJac

boolean indicating if rxode2 will calculate the Jacobain according to the specified ODEs.

calcSens

boolean indicating if rxode2 will calculate the sensitivities according to the specified ODEs. May also be a character vector of the states/parameters whose first-order sensitivities (rx__sens_<state>_BY_<param>__) should be generated.

calcSens2

character vector (or NULL) requesting second-order sensitivities in addition to the first-order ones from calcSens. When supplied, rxode2 also generates the rx__sens_<state>_BY_<p>_BY_<q>__ compartments (the Hessian path), where p ranges over calcSens and q over calcSens2. Used, for example, for population (THETA) second-order event sensitivities. NULL (the default) skips the second-order generation.

calcSens3

character vector (or NULL) requesting third-order sensitivities in addition to the first- and second-order ones. Requires calcSens2 to also be supplied (every calcSens3 parameter needs its own already-built second-order sensitivity compartment, from the pairing of calcSens2 with calcSens3, to reference – so calcSens3 should be a subset of calcSens2, which itself should be a subset of calcSens, mirroring how calcSens2 is used everywhere else in rxode2 today). When supplied, generates the rx__sens_<state>_BY_<p>_BY_<q>_BY_<r>__ compartments, where p ranges over calcSens, q over calcSens2, and r over calcSens3, via rxExpandSens3_(). NULL (the default) skips the third-order generation.

collapseModel

boolean indicating if rxode2 will remove all LHS variables when calculating sensitivities.

package

Package name for pre-compiled binaries.

...

ignored arguments.

linCmtSens

The method to calculate the linCmt() solutions

indLin

Calculate inductive linearization matrices and compile with inductive linearization support.

eventSens

controls how dosing/event-parameter (alag, F, rate, dur, amt) sensitivities are computed when sensitivities are generated: "jump" injects the analytic event ("jump") sensitivities into the sensitivity states at each dosing event, "fd" keeps the legacy finite-difference behavior (the backward-compatible opt-out), and "both" computes both for cross-checking. NULL (the default) uses getOption("rxode2.eventSens", "fd"). When not "fd" and calcSens is supplied, calcJac is forced to TRUE so the Jacobian is available for the jump injection.

verbose

When TRUE be verbose with the linear compartmental model

fullPrint

When using printf within the model, if TRUE print on every step (except ME/indLin), otherwise when FALSE print only when calculating the d/dt

envir

is the environment to look for R user functions (defaults to parent environment)

Value

An object (environment) of class rxode2 (see Chambers and Temple Lang (2001)) consisting of the following list of strings and functions:

* `model` a character string holding the source model specification.
* `get.modelVars`a function that returns a list with 3 character
    vectors, `params`, `state`, and `lhs` of variable names used in the model
    specification. These will be output when the model is computed (i.e., the ODE solved by integration).

  * `solve`{this function solves (integrates) the ODE. This
      is done by passing the code to [rxSolve()].
      This is as if you called `rxSolve(rxode2object, ...)`,
      but returns a matrix instead of a rxSolve object.

      `params`: a numeric named vector with values for every parameter
      in the ODE system; the names must correspond to the parameter
      identifiers used in the ODE specification;

      `events`: an `eventTable` object describing the
      input (e.g., doses) to the dynamic system and observation
      sampling time points (see  [eventTable()]);

      `inits`: a vector of initial values of the state variables
      (e.g., amounts in each compartment), and the order in this vector
      must be the same as the state variables (e.g., PK/PD compartments);


      `stiff`: a logical (`TRUE` by default) indicating whether
      the ODE system is stiff or not.

      For stiff ODE systems (`stiff = TRUE`), `rxode2` uses
      the LSODA (Livermore Solver for Ordinary Differential Equations)
      Fortran package, which implements an automatic method switching
      for stiff and non-stiff problems along the integration interval,
      authored by Hindmarsh and Petzold (2003).

      For non-stiff systems (`stiff = FALSE`), `rxode2` uses `DOP853`,
      an explicit Runge-Kutta method of order 8(5, 3) of Dormand and Prince
      as implemented in C by Hairer and Wanner (1993).

      `trans_abs`: a logical (`FALSE` by default) indicating
      whether to fit a transit absorption term
      (TODO: need further documentation and example);

      `atol`: a numeric absolute tolerance (1e-08 by default);

      `rtol`: a numeric relative tolerance (1e-06 by default).

      The output of \dQuote{solve} is a matrix with as many rows as there
      are sampled time points and as many columns as system variables
      (as defined by the ODEs and additional assignments in the rxode2 model
          code).}

  * `isValid` a function that (naively) checks for model validity,
      namely that the C object code reflects the latest model
      specification.
  * `version` a string with the version of the `rxode2`
      object (not the package).
  * `dynLoad` a function with one `force = FALSE` argument
      that dynamically loads the object code if needed.
  * `dynUnload` a function with no argument that unloads
      the model object code.
  * `delete` removes all created model files, including C and DLL files.
      The model object is no longer valid and should be removed, e.g.,
      `rm(m1)`.
  * `run` deprecated, use `solve`.
  * `get.index` deprecated.
  * `getObj` internal (not user callable) function.

Details

The Rx in the name rxode2 is meant to suggest the abbreviation Rx for a medical prescription, and thus to suggest the package emphasis on pharmacometrics modeling, including pharmacokinetics (PK), pharmacodynamics (PD), disease progression, drug-disease modeling, etc.

References

Chamber, J. M. and Temple Lang, D. (2001) Object Oriented Programming in R. R News, Vol. 1, No. 3, September 2001. https://cran.r-project.org/doc/Rnews/Rnews_2001-3.pdf.

Hindmarsh, A. C. ODEPACK, A Systematized Collection of ODE Solvers. Scientific Computing, R. S. Stepleman et al. (Eds.), North-Holland, Amsterdam, 1983, pp. 55-64.

Petzold, L. R. Automatic Selection of Methods for Solving Stiff and Nonstiff Systems of Ordinary Differential Equations. Siam J. Sci. Stat. Comput. 4 (1983), pp. 136-148.

Hairer, E., Norsett, S. P., and Wanner, G. Solving ordinary differential equations I, nonstiff problems. 2nd edition, Springer Series in Computational Mathematics, Springer-Verlag (1993).

Plevyak, J. dparser, https://dparser.sourceforge.net/. Web. 12 Oct. 2015.

Author

Melissa Hallow, Wenping Wang and Matthew Fidler

Examples

# \donttest{

mod <- function() {
  ini({
    KA   <- .291
    CL   <- 18.6
    V2   <- 40.2
    Q    <- 10.5
    V3   <- 297.0
    Kin  <- 1.0
    Kout <- 1.0
    EC50 <- 200.0
  })
  model({
    # A 4-compartment model, 3 PK and a PD (effect) compartment
    # (notice state variable names 'depot', 'centr', 'peri', 'eff')
    C2 <- centr/V2
    C3 <- peri/V3
    d/dt(depot) <- -KA*depot;
    d/dt(centr) <- KA*depot - CL*C2 - Q*C2 + Q*C3;
    d/dt(peri)  <-                    Q*C2 - Q*C3;
    d/dt(eff)   <- Kin - Kout*(1-C2/(EC50+C2))*eff;
    eff(0)      <- 1
  })
}

m1 <- rxode2(mod)
#>  
#>  
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
print(m1)
#>  ── rxode2-based free-form 4-cmt ODE model ────────────────────────────────────── 
#>  ── Initalization: ──  
#> Fixed Effects ($theta): 
#>      KA      CL      V2       Q      V3     Kin    Kout    EC50 
#>   0.291  18.600  40.200  10.500 297.000   1.000   1.000 200.000 
#> 
#> States ($state or $stateDf): 
#>   Compartment Number Compartment Name
#> 1                  1            depot
#> 2                  2            centr
#> 3                  3             peri
#> 4                  4              eff
#>  ── Model (Normalized Syntax): ── 
#> function() {
#>     ini({
#>         KA <- 0.291
#>         CL <- 18.6
#>         V2 <- 40.2
#>         Q <- 10.5
#>         V3 <- 297
#>         Kin <- 1
#>         Kout <- 1
#>         EC50 <- 200
#>     })
#>     model({
#>         C2 <- centr/V2
#>         C3 <- peri/V3
#>         d/dt(depot) <- -KA * depot
#>         d/dt(centr) <- KA * depot - CL * C2 - Q * C2 + Q * C3
#>         d/dt(peri) <- Q * C2 - Q * C3
#>         d/dt(eff) <- Kin - Kout * (1 - C2/(EC50 + C2)) * eff
#>         eff(0) <- 1
#>     })
#> }

# Step 2 - Create the model input as an EventTable,
# including dosing and observation (sampling) events

# QD (once daily) dosing for 5 days.

qd <- et(amountUnits = "ug", timeUnits = "hours") |>
  et(amt = 10000, addl = 4, ii = 24)

# Sample the system hourly during the first day, every 8 hours
# then after
qd <- qd |> et(0:24) |>
  et(from = 24 + 8, to = 5 * 24, by = 8)

# Step 3 - solve the system

qd.cp <- rxSolve(m1, qd)
#>  
#>  

head(qd.cp)
#>    time       C2        C3     depot    centr      peri      eff
#> 1 0 [h]  0.00000 0.0000000 10000.000    0.000    0.0000 1.000000
#> 2 1 [h] 43.99334 0.9113641  7475.157 1768.532  270.6751 1.083968
#> 3 2 [h] 54.50866 2.6510696  5587.797 2191.248  787.3677 1.179529
#> 4 3 [h] 51.65163 4.4243597  4176.966 2076.396 1314.0348 1.227523
#> 5 4 [h] 44.37513 5.9432612  3122.347 1783.880 1765.1486 1.233503
#> 6 5 [h] 36.46382 7.1389804  2334.004 1465.845 2120.2772 1.214084

# }