RxODE (version 0.7.2-1)

RxODE: Create an ODE-based model specification


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


RxODE(model, modName = basename(wd), wd = ifelse(RxODE.cache.directory ==
  ".", getwd(), RxODE.cache.directory), filename = NULL, do.compile = NULL,
  extraC = NULL, debug = FALSE, calcJac = NULL, calcSens = NULL,
  collapseModel = FALSE, ...)



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.


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.


character string with a working directory where to create a subdirectory according to modName. When specified, a subdirectoy 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


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


logical specifying whether to carry out the parsing of the ODE into C code and its compilation. Default is TRUE


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.


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


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


boolean indicating if RxODE will calculate the sennsitivities according to the specified ODEs.


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


any other arguments are passed to the function readLines, (e.g., encoding).

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 (closure) 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.


An object (closure) of class “RxODE” (see Chambers and Temple Lang (2001)) consisting of the following list of strings and functions:


the name of the model (a copy of the input argument).


a character string holding the source model specification.


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).


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 stifff or not.

For stiff ODE sytems (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 “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 assigments in the RxODE model code).


a function that (naively) checks for model validity, namely that the C object code reflects the latest model specification.


a string with the version of the RxODE object (not the package).


a function with one force = FALSE argument that dynamically loads the object code if needed.


a function with no argument that unloads the model object code.


removes all created model files, including C and DDL files. The model object is no longer valid and should be removed, e.g., rm(m1).


deprecated, use solve.








internal (not user callable) function.

RxODE Syntax

An RxODE model specification consists of one or more statements terminated by semi-colons, ‘;’, and optional comments (comments are delimited by # and an end-of-line marker). NB: Comments are not allowed inside statements.

A block of statements is a set of statements delimeted by curly braces, ‘{ ... }’. Statements can be either assignments or conditional if statements. Assignment statements can be: (1) “simple” assignmets, where the left hand is an identifier (i.e., variable), (2) special “time-derivative” assignments, where the left hand specifies the change of that variable with respect to time e.g., d/dt(depot), or (3) special “jacobian” assignments, where the left hand specifies the change of of the ODE with respect to one of the parameters, e.g. df(depot)/dy(kel). The “jacobian” assignments are not required, and are only useful for very stiff differential systems.

Expressions in assignment and ‘if’ statements can be numeric or logical (no character 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). (Notice that the modulo operator ‘%’ is currently not supported.)

Identifiers in an RxODE model specification can refer to:

  • state variables in the dynamic system (e.g., compartments in a pharmacokinetics/pharamcodynamics model);

  • implied input variable, t (time), podo (oral dose, for absorption models), and tlast (last time point);

  • model parameters, (ka rate of absorption, CL clearance, etc.);

  • pi, for the constant pi.

  • others, as created by assignments as part of the model specification.

Identifiers consists of case-sensitive alphanumeric characters, plus the underscore ‘_’ character. NB: the dot ‘.’ character is not a valid character identifier.

The values of these variables at pre-specified time points are saved as part of the fitted/integrated/solved model (see eventTable, in particular its member function add.sampling that defines a set of time points at which to capture a snapshot of the syste via the values of these variables).

The ODE specification mini-language is parsed with the help of the open source tool DParser, Plevyak (2015).


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).

Plevyek, J. Dparser, http://dparser.sourceforge.net. Web. 12 Oct. 2015.

See Also



Run this code
# 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, modName = "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(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)


# 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);

# }

Run the code above in your browser using DataLab