RxODE: Create an ODE-based model specification

In RxODE: Facilities for Simulating from ODE-Based Models

Create an ODE-based model specification

Description

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

Usage

RxODE(
  model,
  modName = basename(wd),
  wd = getwd(),
  filename = NULL,
  extraC = NULL,
  debug = FALSE,
  calcJac = NULL,
  calcSens = NULL,
  collapseModel = FALSE,
  package = NULL,
  ...,
  linCmtSens = c("linCmtA", "linCmtB", "linCmtC"),
  indLin = FALSE,
  verbose = FALSE
)

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 RxODE 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 RxODE 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 RxODE will calculate the Jacobain according to the specified ODEs.
calcSens	boolean indicating if RxODE will calculate the sensitivities according to the specified ODEs.
collapseModel	boolean indicating if RxODE 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.
verbose	When TRUE be verbose with the linear compartmental model

Details

The Rx in the name RxODE 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.

The ODE-based model specification may be coded inside a character string or in a text file, see Section RxODE Syntax below for coding details. An internal RxODE compilation manager object translates the ODE system into C, compiles it, and dynamically loads the object code into the current R session. The call to RxODE produces an object of class RxODE which consists of a list-like structure (environment) with various member functions (see Section Value below).

For evaluating RxODE models, two types of inputs may be provided: a required set of time points for querying the state of the ODE system and an optional set of doses (input amounts). These inputs are combined into a single event table object created with the function eventTable() or et().

An RxODE model specification consists of one or more statements optionally terminated by semi-colons ; and optional comments (comments are delimited by # and an end-of-line).

A block of statements is a set of statements delimited by curly braces, { ... }.

Statements can be either assignments, conditional if/else if/else, while loops (can be exited by break), special statements, or printing statements (for debugging/testing)

Assignment statements can be:

• simple assignments, where the left hand is an identifier (i.e., variable)

• special time-derivative assignments, where the left hand specifies the change of the amount in the corresponding state variable (compartment) with respect to time e.g., d/dt(depot):

• special initial-condition assignments where the left hand specifies the compartment of the initial condition being specified, e.g. depot(0) = 0

• special model event changes including bioavailability (f(depot)=1), lag time (alag(depot)=0), modeled rate (rate(depot)=2) and modeled duration (dur(depot)=2). An example of these model features and the event specification for the modeled infusions the RxODE data specification is found in RxODE events vignette.

• special change point syntax, or model times. These model times are specified by mtime(var)=time

• special Jacobian-derivative assignments, where the left hand specifies the change in the compartment ode with respect to a variable. For example, if d/dt(y) = dy, then a Jacobian for this compartment can be specified as df(y)/dy(dy) = 1. There may be some advantage to obtaining the solution or specifying the Jacobian for very stiff ODE systems. However, for the few stiff systems we tried with LSODA, this actually slightly slowed down the solving.

Note that assignment can be done by =, <- or ~.

When assigning with the ~ operator, the simple assignments and time-derivative assignments will not be output.

Special statements can be:

• Compartment declaration statements, which can change the default dosing compartment and the assumed compartment number(s) as well as add extra compartment names at the end (useful for multiple-endpoint nlmixr models); These are specified by cmt(compartmentName)

• Parameter declaration statements, which can make sure the input parameters are in a certain order instead of ordering the parameters by the order they are parsed. This is useful for keeping the parameter order the same when using 2 different ODE models. These are specified by param(par1, par2,...)

An example model is shown below:

   # simple assignment
   C2 = centr/V2;

   # time-derivative assignment
   d/dt(centr) = F*KA*depot - CL*C2 - Q*C2 + Q*C3; 

Expressions in assignment and if statements can be numeric or logical, however, no character nor integer expressions are currently supported.

Numeric expressions can include the following numeric operators +, -, *, /, ^ and those mathematical functions defined in the C or the R math libraries (e.g., fabs, exp, log, sin, abs).

You may also access the R’s functions in the R math libraries, like lgammafn for the log gamma function.

The RxODE syntax is case-sensitive, i.e., ABC is different than abc, Abc, ABc, etc.

Identifiers

Like R, Identifiers (variable names) may consist of one or more alphanumeric, underscore _ or period . characters, but the first character cannot be a digit or underscore _.

Identifiers in a model specification can refer to:

• State variables in the dynamic system (e.g., compartments in a pharmacokinetics model).

• Implied input variable, t (time), tlast (last time point), and podo (oral dose, in the undocumented case of absorption transit models).

• Special constants like pi or R’s predefined constants.

• Model parameters (e.g., ka rate of absorption, CL clearance, etc.)

• Others, as created by assignments as part of the model specification; these are referred as LHS (left-hand side) variable.

Currently, the RxODE modeling language only recognizes system state variables and “parameters”, thus, any values that need to be passed from R to the ODE model (e.g., age) should be either passed in the params argument of the integrator function rxSolve() or be in the supplied event data-set.

There are certain variable names that are in the RxODE event tables. To avoid confusion, the following event table-related items cannot be assigned, or used as a state but can be accessed in the RxODE code:

• cmt

• dvid

• addl

• ss

• rate

• id

However the following variables are cannot be used in a model specification:

• evid

• ii

Sometimes RxODE generates variables that are fed back to RxODE. Similarly, nlmixr generates some variables that are used in nlmixr estimation and simulation. These variables start with the either the rx or nlmixr prefixes. To avoid any problems, it is suggested to not use these variables starting with either the rx or nlmixr prefixes.

Logical Operators

Logical operators support the standard R operators ==, != >= <= > and <. Like R these can be in if() or while() statements, ifelse() expressions. Additionally they can be in a standard assignment. For instance, the following is valid:

cov1 = covm*(sexf == "female") + covm*(sexf != "female")

Notice that you can also use character expressions in comparisons. This convenience comes at a cost since character comparisons are slower than numeric expressions. Unlike R, as.numeric or as.integer for these logical statements is not only not needed, but will cause an syntax error if you try to use the function.

Value

An object (environment) of class RxODE (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(RxODEobject, ...)`,
      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`), `RxODE` 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`), `RxODE` 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).e

      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 RxODE 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 `RxODE`
      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.

Author(s)

Melissa Hallow, Wenping Wang and Matthew Fidler

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, http://dparser.sourceforge.net. Web. 12 Oct. 2015.

See Also

eventTable(), et(), add.sampling(), add.dosing()

Examples

# Step 1 - Create a model specification
ode <- "
   # 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;
"

m1 <- RxODE(model = ode)
print(m1)

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

# QD (once daily) dosing for 5 days.

qd <- eventTable(amount.units = "ug", time.units = "hours")
qd$add.dosing(dose = 10000, nbr.doses = 5, dosing.interval = 24)

# Sample the system hourly during the first day, every 8 hours
# then after

qd$add.sampling(0:24)
qd$add.sampling(seq(from = 24 + 8, to = 5 * 24, by = 8))

# Step 3 - set starting parameter estimates and initial
# values of the state

theta <-
  c(
    KA = .291, CL = 18.6,
    V2 = 40.2, Q = 10.5, V3 = 297.0,
    Kin = 1.0, Kout = 1.0, EC50 = 200.0
  )

# init state variable
inits <- c(0, 0, 0, 1)
# Step 4 - Fit the model to the data

qd.cp <- m1$solve(theta, events = qd, inits)

head(qd.cp)

# This returns a matrix.  Note that you can also
# solve using name initial values. For example:

inits <- c(eff = 1)
qd.cp <- solve(m1, theta, events = qd, inits)
print(qd.cp)

plot(qd.cp)

RxODE documentation built on March 23, 2022, 9:06 a.m.