07_single_subjects.qmd

# Getting started with `rxode2`: simulating single subjects {#rxode2}

## Basic concepts of `rxode2` syntax

Now that we've got to grips with events, let's have a look at the models themselves!

### Writing models for `rxode2`

Writing models for `rxode2` (and by extension `nlmixr2`) is relatively straightforward. Here's a simple example.

```{r}

library(rxode2)
library(tidyverse)

mod1 <- function() {
  ini({
    # central compartment
    KA   = 0.294  # /h
    CL   = 18.6   # L/h
    V2   = 40.2   # L
    
    # peripheral compartment
    Q    = 10.5   # L/h
    V3   = 297    # L
    
    # effects
    Kin  = 1      # Effect generation rate constant (/h)
    Kout = 1      # Effect elimination rate constant (/h)
    EC50 = 200    # EC50 (ug/ml)
  })
  model({
    C2          <- centr/V2   # concentration in the central compartment
    C3          <- peri/V3    # concentration in the peripheral compartment
    
    d/dt(depot) <- -KA*depot                         # depot compartment
    d/dt(centr) <-  KA*depot - CL*C2 - Q*C2 + Q*C3   # central compartment
    d/dt(peri)  <-                     Q*C2 - Q*C3   # peripheral compartment
    
    eff(0)      <- 1          # amount in the effect compartment at time 0
    d/dt(eff)   <- Kin - Kout*(1-C2/(EC50+C2))*eff   # effect compartment
  })
}
```

Yes, yes, I know we said simple. There's a lot to unpack here! This code represents a two-compartmental PK model with first-order absorption and linear elimination.

```{dot}

digraph Model {
  overlap=False
  
  node [shape=box]; Depot [label="Gut (depot)"]; 
    V2 [label="Central compartment\n(centr)"]; 
    V3 [label="Peripheral compartment\n(peri)"];
    Excreted [color=white, fontcolor=blue, style=filled];
    In [color=white, fontcolor=white, style=filled];
    Out [color=white, fontcolor=white, style=filled];
    Effect [label="Effect\n(eff)", color=orange, style=filled];
    ab[label="", fixedsize="false", width=0, height=0, shape=none];
    za[label="", fixedsize="false", width=0, height=0, shape=none];
    
  node [shape=circle]; Dose [color=red, fontcolor=white, style=filled]; 

  In -> Effect [label="kin"];
  Effect -> ab [arrowhead=None, label="kout"];
  ab -> Out;
  Dose -> Depot; 
  Depot -> V2 [label=" ka ", fontsize=10]; 
  V2 -> V3 [label=" k23 =\nQ/V2 ", fontsize=10];
  V3 -> V2 [label=" k32 =\nQ/V3", fontsize=10];
  V2 -> Excreted [label="k = CL/V2 ", fontsize=10];
  V2 -> ab [style=dashed, color=orange, arrowhead=tee, label=" C2=centr/V2", fontcolor=orange];

  
  {rank=same; V2; V3; Excreted};
  {rank=same; za; In; Effect; ab; Out};
}
  
```

An `rxode2` model specification consists of a series of one or more statements, optionally terminated by semicolons, and comments (delimited by \# and an end-of-line). Comments are also optional.

### Model blocks

`rxode2` models are divided into two discrete blocks: `ini()`, which sets initial conditions, and `model()`, which defines the model and other aspects of the system that change with time. Blocks of statements are delimited by curly braces. The `ini()` section is pretty straightforward, as you see, with important variables in the system being assigned their starting values. You can also include a variety of other statements (such as conditionals, `while` loops and print commands) if you wish.

The `model()` block is where the action is. Here, as you might imagine, we define the model: the components of the system that change with the independent variable (often this is time). In our example, we are doing the following:

-   Defining concentrations `C2` and `C3` in the central (`centr`) and peripheral (`peri`) PK compartments, respectively
-   Defining the differential equations for the depot (`depot`), `centr`, `peri` and effect (`eff`) compartments
-   Defining the starting amount in `eff` (1 at time 0)

### Making statements: `rxode2` nomenclature and syntax

Before we go any further, it's probably useful to spend a bit of time talking about assignments, nomenclature and syntax in `rxode2`.

#### Assignments and operators

Assignment statements can be:

-   *simple assignments*, in which the left-hand side is an identifier (a variable)
-   *time-derivative assignments*, where the left-hand side 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* in which the left-hand side specifies the compartment of the initial condition being specified, e.g. `depot(0) = 0`
-   special *model event changes* such as bioavailability (e.g. `f(depot) = 1`), lag time (e.g. `alag(depot) = 0`), modeled rate (e.g. `rate(depot) = 2`) and modeled duration of infusion (e.g. `dur(depot) = 2`)
-   special *change-point syntax*, or *modeled event times*, e.g `mtime(var) = time`
-   *Jacobian-derivative assignments*, in which 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 slowed things down.

Assignments can be made using `=`, `<-` or `~`. When using the `~` operator, 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 `nlmixr2` models); these can be specified using `cmt(compartmentName)`
-   *Parameter declaration statements*, which can be used to ensure the input parameters are kept in a certain order instead of ordering the parameters by the order in which they are parsed. This is useful for keeping the parameter order the same when using different ODE models, for example (e.g. `param(par1, par2,...)`)

Expressions in assignment and in conditional (`if`) statements can be numeric or logical.

Numeric expressions can include the standard numeric operators (`+`, `-`, `*`, `/`, `^`) as well as mathematical functions defined in the C or the R math libraries (e.g. `fabs`, `exp`, `log`, `sin`, `abs`). You may also access R's math functions, like `lgammafn` for the log gamma function.

`rxode2` syntax is case-sensitive, like the rest of R: `ABC` is different from `abc`, `Abc`, `ABc`, and so forth.

#### Identifiers

As in R, identifiers (variable names) may consist of one or more alphanumeric, underscore (`_`) or period (`.`) characters, although the first character cannot be a digit or underscore.

Identifiers in a model specification can take the following forms:

-   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 case of absorption transit models)
-   Special constants such as `pi` or R's predefined constants
-   Model parameters (such as rate of absorption and clearance)
-   Other left-hand side (LHS) variables created by assignments as part of the model specification

Currently, the `rxode2` modeling language only recognizes system state variables and parameters. Any values that need to be passed from R to the ODE model (such as covariates) must be passed in the `params` argument of the integrator function (`rxSolve()`) or be available in the supplied event dataset, which we'll get to in a bit.

Some variable names are reserved for use in `rxode2` event tables. The following items cannot be assigned, or used as a state, but can be accessed in the `rxode2` code:

-   `cmt`: compartment
-   `dvid`: dependent variable ID
-   `addl`: number of additional doses
-   `ss`: steady state
-   `rate`: infusion rate
-   `id`: unique subject identifier

The following variables, however, cannot be used in a model specification in any way:

-   `evid`: event type
-   `ii`: interdose interval

`rxode2` generates variables which are used internally. These variables start with an `rx` prefix. To avoid any problems, it is *strongly* suggested not to use an `rx` prefix when writing model code, since all kinds of unpleasant and unpredictable things may happen.

#### Logical Operators

Logical operators support the standard R operators (`==`, `!=`, `>=`, `<=`, `>` and `<`). As in R, these can be used in `if()`, `while()` and `ifelse()` expressions, as well as in standard assignments. 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, however, since character comparisons are slower than numeric expressions. Unlike R, `as.numeric()` or `as.integer()` for logical statements are not only not needed, but not permitted - they will throw an error if you try to use them.

### Interface and data handling between R and the generated C code

Users define the dynamic system's state variables via the `d/dt(identifier)` statements as part of the model specification, and model parameters via the `params` argument in the `rxode2` `solve()` method:

```{r eval=FALSE}

m1 <- rxode2(model = ode, modName = "m1")

# model parameters -- a named vector is required

theta <- c(KA=0.29, CL=18.6, V2=40.2, Q=10.5, V3=297, Kin=1, Kout=1, EC50=200)

# state variables and their amounts at time 0 (the use of names is
# encouraged, but not required)

inits <- c(depot=0, centr=0, peri=0, eff=1)

# qd1 is an eventTable specification with a set of dosing and sampling
# records (code not shown here)

solve(theta, event = qd1, inits = inits) 

```

The values of these variables at pre-specified time points are saved during model fitting/integration and returned with the fitted values as part of the modeling output (see the function `eventTable()`, and in particular its member function `add.sampling()` for further information on defining a set of time points at which to capture the values of these variables).

The ODE specification mini-language is parsed with the help of the open source tool `DParser` [@Plevyak].

### Supported functions

All the supported functions in `rxode2` can be seen with the function `rxSupportedFuns()`. There are a lot of them.

## Other model types

`rxode2` supports a range of model types. So far, we've just looked at ODE systems.

### "Prediction-only" models

"Prediction-only" models are analogous to \$PRED models in NONMEM, which some readers may be familiar with. Here's a very simple example - a one-compartment model with bolus dosing.

```{r pred1}

mod <- rxode2({
    ipre <- 10 * exp(-ke * t);
})

```

Solving prediction-only models is done just the same way as for ODE systems, but is faster.

```{r pred2}

et   <- et(seq(0, 24, length.out=50))
cmt1 <- rxSolve(mod, et, params = c(ke=0.5))
cmt1

```

```{r pred3}

library(ggplot2)
library(patchwork)
ggplot(cmt1, aes(time, ipre)) +
  geom_line(col="red") +
  theme_light()

```

### Solved systems

`rxode2` has a library of models which have been pre-solved, and thus do not need to be defined in terms of ODEs. These can be used by including `linCmt()` function in model code, along with a set of model parameters that fits the model you wish to use. For example:

| *Model*                             | *Parameters*                              | *Microconstants*                              | *Hybrid constants*                      |
|------------------|------------------|------------------|------------------|
| One-compartment, bolus or IV dose   | `cl`, `v`                                 | `v`, `ke`                                     | `v`, `alpha`                            |
| One-compartment, oral dose          | `cl`, `v`, `ka`                           | `v`, `ke`, `ka`                               | `v`, `alpha`, `ka`                      |
| Two-compartment, bolus or IV dose   | `cl`, `v1`, `v2`, `q`                     | `v1`, `k12`, `k21`, `ke`                      | `v`, `alpha`, `beta`, `aob`             |
| Two-compartment, oral dose          | `cl`, `v1`, `v2`, `q`, `ka`               | `v1`, `k12`, `k21`, `ke`, `ka`                | `v`, `alpha`, `beta`, `aob`, `ka`       |
| Three-compartment, bolus or IV dose | `cl`, `vc`, `vp`, `vp2`, `q1`, `q2`       | `v1`, `k12`, `k21`, `k13`, `k31`,`ke`         | `a`, `alpha`, `b`, `beta`, `c`, `gamma` |
| Three-compartment, oral dose        | `cl`, `vc`, `vp`, `vp2`, `q1`, `q2`, `ka` | `a`, `alpha`, `b`, `beta`, `c`, `gamma`, `ka` |                                         |

Most parameters can be specified in a number of ways. For example, central volume in a one-compartment model can be called `v`, `v1` or `vc`, and will be understood as the same parameter by `linCmt()`. Here's a very simple example...

```{r solved}

mod_solved <- rxode2({
    ke <- 0.5
    V <- 1
    ipre <- linCmt();
})

mod_solved


```

We can treat this the same way as an ODE model:

```{r solved2}

et  <- et(amt=10, time=0, cmt=depot) %>%
    et(seq(0, 24, length.out=50))
cmt_solved <- rxSolve(mod_solved, et, params=c(ke=0.5))
cmt_solved

```

### Combining ODEs and solved systems

Solved systems and ODEs can be combined. Let's look at our two-compartment indirect-effect model again.

```{r, warning=FALSE,message=FALSE}

## Set up parameters and initial conditions

theta <- c(KA   = 0.294, 
           CL   = 18.6, 
           V2   = 40.2, 
           Q    = 10.5,  
           V3   = 297,
           Kin  = 1, 
           Kout = 1, 
           EC50 = 200)

inits <- c(eff = 1);

## Set up dosing event information

ev <- eventTable(amount.units='mg', time.units='hours') %>%
    add.dosing(dose=10000, nbr.doses=10, dosing.interval=12, dosing.to=1) %>%
    add.dosing(dose=20000, nbr.doses=5, start.time=120,dosing.interval=24, dosing.to=1) %>%
    add.sampling(0:240);

## Set up a mixed solved/ODE system

mod2 <- rxode2({
  ## the order of variables do not matter
  ## the type of compartmental model is determined by the parameters specified
  
    C2        = linCmt(KA, CL, V2, Q, V3);  # you don't need to provide the parameters, but you can if you want to
    eff(0)    = 1  ## The initial amount in the effect compartment is 1
    d/dt(eff) =  Kin - Kout*(1 - C2/(EC50 + C2))*eff;
})

```

```{r}

mod2

```

Concentration output from the 2-compartment model is assigned to the `C2` variable and is subsequently used in the indirect response model.

Ntoe that when mixing solved systems and ODEs, the solved system's "compartment" is always the last one. This is because the solved system technically isn't a compartment as such. Adding the dosing compartment to the end will not interfere with the actual ODE to be solved.

In this example, therefore, the effect compartment is compartment #1 while the PK dosing compartment for the depot is compartment #2.

Let's solve the system and see what it looks like.

```{r}

x <- mod2 %>%  solve(theta, ev)
print(x)

```

```{r}

p1 <- ggplot(x, aes(time, C2)) + 
  geom_line(col="red") +
  theme_light() +
  labs(title="Concentration")

p2 <- ggplot(x, aes(time, eff)) +   
  geom_line(col="red") +
  theme_light()+
  labs(title="Effect")

p1 + p2 + plot_layout(nrow=2)

```

## Compartment numbers

rxode2 automatically assigns compartment numbers when parsing. For illustrative purposes, elt's look at something a bit more compex than we have so far: a PBPK model for mavoglurant published by Wendling and colleagues [@Wendling2016].

```{r}
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)
  })
}
```

This is quite a meaty model, with 16 compartments linked by ODEs.

![](Wendling2016.jpg)

### How `rxode2` assigns compartment numbers

```{r}
pbpk <- pbpk()
print(pbpk)
```

Here, `Venous_Blood` is assigned to compartment 15. Keeping track of compartment numbers in large models like this can be inconvenient and challenging, and more importantly, can lead to mistakes. While it is easy, and probably clearer, to specify the compartments by name rather than number, other pharmacometric software in common use only supports compartment numbers. Having a way to number compartments easily can therefore be handy when moving between different tools.

### Pre-declaring the compartments

Pre-assigning compartment numbers can be helpful in this situation.

To add the compartments to the model in the order desired, we can pre-declare them with `cmt`. For example, specifying `Venous_Blood` and `Skin` as the first and second compartments, respectively, is pretty straightforward:

```{r}
pbpk2 <- function() {
  model({
    
    cmt(Venous_Blood)  ## Now this is the first compartment, ie cmt=1
    cmt(Skin)          ## Now this is the second compartment, ie cmt=2
    
    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);
  })
}
```

```{r}
pbpk2 <- pbpk2()
pbpk2
```

Now `Venous_Blood` and `Skin` are where we want them.

### Appending compartments

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:

```{r}
mod2 <- function(){
  model({
    C2 = center/V
    d / dt(depot) = -KA * depot
    d/dt(center) = KA * depot - CL*C2
    cmt(eff)
  })
}

mod2 <- mod2()
print(mod2)
```

You can see this more clearly by querying the `simulationModel` property:

```{r}
mod2$simulationModel
```

Defining "extra" compartments before the differential equations is not supported. The model below will throw an error if executed.

```         
mod2 <- rxode2({
    cmt(eff)
    C2 = center/V;
    d / dt(depot) = -KA * depot
    d/dt(center) = KA * depot - CL*C2
})
```

## Transit compartments

Transit compartments provide a useful way to better approximate absorption lag times, without the numerical problems and stiffness associated with the conventional way of modeling these [@Savic2007; @Wilkins2008]. `rxode2` has them built in.

The transit compartment function (`transit`) can be used to specify the model without having to write it out in full (although you could do that too, if you wanted to). `transit()` takes parameters corresponding to the number of transit compartments (`n` in the code below), the mean transit time (`mtt`), and bioavailability (`bio`, which is optional).

```{r}
mod <- function() {
  ini({
    ## Table 3 from Savic 2007
    cl  <- 17.2 # (L/hr)
    vc  <- 45.1 # L
    ka  <- 0.38 # 1/hr
    mtt <- 0.37 # hr
    bio <- 1
    n   <- 20.1
  })
  model({
    k           <- cl/vc
    ktr         <- (n+1)/mtt
    d/dt(depot) <- transit(n,mtt,bio)-ka*depot
    # or alternately -
    # d/dt(depot) <- exp(log(bio*podo(depot))+log(ktr)+n*log(ktr*tad(depot))-
    #                    ktr*tad(depot)-lgammafn(n+1))-ka*depot
    d/dt(cen)   <- ka*depot-k*cen
  })
}

et <- et(0, 7, length.out=200) %>%
  et(amt=20, evid=7)

transit <- rxSolve(mod, et)

ggplot(transit, aes(time, cen)) +
  geom_line(col="red") +
  xlab("Time") + ylab("Concentration") +
  theme_light()

```

A couple of things to keep in mind when using this approach:

-   This approach implicitly assumes that the absorption through the transit compartment is completed before the next dose is given If this isn't the case, some additional code is needed to correct for this [@Wilkins2008].

-   Different dose types (bolus or infusion) to the `depot` compartment affect the time after dose calculation (`tad`) which is used in the transit compartment code. Direct doses into compartments are therefore not currently supported. The most stable way around this is to use `tad(cmt)` and `podo(cmt)` - this way, doses to other compartments do not affect the transit compartment machinery.

-   Internally, the `transit` syntax uses either the currently defined compartment, `d/dt(cmt)=transit(...)`, or `cmt`. If the transit compartment is used outside of a `d/dt()` (not recommended), the `cmt` that is used is the last `d/dt(cmt)` defined it the model. This also means compartments do not affect one another (ie a oral, transit compartment drug dosed immediately with an IV infusion)

## Covariates

## Multiple subjects

## Working with `rxode2` output

## Piping

## Special cases

### Jacobian solving

### `rxode2` models in `shiny`

### Precompiled `rxode2` models in other R packages