`RxODE-sim-var.Rmd`

In pharmacometrics the nonlinear-mixed effect modeling software (like nlmixr) characterizes the between-subject variability. With this between subject variability you can simulate new subjects.

Assuming that you have a 2-compartment, indirect response model, you can set create an RxODE model describing this system below:

The next step is to get the parameters into R so that you can start the simulation:

```
theta <- c(KA=2.94E-01, TCl=1.86E+01, V2=4.02E+01, # central
Q=1.05E+01, V3=2.97E+02, # peripheral
Kin=1, Kout=1, EC50=200, prop.err=0) # effects
```

In this case, I use `lotri`

to specify the omega since it uses similar lower-triangular matrix specification as nlmixr (also similar to NONMEM):

```
## the column names of the omega matrix need to match the parameters specified by RxODE
omega <- lotri(eta.Cl ~ 0.4^2)
omega
#> eta.Cl
#> eta.Cl 0.16
```

The next step to simulate is to create the dosing regimen for overall simulation:

If you wish, you can also add sampling times (though now RxODE can fill these in for you):

`ev <- ev %>% et(0,48, length.out=100)`

Note the `et`

takes similar arguments as `seq`

when adding sampling times. There are more methods to adding sampling times and events to make complex dosing regimens (See the event vignette). This includes ways to add variability to the both the sampling and dosing times).

Once this is complete you can simulate using the `rxSolve`

routine:

`sim <- rxSolve(mod,theta,ev,omega=omega,nSub=100)`

To quickly look and customize your simulation you use the default `plot`

routine. Since this is an RxODE object, it will create a `ggplot2`

object that you can modify as you wish. The extra parameter to the `plot`

tells `RxODE`

/`R`

what piece of information you are interested in plotting. In this case, we are interested in looking at the derived parameter `C2`

:

`plot`

```
library(ggplot2)
## The plots from RxODE are ggplots so they can be modified with
## standard ggplot commands.
plot(sim, C2, log="y") +
ylab("Central Compartment")
```

Of course this additional parameter could also be a state value, like `eff`

:

```
## They also takes many of the standard plot arguments; See ?plot
plot(sim, eff, log="y", ylab="Effect")
```

Or you could even look at the two side-by-side:

`plot(sim, C2, eff)`

Usually in pharmacometric simulations it is not enough to simply simulate the system. We have to do something easier to digest, like look at the central and extreme tendencies of the simulation.

Since the `RxODE`

solve object is a type of data frame

It is now straightforward to perform calculations and generate plots with the simulated data. You can

Below, the 5th, 50th, and 95th percentiles of the simulated data are plotted.

```
confint(sim, "C2", level=0.95) %>%
plot(ylab="Central Concentration", log="y")
#> Warning in confint.rxSolve(sim, "C2", level = 0.95): In order to put confidence
#> bands around the intervals, you need at least 2500 simulations.
#> summarizing data...done.
```

```
confint(sim, "eff", level=0.95) %>%
plot(ylab="Effect")
#> Warning in confint.rxSolve(sim, "eff", level = 0.95): In order to put confidence
#> bands around the intervals, you need at least 2500 simulations.
#> summarizing data...done.
```

Note that you can see the parameters that were simulated for the example

```
head(sim$param)
#> sim.id V2 prop.err V3 TCl eta.Cl KA Q Kin Kout EC50
#> 1 1 40.2 0 297 18.6 0.34990337 0.294 10.5 1 1 200
#> 2 2 40.2 0 297 18.6 0.17710829 0.294 10.5 1 1 200
#> 3 3 40.2 0 297 18.6 0.12633042 0.294 10.5 1 1 200
#> 4 4 40.2 0 297 18.6 0.05145749 0.294 10.5 1 1 200
#> 5 5 40.2 0 297 18.6 -0.04857253 0.294 10.5 1 1 200
#> 6 6 40.2 0 297 18.6 -0.91570711 0.294 10.5 1 1 200
```

In addition to conveniently simulating between subject variability, you can also easily simulate unexplained variability.

```
mod <- RxODE({
eff(0) = 1
C2 = centr/V2;
C3 = peri/V3;
CL = TCl*exp(eta.Cl) ## This is coded as a variable in the model
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;
e = eff+eff.err
cp = centr*(1+cp.err)
})
theta <- c(KA=2.94E-01, TCl=1.86E+01, V2=4.02E+01, # central
Q=1.05E+01, V3=2.97E+02, # peripheral
Kin=1, Kout=1, EC50=200) # effects
sigma <- lotri(eff.err ~ 0.1, cp.err ~ 0.1)
sim <- rxSolve(mod, theta, ev, omega=omega, nSub=100, sigma=sigma)
s <- confint(sim, c("eff", "centr"));
#> Warning in confint.rxSolve(sim, c("eff", "centr")): In order to put confidence
#> bands around the intervals, you need at least 2500 simulations.
#> summarizing data...done.
plot(s)
```

Sometimes you may want to match the dosing and observations of individuals in a clinical trial. To do this you will have to create a data.frame using the `RxODE`

event specification as well as an `ID`

column to indicate an individual. The RxODE event vignette talks more about how these datasets should be created.

```
library(dplyr)
#>
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#>
#> filter, lag
#> The following objects are masked from 'package:base':
#>
#> intersect, setdiff, setequal, union
ev1 <- eventTable(amount.units="mg", time.units="hours") %>%
add.dosing(dose=10000, nbr.doses=1, dosing.to=2) %>%
add.sampling(seq(0,48,length.out=10));
ev2 <- eventTable(amount.units="mg", time.units="hours") %>%
add.dosing(dose=5000, nbr.doses=1, dosing.to=2) %>%
add.sampling(seq(0,48,length.out=8));
dat <- rbind(data.frame(ID=1, ev1$get.EventTable()),
data.frame(ID=2, ev2$get.EventTable()))
## Note the number of subject is not needed since it is determined by the data
sim <- rxSolve(mod, theta, dat, omega=omega, sigma=sigma)
sim %>% select(id, time, e, cp)
#> id time e cp
#> 1 1 0.000000 [h] 0.8144742 7427.16624
#> 2 1 5.333333 [h] 0.9088491 733.04817
#> 3 1 10.666667 [h] 1.0700593 259.92205
#> 4 1 16.000000 [h] 0.6200344 223.39023
#> 5 1 21.333333 [h] 0.5376364 116.94388
#> 6 1 26.666667 [h] 1.5582777 178.91611
#> 7 1 32.000000 [h] 1.2177398 160.43714
#> 8 1 37.333333 [h] 1.0438762 92.62890
#> 9 1 42.666667 [h] 0.7487209 116.32714
#> 10 1 48.000000 [h] 0.4937465 107.91025
#> 11 2 0.000000 [h] 0.6631636 7416.25090
#> 12 2 6.857143 [h] 1.7473397 106.27233
#> 13 2 13.714286 [h] 1.6844069 69.48003
#> 14 2 20.571429 [h] 1.6304021 50.19524
#> 15 2 27.428571 [h] 0.8856005 69.30219
#> 16 2 34.285714 [h] 0.5952504 50.20523
#> 17 2 41.142857 [h] 0.8388760 55.47943
#> 18 2 48.000000 [h] 0.8981448 45.01631
```

By either using a simple single event table, or data from a clinical trial as described above, a complete clinical trial simulation can be performed.

Typically in clinical trial simulations you want to account for the uncertainty in the fixed parameter estimates, and even the uncertainty in both your between subject variability as well as the unexplained variability.

`RxODE`

allows you to account for these uncertainties by simulating multiple virtual “studies,” specified by the parameter `nStud`

. In a single virtual study:

A Population effect parameter is sampled from a multivariate normal distribution with mean given by the parameter estimates and the variance specified by the named matrix

`thetaMat`

.A between subject variability/covariance matrix is sampled from either a scaled inverse chi-squared distribution (for the univariate case) or a inverse Wishart that is parameterized to scale to the conjugate prior covariance term, as described by the wikipedia article. (This is not the same as the scaled inverse Wishart distribution ). In the case of the between subject variability, the variance/covariance matrix is given by the ‘omega’ matrix parameter and the degrees of freedom is the number of subjects in the simulation.

Unexplained variability is also simulated from the scaled inverse chi squared distribution or inverse Wishart distribution with the variance/covariance matrix given by the ‘sigma’ matrix parameter and the degrees of freedom given by the number of observations being simulated.

The covariance/variance prior is simulated from `RxODE`

s `cvPost()`

function.

An example of this simulation is below:

```
## Creating covariance matrix
tmp <- matrix(rnorm(8^2), 8, 8)
tMat <- tcrossprod(tmp, tmp) / (8 ^ 2)
dimnames(tMat) <- list(NULL, names(theta))
sim <- rxSolve(mod, theta, ev, omega=omega, nSub=100, sigma=sigma, thetaMat=tMat, nStud=10,
dfSub=10, dfObs=100)
s <-sim %>% confint(c("centr", "eff"))
#> Summarizing data
#> done.
plot(s)
```

If you wish you can see what `omega`

and `sigma`

was used for each virtual study by accessing them in the solved data object with `$omega.list`

and `$sigma.list`

:

```
head(sim$omega.list)
#> [[1]]
#> [,1]
#> [1,] 0.1375783
#>
#> [[2]]
#> [,1]
#> [1,] 0.219712
#>
#> [[3]]
#> [,1]
#> [1,] 0.07645734
#>
#> [[4]]
#> [,1]
#> [1,] 0.1357399
#>
#> [[5]]
#> [,1]
#> [1,] 0.2711316
#>
#> [[6]]
#> [,1]
#> [1,] 0.1323411
```

```
head(sim$sigma.list)
#> [[1]]
#> [,1] [,2]
#> [1,] 0.115063832 0.003474916
#> [2,] 0.003474916 0.112844178
#>
#> [[2]]
#> [,1] [,2]
#> [1,] 0.09346365 0.01247610
#> [2,] 0.01247610 0.08278152
#>
#> [[3]]
#> [,1] [,2]
#> [1,] 0.097531442 0.000516414
#> [2,] 0.000516414 0.123232919
#>
#> [[4]]
#> [,1] [,2]
#> [1,] 0.088659267 0.009842486
#> [2,] 0.009842486 0.089518635
#>
#> [[5]]
#> [,1] [,2]
#> [1,] 0.09830986 -0.02394862
#> [2,] -0.02394862 0.11792143
#>
#> [[6]]
#> [,1] [,2]
#> [1,] 0.093216731 -0.001967668
#> [2,] -0.001967668 0.106582003
```

You can also see the parameter realizations from the `$params`

data frame.

If you do not wish to sample from the prior distributions of either the `omega`

or `sigma`

matrices, you can turn off this feature by specifying the `simVariability = FALSE`

option when solving:

```
sim <- rxSolve(mod, theta, ev, omega=omega, nSub=100, sigma=sigma, thetaMat=tMat, nStud=10,
simVariability=FALSE);
s <-sim %>% confint(c("centr", "eff"))
#> Summarizing data
#> done.
plot(s)
```

Note since realizations of `omega`

and `sigma`

were not simulated, `$omega.list`

and `$sigma.list`

both return `NULL`

.

RxODE now supports multi-threaded solving on OpenMP supported compilers, including linux and windows. Mac OSX can also be supported By default it uses all your available cores for solving as determined by `rxCores()`

. This may be overkill depending on your system, at a certain point the speed of solving is limited by things other than computing power.

You can also speed up simulation by using the multi-cores to generate random deviates with `mvnfast`

(either `mvnfast::rmvn()`

or `mvnfast::rmvt()`

). This is controlled by the `nCoresRV`

parameter. For example:

```
sim <- rxSolve(mod, theta, ev, omega=omega, nSub=100, sigma=sigma, thetaMat=tMat, nStud=10,
nCoresRV=2);
s <-sim %>% confint(c("eff", "centr"))
#> Summarizing data
#> done.
```

The default for this is `1`

core since the result depends on the number of cores and the random seed you use in your simulation. However, you can always speed up this process with more cores if you are sure your collaborators have the same number of cores available to them and have OpenMP thread-capable compile.