Chapter 6 rxode2 syntax

This briefly describes the syntax used to define models that rxode2 will translate into R-callable compiled code. It also describes the communication of variables between R and the rxode2 modeling specification.

6.1 Example

   # An rxode2 model specification (this line is a comment).

   if(comed==0){   # concomitant medication (con-med)?
      F = 1.0;     # full bioavailability w.o. con-med
   }
   else {
      F = 0.80;    # 20% reduced bioavailability
   }

   C2 = centr/V2;  # concentration in the central compartment
   C3 = peri/V3;   # concentration in the peripheral compartment

   # ODE describing the PK and PD

   d/dt(depot) = -KA*depot;
   d/dt(centr) = F*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;

6.2 Syntax

This basic 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 rxode2 data specification is found in rxode2 events section.

  • 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. Note that with the rxode2 model functions assignment with ~ can also be overloaded with a residual distribution specification.

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.

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 rxode2 syntax is case-sensitive, i.e., ABC is different than abc, Abc, ABc, etc.

6.2.1 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 rxode2 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 rxode2 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 rxode2 code:

  • cmt
  • dvid
  • addl
  • ss
  • rate
  • id

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

  • evid
  • ii

Sometimes rxode2 generates variables that are fed back to rxode2. Similarly, nlmixr2 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.

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

6.4 Supported functions

All the supported functions in rxode2 can be seen with the rxSupportedFuns().

A brief description of the built-in functions are in the following table:

Note that lag(cmt) = is equivalent to alag(cmt) = and not the same as = lag(wt)

6.5 Reserved keywords

There are a few reserved keywords in a rxode2 model. They are in the following table:

Note that rxFlag will always output 11 or calc_lhs since that is where the final variables are calculated, though you can tweak or test certain parts of rxode2 by using this flag.

6.6 Residual functions when using rxode2 functions

In addition to ~ hiding output for certain types of output, it also is used to specify a residual output or endpoint when the input is an rxode2 model function (that includes the residual in the model({}) block).

These specifications are of the form:

var ~ add(add.sd)

Indicating the variable var is the variable that represents the individual central tendencies of the model and it also represents the compartment specification in the data-set.

You can also change the compartment name using the | syntax, that is:

var ~ add(add.sd) | cmt

In the above case var represents the central tendency and cmt represents the compartment or dvid specification.

6.6.1 Transformations

For normal and related distributions, you can apply the transformation on both sides by using some keywords/functions to apply these transformations.

Transformation rxode2/nlmixr2 code
Box-Cox +coxBox(lambda)
Yeo-Johnson +yeoJohnson(lambda)
logit-normal +logitNorm(logit.sd, low, hi)
probit-normal +probitNorm(probid.sd, low, hi)
log-normal +lnorm(lnorm.sd)

By default for the likelihood for all of these transformations is calculated on the untransformed scale.

For bounded variables like logit-normal or probit-normal the low and high values are defaulted to 0 and 1 if missing.

For models where you wish to have a proportional model on one of these transformation you can replace the standard deviation with NA

To allow for more transformations, lnorm(), probitNorm() and logitNorm() can be combined the variance stabilizing yeoJohnson() transformation.

6.6.3 Notes on additive + proportional models

There are two different ways to specify additive and proportional models, which we will call combined1 and combined2, the same way that Monolix calls the two distributions (to avoid between software differences in naming).

The first, combined1, assumes that the additive and proportional differences are on the standard deviation scale, or:

\[ y=f+(a+b\times f^c)\times\varepsilon \]

The second, combined2, assumes that the additive and proportional differences are combined on a variance scale:

\[ y=f+\left(\sqrt{a^2+b^2\times f^{2c}}\right)\times\varepsilon \]

The default in nlmixr2/rxode2 if not otherwise specified is combined2 since it mirrors how adding 2 normal distributions in statistics will add their variances (not the standard deviations). However, the combined1 can describe the data possibly even better than combined2 so both are possible options in rxode2/nlmixr2.

6.6.4 Distributions of known likelihoods

For residuals that are not related to normal, t-distribution or cauchy, often the residual specification is of the form:

cmt ~ dbeta(alpha, beta)

Where the compartment specification is on the left handed side of the specification.

For generalized likelihood you can specify:

ll(cmt) ~ llik specification

6.6.5 Ordinal likelihoods

Finally, ordinal likelihoods/simulations can be specified in 2 ways. The first is:

err ~ c(p0, p1, p2)

Here err represents the compartment and p0 is the probability of being in a specific category:

Category Probability
1 p0
2 p1
3 p2
4 1-p0-p1-p2

It is up to the model to ensure that the sum of the p values are less than 1. Additionally you can write an arbitrary number of categories in the ordinal model described above.

It seems a little off that p0 is the probability for category 1 and sometimes scores are in non-whole numbers. This can be modeled as follows:

err ~ c(p0=0, p1=1, p2=2, 3)

Here the numeric categories are specified explicitly, and the probabilities remain the same:

Category Probability
0 p0
1 p1
2 p2
3 1-p0-p1-p2

6.6.6 General table of supported residual distributions

In general all the that are supported are in the following table (available in rxode2::rxResidualError)

6.7 cmt() changing compartment numbers for states

The compartment order can be changed with the cmt() syntax in the model. To understand what the cmt() can do you need to understand how rxode2 numbers the compartments.

Below is an example of how rxode2 numbers compartments

6.7.1 How rxode2 numbers compartments

rxode2 automatically assigns compartment numbers when parsing. For example, with the Mavoglurant PBPK model the following model may be used:

library(rxode2)

pbpk <- function() {
  model({
    KbBR = exp(lKbBR)
    KbMU = exp(lKbMU)
    KbAD = exp(lKbAD)
    CLint= exp(lCLint + eta.LClint)
    KbBO = exp(lKbBO)
    KbRB = exp(lKbRB)

    ## Regional blood flows
    # Cardiac output (L/h) from White et al (1968)
    CO  = (187.00*WT^0.81)*60/1000 
    QHT = 4.0 *CO/100
    QBR = 12.0*CO/100
    QMU = 17.0*CO/100
    QAD = 5.0 *CO/100
    QSK = 5.0 *CO/100
    QSP = 3.0 *CO/100
    QPA = 1.0 *CO/100
    QLI = 25.5*CO/100
    QST = 1.0 *CO/100
    QGU = 14.0*CO/100
    # Hepatic artery blood flow
    QHA = QLI - (QSP + QPA + QST + QGU) 
    QBO = 5.0 *CO/100
    QKI = 19.0*CO/100
    QRB = CO - (QHT + QBR + QMU + QAD + QSK + QLI + QBO + QKI)
    QLU = QHT + QBR + QMU + QAD + QSK + QLI + QBO + QKI + QRB

    ## Organs' volumes = organs' weights / organs' density
    VLU = (0.76 *WT/100)/1.051
    VHT = (0.47 *WT/100)/1.030
    VBR = (2.00 *WT/100)/1.036
    VMU = (40.00*WT/100)/1.041
    VAD = (21.42*WT/100)/0.916
    VSK = (3.71 *WT/100)/1.116
    VSP = (0.26 *WT/100)/1.054
    VPA = (0.14 *WT/100)/1.045
    VLI = (2.57 *WT/100)/1.040
    VST = (0.21 *WT/100)/1.050
    VGU = (1.44 *WT/100)/1.043
    VBO = (14.29*WT/100)/1.990
    VKI = (0.44 *WT/100)/1.050
    VAB = (2.81 *WT/100)/1.040
    VVB = (5.62 *WT/100)/1.040
    VRB = (3.86 *WT/100)/1.040

    ## Fixed parameters
    BP = 0.61      # Blood:plasma partition coefficient
    fup = 0.028    # Fraction unbound in plasma
    fub = fup/BP   # Fraction unbound in blood

    KbLU = exp(0.8334)
    KbHT = exp(1.1205)
    KbSK = exp(-.5238)
    KbSP = exp(0.3224)
    KbPA = exp(0.3224)
    KbLI = exp(1.7604)
    KbST = exp(0.3224)
    KbGU = exp(1.2026)
    KbKI = exp(1.3171)

    ##-----------------------------------------
    S15 = VVB*BP/1000
    C15 = Venous_Blood/S15

    ##-----------------------------------------
    d/dt(Lungs) = QLU*(Venous_Blood/VVB - Lungs/KbLU/VLU)
    d/dt(Heart) = QHT*(Arterial_Blood/VAB - Heart/KbHT/VHT)
    d/dt(Brain) = QBR*(Arterial_Blood/VAB - Brain/KbBR/VBR)
    d/dt(Muscles) = QMU*(Arterial_Blood/VAB - Muscles/KbMU/VMU)
    d/dt(Adipose) = QAD*(Arterial_Blood/VAB - Adipose/KbAD/VAD)
    d/dt(Skin) = QSK*(Arterial_Blood/VAB - Skin/KbSK/VSK)
    d/dt(Spleen) = QSP*(Arterial_Blood/VAB - Spleen/KbSP/VSP)
    d/dt(Pancreas) = QPA*(Arterial_Blood/VAB - Pancreas/KbPA/VPA)
    d/dt(Liver) = QHA*Arterial_Blood/VAB + QSP*Spleen/KbSP/VSP +
      QPA*Pancreas/KbPA/VPA + QST*Stomach/KbST/VST +
      QGU*Gut/KbGU/VGU - CLint*fub*Liver/KbLI/VLI - QLI*Liver/KbLI/VLI
    d/dt(Stomach) = QST*(Arterial_Blood/VAB - Stomach/KbST/VST)
    d/dt(Gut) = QGU*(Arterial_Blood/VAB - Gut/KbGU/VGU)
    d/dt(Bones) = QBO*(Arterial_Blood/VAB - Bones/KbBO/VBO)
    d/dt(Kidneys) = QKI*(Arterial_Blood/VAB - Kidneys/KbKI/VKI)
    d/dt(Arterial_Blood) = QLU*(Lungs/KbLU/VLU - Arterial_Blood/VAB)
    d/dt(Venous_Blood) = QHT*Heart/KbHT/VHT + QBR*Brain/KbBR/VBR +
      QMU*Muscles/KbMU/VMU + QAD*Adipose/KbAD/VAD + QSK*Skin/KbSK/VSK +
      QLI*Liver/KbLI/VLI + QBO*Bones/KbBO/VBO + QKI*Kidneys/KbKI/VKI +
      QRB*Rest_of_Body/KbRB/VRB - QLU*Venous_Blood/VVB
    d/dt(Rest_of_Body) = QRB*(Arterial_Blood/VAB - Rest_of_Body/KbRB/VRB)
  })
}

If you look at the printout, you can see where rxode2 assigned the compartment number(s)

pbpk <- pbpk()
print(pbpk)
#>  ── rxode2-based free-form 16-cmt ODE model ───────────────── 
#> 
#> States ($state or $stateDf): 
#>    Compartment Number Compartment Name
#> 1                   1            Lungs
#> 2                   2            Heart
#> 3                   3            Brain
#> 4                   4          Muscles
#> 5                   5          Adipose
#> 6                   6             Skin
#> 7                   7           Spleen
#> 8                   8         Pancreas
#> 9                   9            Liver
#> 10                 10          Stomach
#> 11                 11              Gut
#> 12                 12            Bones
#> 13                 13          Kidneys
#> 14                 14   Arterial_Blood
#> 15                 15     Venous_Blood
#> 16                 16     Rest_of_Body
#>  ── Model (Normalized Syntax): ── 
#> function() {
#>     model({
#>         KbBR = exp(lKbBR)
#>         KbMU = exp(lKbMU)
#>         KbAD = exp(lKbAD)
#>         CLint = exp(lCLint + eta.LClint)
#>         KbBO = exp(lKbBO)
#>         KbRB = exp(lKbRB)
#>         CO = (187 * WT^0.81) * 60/1000
#>         QHT = 4 * CO/100
#>         QBR = 12 * CO/100
#>         QMU = 17 * CO/100
#>         QAD = 5 * CO/100
#>         QSK = 5 * CO/100
#>         QSP = 3 * CO/100
#>         QPA = 1 * CO/100
#>         QLI = 25.5 * CO/100
#>         QST = 1 * CO/100
#>         QGU = 14 * CO/100
#>         QHA = QLI - (QSP + QPA + QST + QGU)
#>         QBO = 5 * CO/100
#>         QKI = 19 * CO/100
#>         QRB = CO - (QHT + QBR + QMU + QAD + QSK + QLI + QBO + 
#>             QKI)
#>         QLU = QHT + QBR + QMU + QAD + QSK + QLI + QBO + QKI + 
#>             QRB
#>         VLU = (0.76 * WT/100)/1.051
#>         VHT = (0.47 * WT/100)/1.03
#>         VBR = (2 * WT/100)/1.036
#>         VMU = (40 * WT/100)/1.041
#>         VAD = (21.42 * WT/100)/0.916
#>         VSK = (3.71 * WT/100)/1.116
#>         VSP = (0.26 * WT/100)/1.054
#>         VPA = (0.14 * WT/100)/1.045
#>         VLI = (2.57 * WT/100)/1.04
#>         VST = (0.21 * WT/100)/1.05
#>         VGU = (1.44 * WT/100)/1.043
#>         VBO = (14.29 * WT/100)/1.99
#>         VKI = (0.44 * WT/100)/1.05
#>         VAB = (2.81 * WT/100)/1.04
#>         VVB = (5.62 * WT/100)/1.04
#>         VRB = (3.86 * WT/100)/1.04
#>         BP = 0.61
#>         fup = 0.028
#>         fub = fup/BP
#>         KbLU = exp(0.8334)
#>         KbHT = exp(1.1205)
#>         KbSK = exp(-0.5238)
#>         KbSP = exp(0.3224)
#>         KbPA = exp(0.3224)
#>         KbLI = exp(1.7604)
#>         KbST = exp(0.3224)
#>         KbGU = exp(1.2026)
#>         KbKI = exp(1.3171)
#>         S15 = VVB * BP/1000
#>         C15 = Venous_Blood/S15
#>         d/dt(Lungs) = QLU * (Venous_Blood/VVB - Lungs/KbLU/VLU)
#>         d/dt(Heart) = QHT * (Arterial_Blood/VAB - Heart/KbHT/VHT)
#>         d/dt(Brain) = QBR * (Arterial_Blood/VAB - Brain/KbBR/VBR)
#>         d/dt(Muscles) = QMU * (Arterial_Blood/VAB - Muscles/KbMU/VMU)
#>         d/dt(Adipose) = QAD * (Arterial_Blood/VAB - Adipose/KbAD/VAD)
#>         d/dt(Skin) = QSK * (Arterial_Blood/VAB - Skin/KbSK/VSK)
#>         d/dt(Spleen) = QSP * (Arterial_Blood/VAB - Spleen/KbSP/VSP)
#>         d/dt(Pancreas) = QPA * (Arterial_Blood/VAB - Pancreas/KbPA/VPA)
#>         d/dt(Liver) = QHA * Arterial_Blood/VAB + QSP * Spleen/KbSP/VSP + 
#>             QPA * Pancreas/KbPA/VPA + QST * Stomach/KbST/VST + 
#>             QGU * Gut/KbGU/VGU - CLint * fub * Liver/KbLI/VLI - 
#>             QLI * Liver/KbLI/VLI
#>         d/dt(Stomach) = QST * (Arterial_Blood/VAB - Stomach/KbST/VST)
#>         d/dt(Gut) = QGU * (Arterial_Blood/VAB - Gut/KbGU/VGU)
#>         d/dt(Bones) = QBO * (Arterial_Blood/VAB - Bones/KbBO/VBO)
#>         d/dt(Kidneys) = QKI * (Arterial_Blood/VAB - Kidneys/KbKI/VKI)
#>         d/dt(Arterial_Blood) = QLU * (Lungs/KbLU/VLU - Arterial_Blood/VAB)
#>         d/dt(Venous_Blood) = QHT * Heart/KbHT/VHT + QBR * Brain/KbBR/VBR + 
#>             QMU * Muscles/KbMU/VMU + QAD * Adipose/KbAD/VAD + 
#>             QSK * Skin/KbSK/VSK + QLI * Liver/KbLI/VLI + QBO * 
#>             Bones/KbBO/VBO + QKI * Kidneys/KbKI/VKI + QRB * Rest_of_Body/KbRB/VRB - 
#>             QLU * Venous_Blood/VVB
#>         d/dt(Rest_of_Body) = QRB * (Arterial_Blood/VAB - Rest_of_Body/KbRB/VRB)
#>     })
#> }

You can also see this with the classic rxode2 model. In that case you use the summary() function:

pbpk <- pbpk$simulationModel
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
summary(pbpk)
#> rxode2 2.0.13.9000 model named rx_81d265cd14191d7bcaa87474579f4737 model (✔ ready). 
#> DLL: /tmp/RtmpbXPxP7/rxode2/rx_81d265cd14191d7bcaa87474579f4737__.rxd/rx_81d265cd14191d7bcaa87474579f4737_.so
#> NULL
#> 
#> Calculated Variables:
#>  [1] "KbBR"  "KbMU"  "KbAD"  "CLint" "KbBO"  "KbRB"  "CO"    "QHT"   "QBR"  
#> [10] "QMU"   "QAD"   "QSK"   "QSP"   "QPA"   "QLI"   "QST"   "QGU"   "QHA"  
#> [19] "QBO"   "QKI"   "QRB"   "QLU"   "VLU"   "VHT"   "VBR"   "VMU"   "VAD"  
#> [28] "VSK"   "VSP"   "VPA"   "VLI"   "VST"   "VGU"   "VBO"   "VKI"   "VAB"  
#> [37] "VVB"   "VRB"   "fub"   "KbLU"  "KbHT"  "KbSK"  "KbSP"  "KbPA"  "KbLI" 
#> [46] "KbST"  "KbGU"  "KbKI"  "S15"   "C15"  
#> ── rxode2 Model Syntax ──
#> rxode2({
#>     param(lKbBR, lKbMU, lKbAD, lCLint, eta.LClint, lKbBO, lKbRB, 
#>         WT)
#>     KbBR = exp(lKbBR)
#>     KbMU = exp(lKbMU)
#>     KbAD = exp(lKbAD)
#>     CLint = exp(lCLint + eta.LClint)
#>     KbBO = exp(lKbBO)
#>     KbRB = exp(lKbRB)
#>     CO = (187 * WT^0.81) * 60/1000
#>     QHT = 4 * CO/100
#>     QBR = 12 * CO/100
#>     QMU = 17 * CO/100
#>     QAD = 5 * CO/100
#>     QSK = 5 * CO/100
#>     QSP = 3 * CO/100
#>     QPA = 1 * CO/100
#>     QLI = 25.5 * CO/100
#>     QST = 1 * CO/100
#>     QGU = 14 * CO/100
#>     QHA = QLI - (QSP + QPA + QST + QGU)
#>     QBO = 5 * CO/100
#>     QKI = 19 * CO/100
#>     QRB = CO - (QHT + QBR + QMU + QAD + QSK + QLI + QBO + QKI)
#>     QLU = QHT + QBR + QMU + QAD + QSK + QLI + QBO + QKI + QRB
#>     VLU = (0.76 * WT/100)/1.051
#>     VHT = (0.47 * WT/100)/1.03
#>     VBR = (2 * WT/100)/1.036
#>     VMU = (40 * WT/100)/1.041
#>     VAD = (21.42 * WT/100)/0.916
#>     VSK = (3.71 * WT/100)/1.116
#>     VSP = (0.26 * WT/100)/1.054
#>     VPA = (0.14 * WT/100)/1.045
#>     VLI = (2.57 * WT/100)/1.04
#>     VST = (0.21 * WT/100)/1.05
#>     VGU = (1.44 * WT/100)/1.043
#>     VBO = (14.29 * WT/100)/1.99
#>     VKI = (0.44 * WT/100)/1.05
#>     VAB = (2.81 * WT/100)/1.04
#>     VVB = (5.62 * WT/100)/1.04
#>     VRB = (3.86 * WT/100)/1.04
#>     BP = 0.61
#>     fup = 0.028
#>     fub = fup/BP
#>     KbLU = exp(0.8334)
#>     KbHT = exp(1.1205)
#>     KbSK = exp(-0.5238)
#>     KbSP = exp(0.3224)
#>     KbPA = exp(0.3224)
#>     KbLI = exp(1.7604)
#>     KbST = exp(0.3224)
#>     KbGU = exp(1.2026)
#>     KbKI = exp(1.3171)
#>     S15 = VVB * BP/1000
#>     C15 = Venous_Blood/S15
#>     d/dt(Lungs) = QLU * (Venous_Blood/VVB - Lungs/KbLU/VLU)
#>     d/dt(Heart) = QHT * (Arterial_Blood/VAB - Heart/KbHT/VHT)
#>     d/dt(Brain) = QBR * (Arterial_Blood/VAB - Brain/KbBR/VBR)
#>     d/dt(Muscles) = QMU * (Arterial_Blood/VAB - Muscles/KbMU/VMU)
#>     d/dt(Adipose) = QAD * (Arterial_Blood/VAB - Adipose/KbAD/VAD)
#>     d/dt(Skin) = QSK * (Arterial_Blood/VAB - Skin/KbSK/VSK)
#>     d/dt(Spleen) = QSP * (Arterial_Blood/VAB - Spleen/KbSP/VSP)
#>     d/dt(Pancreas) = QPA * (Arterial_Blood/VAB - Pancreas/KbPA/VPA)
#>     d/dt(Liver) = QHA * Arterial_Blood/VAB + QSP * Spleen/KbSP/VSP + 
#>         QPA * Pancreas/KbPA/VPA + QST * Stomach/KbST/VST + QGU * 
#>         Gut/KbGU/VGU - CLint * fub * Liver/KbLI/VLI - QLI * Liver/KbLI/VLI
#>     d/dt(Stomach) = QST * (Arterial_Blood/VAB - Stomach/KbST/VST)
#>     d/dt(Gut) = QGU * (Arterial_Blood/VAB - Gut/KbGU/VGU)
#>     d/dt(Bones) = QBO * (Arterial_Blood/VAB - Bones/KbBO/VBO)
#>     d/dt(Kidneys) = QKI * (Arterial_Blood/VAB - Kidneys/KbKI/VKI)
#>     d/dt(Arterial_Blood) = QLU * (Lungs/KbLU/VLU - Arterial_Blood/VAB)
#>     d/dt(Venous_Blood) = QHT * Heart/KbHT/VHT + QBR * Brain/KbBR/VBR + 
#>         QMU * Muscles/KbMU/VMU + QAD * Adipose/KbAD/VAD + QSK * 
#>         Skin/KbSK/VSK + QLI * Liver/KbLI/VLI + QBO * Bones/KbBO/VBO + 
#>         QKI * Kidneys/KbKI/VKI + QRB * Rest_of_Body/KbRB/VRB - 
#>         QLU * Venous_Blood/VVB
#>     d/dt(Rest_of_Body) = QRB * (Arterial_Blood/VAB - Rest_of_Body/KbRB/VRB)
#> })

In this case, Venous_Blood is assigned to compartment 15. Figuring this out can be inconvenient and also lead to re-numbering compartment in simulation or estimation datasets. While it is easy and probably clearer to specify the compartment by name, other tools only support compartment numbers. Therefore, having a way to number compartment easily can lead to less data modification between multiple tools.

6.7.2 Changing compartments by pre-declaring with cmt()

To add the compartments to the rxode2 model in the order you desire you simply need to pre-declare the compartments with cmt. For example specifying is Venous_Blood and Skin to be the 1st and 2nd compartments, respectively, is simple:

pbpk2 <- function() {
  model({
    ## Now this is the first compartment, ie cmt=1
    cmt(Venous_Blood)
    ## Skin may be a compartment you wish to dose to as well,
    ##  so it is now cmt=2
    cmt(Skin)
    KbBR = exp(lKbBR)
    KbMU = exp(lKbMU)
    KbAD = exp(lKbAD)
    CLint= exp(lCLint + eta.LClint)
    KbBO = exp(lKbBO)
    KbRB = exp(lKbRB)

    ## Regional blood flows
    # Cardiac output (L/h) from White et al (1968)m
    CO  = (187.00*WT^0.81)*60/1000;
    QHT = 4.0 *CO/100;
    QBR = 12.0*CO/100;
    QMU = 17.0*CO/100;
    QAD = 5.0 *CO/100;
    QSK = 5.0 *CO/100;
    QSP = 3.0 *CO/100;
    QPA = 1.0 *CO/100;
    QLI = 25.5*CO/100;
    QST = 1.0 *CO/100;
    QGU = 14.0*CO/100;
    QHA = QLI - (QSP + QPA + QST + QGU); # Hepatic artery blood flow
    QBO = 5.0 *CO/100;
    QKI = 19.0*CO/100;
    QRB = CO - (QHT + QBR + QMU + QAD + QSK + QLI + QBO + QKI);
    QLU = QHT + QBR + QMU + QAD + QSK + QLI + QBO + QKI + QRB;

    ## Organs' volumes = organs' weights / organs' density
    VLU = (0.76 *WT/100)/1.051;
    VHT = (0.47 *WT/100)/1.030;
    VBR = (2.00 *WT/100)/1.036;
    VMU = (40.00*WT/100)/1.041;
    VAD = (21.42*WT/100)/0.916;
    VSK = (3.71 *WT/100)/1.116;
    VSP = (0.26 *WT/100)/1.054;
    VPA = (0.14 *WT/100)/1.045;
    VLI = (2.57 *WT/100)/1.040;
    VST = (0.21 *WT/100)/1.050;
    VGU = (1.44 *WT/100)/1.043;
    VBO = (14.29*WT/100)/1.990;
    VKI = (0.44 *WT/100)/1.050;
    VAB = (2.81 *WT/100)/1.040;
    VVB = (5.62 *WT/100)/1.040;
    VRB = (3.86 *WT/100)/1.040;

    ## Fixed parameters
    BP = 0.61;      # Blood:plasma partition coefficient
    fup = 0.028;    # Fraction unbound in plasma
    fub = fup/BP;   # Fraction unbound in blood

    KbLU = exp(0.8334);
    KbHT = exp(1.1205);
    KbSK = exp(-.5238);
    KbSP = exp(0.3224);
    KbPA = exp(0.3224);
    KbLI = exp(1.7604);
    KbST = exp(0.3224);
    KbGU = exp(1.2026);
    KbKI = exp(1.3171);

    ##-----------------------------------------
    S15 = VVB*BP/1000;
    C15 = Venous_Blood/S15

    ##-----------------------------------------
    d/dt(Lungs) = QLU*(Venous_Blood/VVB - Lungs/KbLU/VLU);
    d/dt(Heart) = QHT*(Arterial_Blood/VAB - Heart/KbHT/VHT);
    d/dt(Brain) = QBR*(Arterial_Blood/VAB - Brain/KbBR/VBR);
    d/dt(Muscles) = QMU*(Arterial_Blood/VAB - Muscles/KbMU/VMU);
    d/dt(Adipose) = QAD*(Arterial_Blood/VAB - Adipose/KbAD/VAD);
    d/dt(Skin) = QSK*(Arterial_Blood/VAB - Skin/KbSK/VSK);
    d/dt(Spleen) = QSP*(Arterial_Blood/VAB - Spleen/KbSP/VSP);
    d/dt(Pancreas) = QPA*(Arterial_Blood/VAB - Pancreas/KbPA/VPA);
    d/dt(Liver) = QHA*Arterial_Blood/VAB + QSP*Spleen/KbSP/VSP +
      QPA*Pancreas/KbPA/VPA + QST*Stomach/KbST/VST + QGU*Gut/KbGU/VGU -
      CLint*fub*Liver/KbLI/VLI - QLI*Liver/KbLI/VLI;
      d/dt(Stomach) = QST*(Arterial_Blood/VAB - Stomach/KbST/VST);
      d/dt(Gut) = QGU*(Arterial_Blood/VAB - Gut/KbGU/VGU);
      d/dt(Bones) = QBO*(Arterial_Blood/VAB - Bones/KbBO/VBO);
      d/dt(Kidneys) = QKI*(Arterial_Blood/VAB - Kidneys/KbKI/VKI);
      d/dt(Arterial_Blood) = QLU*(Lungs/KbLU/VLU - Arterial_Blood/VAB);
      d/dt(Venous_Blood) = QHT*Heart/KbHT/VHT + QBR*Brain/KbBR/VBR +
        QMU*Muscles/KbMU/VMU + QAD*Adipose/KbAD/VAD + QSK*Skin/KbSK/VSK +
        QLI*Liver/KbLI/VLI + QBO*Bones/KbBO/VBO + QKI*Kidneys/KbKI/VKI +
        QRB*Rest_of_Body/KbRB/VRB - QLU*Venous_Blood/VVB;
        d/dt(Rest_of_Body) = QRB*(Arterial_Blood/VAB - Rest_of_Body/KbRB/VRB);
  })
}

You can see this change in the simple printout

pbpk2 <- pbpk2()
pbpk2
#>  ── rxode2-based free-form 16-cmt ODE model ───────────────── 
#> 
#> States ($state or $stateDf): 
#>    Compartment Number Compartment Name
#> 1                   1     Venous_Blood
#> 2                   2             Skin
#> 3                   3            Lungs
#> 4                   4            Heart
#> 5                   5            Brain
#> 6                   6          Muscles
#> 7                   7          Adipose
#> 8                   8           Spleen
#> 9                   9         Pancreas
#> 10                 10            Liver
#> 11                 11          Stomach
#> 12                 12              Gut
#> 13                 13            Bones
#> 14                 14          Kidneys
#> 15                 15   Arterial_Blood
#> 16                 16     Rest_of_Body
#>  ── Model (Normalized Syntax): ── 
#> function() {
#>     model({
#>         cmt(Venous_Blood)
#>         cmt(Skin)
#>         KbBR = exp(lKbBR)
#>         KbMU = exp(lKbMU)
#>         KbAD = exp(lKbAD)
#>         CLint = exp(lCLint + eta.LClint)
#>         KbBO = exp(lKbBO)
#>         KbRB = exp(lKbRB)
#>         CO = (187 * WT^0.81) * 60/1000
#>         QHT = 4 * CO/100
#>         QBR = 12 * CO/100
#>         QMU = 17 * CO/100
#>         QAD = 5 * CO/100
#>         QSK = 5 * CO/100
#>         QSP = 3 * CO/100
#>         QPA = 1 * CO/100
#>         QLI = 25.5 * CO/100
#>         QST = 1 * CO/100
#>         QGU = 14 * CO/100
#>         QHA = QLI - (QSP + QPA + QST + QGU)
#>         QBO = 5 * CO/100
#>         QKI = 19 * CO/100
#>         QRB = CO - (QHT + QBR + QMU + QAD + QSK + QLI + QBO + 
#>             QKI)
#>         QLU = QHT + QBR + QMU + QAD + QSK + QLI + QBO + QKI + 
#>             QRB
#>         VLU = (0.76 * WT/100)/1.051
#>         VHT = (0.47 * WT/100)/1.03
#>         VBR = (2 * WT/100)/1.036
#>         VMU = (40 * WT/100)/1.041
#>         VAD = (21.42 * WT/100)/0.916
#>         VSK = (3.71 * WT/100)/1.116
#>         VSP = (0.26 * WT/100)/1.054
#>         VPA = (0.14 * WT/100)/1.045
#>         VLI = (2.57 * WT/100)/1.04
#>         VST = (0.21 * WT/100)/1.05
#>         VGU = (1.44 * WT/100)/1.043
#>         VBO = (14.29 * WT/100)/1.99
#>         VKI = (0.44 * WT/100)/1.05
#>         VAB = (2.81 * WT/100)/1.04
#>         VVB = (5.62 * WT/100)/1.04
#>         VRB = (3.86 * WT/100)/1.04
#>         BP = 0.61
#>         fup = 0.028
#>         fub = fup/BP
#>         KbLU = exp(0.8334)
#>         KbHT = exp(1.1205)
#>         KbSK = exp(-0.5238)
#>         KbSP = exp(0.3224)
#>         KbPA = exp(0.3224)
#>         KbLI = exp(1.7604)
#>         KbST = exp(0.3224)
#>         KbGU = exp(1.2026)
#>         KbKI = exp(1.3171)
#>         S15 = VVB * BP/1000
#>         C15 = Venous_Blood/S15
#>         d/dt(Lungs) = QLU * (Venous_Blood/VVB - Lungs/KbLU/VLU)
#>         d/dt(Heart) = QHT * (Arterial_Blood/VAB - Heart/KbHT/VHT)
#>         d/dt(Brain) = QBR * (Arterial_Blood/VAB - Brain/KbBR/VBR)
#>         d/dt(Muscles) = QMU * (Arterial_Blood/VAB - Muscles/KbMU/VMU)
#>         d/dt(Adipose) = QAD * (Arterial_Blood/VAB - Adipose/KbAD/VAD)
#>         d/dt(Skin) = QSK * (Arterial_Blood/VAB - Skin/KbSK/VSK)
#>         d/dt(Spleen) = QSP * (Arterial_Blood/VAB - Spleen/KbSP/VSP)
#>         d/dt(Pancreas) = QPA * (Arterial_Blood/VAB - Pancreas/KbPA/VPA)
#>         d/dt(Liver) = QHA * Arterial_Blood/VAB + QSP * Spleen/KbSP/VSP + 
#>             QPA * Pancreas/KbPA/VPA + QST * Stomach/KbST/VST + 
#>             QGU * Gut/KbGU/VGU - CLint * fub * Liver/KbLI/VLI - 
#>             QLI * Liver/KbLI/VLI
#>         d/dt(Stomach) = QST * (Arterial_Blood/VAB - Stomach/KbST/VST)
#>         d/dt(Gut) = QGU * (Arterial_Blood/VAB - Gut/KbGU/VGU)
#>         d/dt(Bones) = QBO * (Arterial_Blood/VAB - Bones/KbBO/VBO)
#>         d/dt(Kidneys) = QKI * (Arterial_Blood/VAB - Kidneys/KbKI/VKI)
#>         d/dt(Arterial_Blood) = QLU * (Lungs/KbLU/VLU - Arterial_Blood/VAB)
#>         d/dt(Venous_Blood) = QHT * Heart/KbHT/VHT + QBR * Brain/KbBR/VBR + 
#>             QMU * Muscles/KbMU/VMU + QAD * Adipose/KbAD/VAD + 
#>             QSK * Skin/KbSK/VSK + QLI * Liver/KbLI/VLI + QBO * 
#>             Bones/KbBO/VBO + QKI * Kidneys/KbKI/VKI + QRB * Rest_of_Body/KbRB/VRB - 
#>             QLU * Venous_Blood/VVB
#>         d/dt(Rest_of_Body) = QRB * (Arterial_Blood/VAB - Rest_of_Body/KbRB/VRB)
#>     })
#> }

The first two compartments are Venous_Blood followed by Skin.

6.7.3 Appending compartments to the model with cmt()

You can also append “compartments” to the model. Because of the ODE solving internals, you cannot add fake compartments to the model until after all the differential equations are defined.

For example this is legal:

ode.1c.ka <- function(){
  model({
    C2 = center/V
    d / dt(depot) = -KA * depot
    d/dt(center) = KA * depot - CL*C2
    cmt(eff)
  })
}

ode.1c.ka <- ode.1c.ka()
print(ode.1c.ka)
#>  ── rxode2-based free-form 2-cmt ODE model ────────────────── 
#> 
#> States ($state or $stateDf): 
#>   Compartment Number Compartment Name
#> 1                  1            depot
#> 2                  2           center
#>  ── Model (Normalized Syntax): ── 
#> function() {
#>     model({
#>         C2 = center/V
#>         d/dt(depot) = -KA * depot
#>         d/dt(center) = KA * depot - CL * C2
#>         cmt(eff)
#>     })
#> }

You can see this more clearly with the underlying classic rxode2 model:

ode.1c.ka$simulationModel
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
#> rxode2 2.0.13.9000 model named rx_5aa9aaa0d2b031c6e0671346ef0361f7 model (✔ ready). 
#> x$state: depot, center
#> x$stateExtra: eff
#> x$params: V, KA, CL
#> x$lhs: C2

But compartments defined before all the differential equations is not supported; So the model below:

ode.1c.ka <- rxode2({
    cmt(eff)
    C2 = center/V;
    d / dt(depot) = -KA * depot
    d/dt(center) = KA * depot - CL*C2
})

will give an error:

Error in rxModelVars_(obj) : 
  Evaluation error: Compartment 'eff' needs differential equations defined.