RxODE 包的使用 (CPT tutorial paper part 1)

本文基于tutorial paper : Wang W, Hallow KM, James DA. A Tutorial on RxODE: Simulating Differential Equation Pharmacometric Models in R. CPT Pharmacometrics Syst Pharmacol. 2016 Jan;5(1):3-10. doi: 10.1002/psp4.12052. Epub 2015 Dec 19. PMID: 26844010; PMCID: PMC4728294.

rxode2包是一个可以在R中模拟药代动力学模型的包，其具体作用总结如下图：

CASE STUDY

THIS CASE STUDY IS FROM THE PAPER
情况概述：
aim: adaptive dosing
PK: 2-comp
PD: indirect response
为了保持有效性同时降低毒副，对目标的抑制需保持 40-60%
具体模型结构如图：

1. 引入包

library(rxode2)
library(ggplot2)

NB：
To build models with rxode2, you need a working c compiler. To use parallel threaded solving in rxode2, this c compiler needs to support open-mp.
You can check to see if R has working c compiler you can check with: 如果返回True就没有问题

## install.packages("pkgbuild")
pkgbuild::has_build_tools(debug = TRUE)

## Found in Rtools 4.3 installation folder
## [1] TRUE

2. 定义模型结构

NB:
(1) 可以使用微分方程和代数式 both differential and algebraic equations are permitted
(2) 微分方程以 ‘d/dt(var_name)=’ 来表达
(3) 每个式子必须以分号 ; 分割。 Each equation is separated by a semicolon.
(4) 所有引用的未定义值都被当做输入参数 All referenced, undefined quantities are assumed to be input parameters.

# define model 
ode <- '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;' #每一个式子后都要加分号，包括最后一个式子

3. 定义模型参数和initial value

NB：
(1) 模型参数必须以命名向量定义。该命名向量中的参数名必须是ODE模型中的参数的超集（即ODE中所有参数都要在这个命名向量中定义，不可缺失）。命名向量中参数的顺序并不重要。
(2) Initial Condition (ic) 以向量定义。ic的数量必须与ode的数量一致，且顺序必须与ode的顺序一致。
原文：Model parameters are defined in a named vector. Names of parameters in the vector must be a superset of parameters in the ODE model, and the order of parameters within the vector is not important. Initial conditions (ICs) are defined through a vector as well. The number of ICs must equal exactly the number of ODEs in the model, and the order must be the same as the order in which the ODEs are listed in the model.

# Define parameters and initial conditions
params <- c(KA=0.3, CL=7, V2=40, Q=10, V3=300, Kin=0.2, Kout=0.2, EC50=8)
inits <- c(0, 0, 0, 1)

4. 定义 event Table

eventTable 由 RxODE 函数 eventTable()生成，用来描述dosing和sampling的事件。
add.dosing() 函数： 向eventTable里增加信息. Specified the dose amount, number of doses, dosing interval, compartment to dose into, rate (if an infusion), and dosing start time. More complex dosing schedules can be simulated by applying the add.dosing() function multiple times. The add.sampling() function allows specification of the time points to be included in the simulation output.
常见dosing 的例子：

在我们这里以简单的单次给药为例

ev <- eventTable()
# Specify dose
ev$add.dosing(dose = 10000, nbr.doses = 1)
# Specify sampling
ev$add.sampling(0:240) #这样采样240个点，ev有240行

head(ev,5)

##   id low time high       cmt   amt rate ii addl evid ss dur
## 1  1  NA    0   NA     (obs)    NA   NA NA   NA    0 NA  NA
## 2  1  NA    0   NA (default) 10000    0  0    0    1  0   0
## 3  1  NA    1   NA     (obs)    NA   NA NA   NA    0 NA  NA
## 4  1  NA    2   NA     (obs)    NA   NA NA   NA    0 NA  NA
## 5  1  NA    3   NA     (obs)    NA   NA NA   NA    0 NA  NA

5. 模拟

一、模拟
Calling the RxODE() function, the model is translated into C code, compiled, and dynamically loaded into the running R process. A simulation can then be performed by calling the R model object’s function run(), with the specified parameter vector, initial conditions vector, and event table as inputs.

# Compile model
mod1 <- RxODE(model=ode, modName='mod1')

## using C compiler: 'gcc.exe (GCC) 12.2.0'

# Run simulation
x <- mod1$run(params, ev, inits)

## Warning: Assumed order of inputs: depot, centr, peri, eff

head(x,5)

##      time       C2        C3     depot    centr      peri
## [1,]    0  0.00000 0.0000000 10000.000    0.000    0.0000
## [2,]    1 52.30679 0.9747541  7408.182 2092.272  292.4262
## [3,]    2 73.30489 3.0656620  5488.116 2932.196  919.6986
## [4,]    3 77.47397 5.4692230  4065.697 3098.959 1640.7669
## [5,]    4 73.25391 7.7752476  3011.942 2930.157 2332.5743
##           eff
## [1,] 1.000000
## [2,] 1.142853
## [3,] 1.315701
## [4,] 1.489110
## [5,] 1.659036

All state variables as well as other variables computed in the model are returned in the output matrix, at the times specified in the eventTable.
二、画图

x_df <- data.frame(x) #得到的x是"matrix" "array" 要转化成dataframe
ggplot(data = x_df, mapping = aes(x = time, y = C2)) + geom_line()

三、ODE solver
The user can also choose to specify the absolute and/or relative tolerance, as well as the type of solver to be used:

X <- m1$run(theta, ev, inits, stiff=F, atol=1e-8,rtol=1e-6)

RxODE uses the LSODA and a Runge-Kutta integrators for stiff and non-stiff equations, respectively. The LSODA (Livermore Solver for Ordinary Differential Equations) Fortran package is an automatic method switching for stiff and non-stiff problems throughout the integration interval. For purely non-stiff systems, RxODE uses DOP853, an explicit Runge-Kutta method of order 8(5,3).

6. variability的模拟

**得到参数矩阵：**RxODE提供了一种直接的方法来执行包含参数可变性和不确定性的模拟，即通过R采样函数（如rnorm和mvrnorm）生成一个参数值矩阵，其中每行表示一组参数值（如下列中params.all矩阵）。用户可以定义不确定性和变异的level。下面的代码为100个在CL和V2上具有相关的个体间变异的受试者参数。
原文：RxODE provides a straightforward way to perform simulations that incorporate parameter variability and uncertainty. A matrix of parameter values can be generated, where each row represents one set of parameter values. R sampling functions such as rnorm and mvrnorm can be utilized to build this matrix, but decisions on the level(s) of uncertainty to include are left to the discretion of the modeler, since this depends on the specific questions a given simulation is designed to address. The following code generates parameters for 100 subjects with correlated interindividual variability on CL and V2.

nsub <- 100# 100 subproblems
sigma <- matrix(c(0.09,0.08,0.08,0.25),2,2) # IIV covariance matrix
mv <- MASS::mvrnorm(n=nsub, rep(0,2), sigma) # Sample from covariance matrix
CL <- 7*exp(mv[,1])
V2 <- 40*exp(mv[,2])
params.all <- cbind(KA=0.3, CL=CL, V2=V2, Q=10, V3=300, Kin=0.2, Kout=0.2, EC50=8)

head(params.all,5)

##       KA        CL       V2  Q  V3 Kin Kout EC50
## [1,] 0.3  6.177585 37.42641 10 300 0.2  0.2    8
## [2,] 0.3  4.496332 47.95191 10 300 0.2  0.2    8
## [3,] 0.3 13.746675 64.89835 10 300 0.2  0.2    8
## [4,] 0.3  5.266878 32.39257 10 300 0.2  0.2    8
## [5,] 0.3  8.377355 89.71746 10 300 0.2  0.2    8

模拟
一旦生成了这个参数矩阵，就可以通过循环遍历参数矩阵来模拟每组参数（每行作为参数组输入）。
原文：Once this parameter matrix is generated, each subproblem can be simulated by looping through the parameter matrix, using each row as an input for the simulation, and collecting the output of each simulation in an output matrix.

res <- NULL
# Loop through each row of parameter values and simulation
for (i in 1:nsub) {
  params <- params.all[i,]
  x <- mod1$solve(params, ev, inits = inits)
  #Store results for effect compartment
  res <- cbind(res, x[, "eff"])
}

#用apply 函数简化代码
res2 <- apply(params.all, 1, function(theta) mod1$run(theta, ev, inits)[, 'eff'])

#画图
左图individual effect plots，右图median, 5th, and 95th percentile effect profile

# Plot results
par(mfrow = c(1,2), mar = c(4,4,1,1))
matplot(res, type = "l",  ylab = "Effect", xlab = "Time (hrs)")

# Calculate and plot quantiles 
res.q.t <- apply(res, 1, quantile, prob = c(.05, .5, .95))
matplot(t(res.q.t), type = "l", lty = c(2,1,2), col = c(2,1,2), ylab = "Effect", xlab = "Time (hrs)")
legend("topright", c("median","5 and 95%"), lty = c(1,2), col = c("black","red"), cex = .8)

软件版本

> devtools::session_info()
─ Session info ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.3.0 (2023-04-21 ucrt)
 os       Windows 10 x64 (build 19045)
 system   x86_64, mingw32
 ui       RStudio
 language (EN)
 collate  English_United States.utf8
 ctype    English_United States.utf8
 tz       Europe/Copenhagen
 date     2023-08-25
 rstudio  2023.06.1+524 Mountain Hydrangea (desktop)
 pandoc   3.1.1 @ C:/XXX

─ Packages ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
 package     * version date (UTC) lib source
 cachem        1.0.8   2023-05-01 [1] CRAN (R 4.3.1)
 callr         3.7.3   2022-11-02 [1] CRAN (R 4.3.1)
 cli           3.6.1   2023-03-23 [1] CRAN (R 4.3.1)
 crayon        1.5.2   2022-09-29 [1] CRAN (R 4.3.1)
 devtools      2.4.5   2022-10-11 [1] CRAN (R 4.3.1)
 digest        0.6.31  2022-12-11 [1] CRAN (R 4.3.1)
 ellipsis      0.3.2   2021-04-29 [1] CRAN (R 4.3.1)
 evaluate      0.21    2023-05-05 [1] CRAN (R 4.3.1)
 fastmap       1.1.1   2023-02-24 [1] CRAN (R 4.3.1)
 fs            1.6.2   2023-04-25 [1] CRAN (R 4.3.1)
 glue          1.6.2   2022-02-24 [1] CRAN (R 4.3.1)
 htmltools     0.5.5   2023-03-23 [1] CRAN (R 4.3.1)
 htmlwidgets   1.6.2   2023-03-17 [1] CRAN (R 4.3.1)
 httpuv        1.6.11  2023-05-11 [1] CRAN (R 4.3.1)
 knitr         1.43    2023-05-25 [1] CRAN (R 4.3.1)
 later         1.3.1   2023-05-02 [1] CRAN (R 4.3.1)
 lifecycle     1.0.3   2022-10-07 [1] CRAN (R 4.3.1)
 magrittr      2.0.3   2022-03-30 [1] CRAN (R 4.3.1)
 memoise       2.0.1   2021-11-26 [1] CRAN (R 4.3.1)
 mime          0.12    2021-09-28 [1] CRAN (R 4.3.0)
 miniUI        0.1.1.1 2018-05-18 [1] CRAN (R 4.3.1)
 pkgbuild      1.4.2   2023-06-26 [1] CRAN (R 4.3.1)
 pkgload       1.3.2   2022-11-16 [1] CRAN (R 4.3.1)
 prettyunits   1.1.1   2020-01-24 [1] CRAN (R 4.3.1)
 processx      3.8.1   2023-04-18 [1] CRAN (R 4.3.1)
 profvis       0.3.8   2023-05-02 [1] CRAN (R 4.3.1)
 promises      1.2.1   2023-08-10 [1] CRAN (R 4.3.1)
 ps            1.7.5   2023-04-18 [1] CRAN (R 4.3.1)
 purrr         1.0.1   2023-01-10 [1] CRAN (R 4.3.1)
 R6            2.5.1   2021-08-19 [1] CRAN (R 4.3.1)
 Rcpp          1.0.11  2023-07-06 [1] CRAN (R 4.3.1)
 remotes       2.4.2.1 2023-07-18 [1] CRAN (R 4.3.1)
 rlang         1.1.1   2023-04-28 [1] CRAN (R 4.3.1)
 rmarkdown     2.22    2023-06-01 [1] CRAN (R 4.3.1)
 rstudioapi    0.14    2022-08-22 [1] CRAN (R 4.3.1)
 sessioninfo   1.2.2   2021-12-06 [1] CRAN (R 4.3.1)
 shiny         1.7.5   2023-08-12 [1] CRAN (R 4.3.1)
 stringi       1.7.12  2023-01-11 [1] CRAN (R 4.3.0)
 stringr       1.5.0   2022-12-02 [1] CRAN (R 4.3.1)
 urlchecker    1.0.1   2021-11-30 [1] CRAN (R 4.3.1)
 usethis       2.2.2   2023-07-06 [1] CRAN (R 4.3.1)
 vctrs         0.6.3   2023-06-14 [1] CRAN (R 4.3.1)
 xfun          0.39    2023-04-20 [1] CRAN (R 4.3.1)
 xtable        1.8-4   2019-04-21 [1] CRAN (R 4.3.1)
 yaml          2.3.7   2023-01-23 [1] CRAN (R 4.3.0)