Introduction

This is a simple example on how to couple PopED with external R based PKPD simulation tools. Typically, these tools might be R packages that can simulate from hierarchical, ordinary differential equation (ODE) based models. In this document you will see how to couple PopED to models, defined with ODEs, implemented using:

  • deSolve (with native R ODE models)
  • deSolve (with compiled C ODE models)
  • deSolve (with compiled C++ ODE models using Rcpp)
  • PKPDsim
  • mrgsolve
  • RxODE

For the future we will also show how to use

  • PKPDmodels
  • dMod

The model

In this example we will use a one-compartment linear absorption model

This model can be described with the following set of ODEs:

\[ \begin{split} \frac{dA_{0}}{dt} &= - k_{a} \cdot A_{0}\\ \frac{dA_{1}}{dt} &=-(CL/V_1)\cdot A_1 + k_{a} \cdot A_{0} \\ f(t) &= A_1/V_1 \end{split} (\#eq:ode) \]

All compartment amounts are assumed to be zero at time zero (\(\boldsymbol{A}[t=0]=0\)). Inputs to the system come in tablet form and are added to the amount in \(A_{0}\) according to

\[ \text{Input}(t,D,t_D) = \begin{cases} D, &\text{if} \quad t = t_D \\ 0, &\text{otherwise} \end{cases} (\#eq:input) \]

Parameter values are defined as:

\[ \begin{split} k_a &= \theta_1 \cdot e^{\eta_1} \\ CL &= \theta_2 \cdot e^{\eta_2} \\ V_1 &= \theta_3 \cdot e^{\eta_3} \\ \end{split} (\#eq:par) \] where elements of the between subject variability (BSV), \(\eta_{j}\), vary across individuals and come from normal distributions with means of zero and variances of \(\omega^2_{j}\).

The residual unexplained variability (RUV) model has a proportional and additive component

\[ y = A_1/V_1 \cdot (1+\varepsilon_1) + \varepsilon_2 (\#eq:ruv) \]

elements of \(\boldsymbol{\varepsilon}_{j}\) vary accross observations and come from normal distributions with means of zero and variances of \(\sigma^2_{j}\).

Parameter values are shown in Table @ref(tab:values).
Model parameter values.
Parameter Value
\(k_a\) 0.2500
CL 3.7500
\(V_1\) 72.8000
\(\omega^2_{k_a}\) 0.0900
\(\omega^2_{CL}\) 0.0625
\(\omega^2_{V_1}\) 0.0900
\(\sigma^2_{prop}\) 0.0400
\(\sigma^2_{add}\) 0.0025

Model implementation

Below we implement this model using a number of different methods. For the ODE solvers, if possible, we set the tuning parameters to be the same values (atol, rtol, etc.).

Analytic solution

First we implement an analytic solution to the model in a function that could be used in PopED.

ff <- function(model_switch,xt,parameters,poped.db){
  with(as.list(parameters),{
    y=xt
    N = floor(xt/TAU)+1
    y=(DOSE/V)*(KA/(KA - CL/V)) * 
      (exp(-CL/V * (xt - (N - 1) * TAU)) * (1 - exp(-N * CL/V * TAU))/(1 - exp(-CL/V * TAU)) - 
         exp(-KA * (xt - (N - 1) * TAU)) * (1 - exp(-N * KA * TAU))/(1 - exp(-KA * TAU)))  
    return(list( y=y,poped.db=poped.db))
  })
}

ODE solution using deSolve

Model using ODEs defined in deSolve

PK_1_comp_oral_ode <- function(Time, State, Pars){
  with(as.list(c(State, Pars)), {    
    dA1 <- -KA*A1
    dA2 <- KA*A1 - (CL/V)*A2
    return(list(c(dA1, dA2)))
  })
}
ff_ode_desolve <- function(model_switch, xt, parameters, poped.db){
  with(as.list(parameters),{
    A_ini <- c(A1=0, A2=0)
    
    #Set up time points for the ODE
    times_xt <- drop(xt)
    times <- c(0,times_xt) ## add extra time for start of integration
    dose_times = seq(from=0,to=max(times_xt),by=TAU)
    times <- c(times,dose_times)
    times <- sort(times) 
    times <- unique(times) # remove duplicates
    
    eventdat <- data.frame(var = c("A1"), 
                           time = dose_times,
                           value = c(DOSE), method = c("add"))
    
    out <- ode(A_ini, times, PK_1_comp_oral_ode, parameters, 
               events = list(data = eventdat),
               atol=1e-8, rtol=1e-8,maxsteps=5000)
    
    # grab timepoint values
    out = out[match(times_xt,out[,"time"]),]
    
    y = out[,"A2"]/V
    
    y=cbind(y) # must be a column matrix 
    return(list(y=y,poped.db=poped.db))
  })
}

ODE solution using deSolve and compiled C code

Model using ODEs defined in deSolve with compiled code

cat(readChar(system.file("examples/one_comp_oral_CL.c", package="PopED"), 1e5))
#> /* file tmdd_qss_one_target.c */
#> #include <R.h>
#> static double parms[3];
#> #define CL parms[0]
#> #define V parms[1]
#> #define KA parms[2]
#> 
#> /* initializer  */
#> void initmod(void (* odeparms)(int *, double *))
#> {
#>   int N=3;
#>   odeparms(&N, parms);
#> }
#> 
#> /* Derivatives and 1 output variable */
#> void derivs (int *neq, double *t, double *y, double *ydot,
#>       double *yout, int *ip)
#> {
#>     
#>   if (ip[0] <1) error("nout should be at least 1");
#>     
#>   ydot[0] = -KA*y[0];
#>   ydot[1] = KA*y[0] - CL/V*y[1];
#>   yout[0] = y[0]+y[1];
#> }
#> 
#> /* END file tmdd_qss_one_target.c */
file.copy(system.file("examples/one_comp_oral_CL.c", package="PopED"),"./one_comp_oral_CL.c")
#> [1] FALSE
system('R CMD SHLIB one_comp_oral_CL.c')
dyn.load(paste("one_comp_oral_CL", .Platform$dynlib.ext, sep = ""))
ff_ode_desolve_c <- function(model_switch, xt, parameters, poped.db){
  with(as.list(parameters),{
    A_ini <- c(A1=0, A2=0)
    
    #Set up time points for the ODE
    times_xt <- drop(xt)
    times <- c(0,times_xt) ## add extra time for start of integration
    dose_times = seq(from=0,to=max(times_xt),by=TAU)
    times <- c(times,dose_times)
    times <- sort(times) 
    times <- unique(times) # remove duplicates
    
    eventdat <- data.frame(var = c("A1"), 
                           time = dose_times,
                           value = c(DOSE), method = c("add"))
    
    out <- ode(A_ini, times, func = "derivs", 
               parms = c(CL,V,KA), 
               dllname = "one_comp_oral_CL",
               initfunc = "initmod", nout = 1, 
               outnames = "Sum",
               events = list(data = eventdat),
               atol=1e-8, rtol=1e-8,maxsteps=5000)

    # grab timepoint values
    out = out[match(times_xt,out[,"time"]),]
    
    y = out[, "A2"]/V
    
    y=cbind(y) # must be a column matrix 
    return(list(y=y,poped.db=poped.db))
  })
}

ODE solution using deSolve and compiled C++ code (via Rcpp)

Here we define the ODE system using inline C++ code via Rcpp

cppFunction('List one_comp_oral_rcpp(double Time, NumericVector A, NumericVector Pars) {
int n = A.size();
NumericVector dA(n);

double CL = Pars[0];
double V = Pars[1];
double KA = Pars[2];

dA[0] = -KA*A[0];
dA[1] = KA*A[0] - (CL/V)*A[1];
return List::create(dA);
}')
ff_ode_desolve_rcpp <- function(model_switch, xt, p, poped.db){
    A_ini <- c(A1=0, A2=0)
    
    #Set up time points for the ODE
    times_xt <- drop(xt)
    times <- c(0,times_xt) ## add extra time for start of integration
    dose_times = seq(from=0,to=max(times_xt),by=p[["TAU"]])
    times <- c(times,dose_times)
    times <- sort(times) 
    times <- unique(times) # remove duplicates
    
    eventdat <- data.frame(var = c("A1"), 
                           time = dose_times,
                           value = c(p[["DOSE"]]), method = c("add"))
    
    out <- ode(A_ini, times, 
               one_comp_oral_rcpp, 
               c(CL=p[["CL"]],V=p[["V"]], KA=p[["KA"]]), 
               events = list(data = eventdat),
               atol=1e-8, rtol=1e-8,maxsteps=5000)
    
    
    # grab timepoint values for central comp
    y = out[match(times_xt,out[,"time"]),"A2",drop=F]/p[["V"]]
    
    return(list(y=y,poped.db=poped.db))
}

ODE solution using PKPDsim

The same model written as a set of ODEs using PKPDsim:

pk1cmtoral <- PKPDsim::new_ode_model("pk_1cmt_oral") # take from library
ff_ode_pkpdsim <- function(model_switch, xt, p, poped.db){
    #Set up time points for the ODE
    times_xt <- drop(xt)  
    dose_times <- seq(from=0,to=max(times_xt),by=p[["TAU"]])
    times <- sort(unique(c(0,times_xt,dose_times)))

    N = length(dose_times)
    regimen = PKPDsim::new_regimen(amt=p[["DOSE"]],n=N,interval=p[["TAU"]])
    design <- PKPDsim::sim(
      ode = pk1cmtoral, 
      parameters = c(CL=p[["CL"]],V=p[["V"]],KA=p[["KA"]]), 
      regimen = regimen,
      only_obs = TRUE,
      t_obs = times,
      checks = FALSE,
      return_design = TRUE)
    tmp <- PKPDsim::sim_core(sim_object = design, ode = pk1cmtoral)
    y <- tmp$y
    m_tmp <- match(round(times_xt,digits = 6),tmp[,"t"])
    if(any(is.na(m_tmp))){
      stop("can't find time points in solution\n", 
           "try changing the digits argument in the match function")
    } 
    
    y <- y[m_tmp]
    return(list(y = y, poped.db = poped.db))
}

ODE solution using mrgsolve

The same model written as a set of ODEs using mrgsolve:

code <- '
$PARAM CL=3.75, V=72.8, KA=0.25
$CMT DEPOT CENT
$ODE
dxdt_DEPOT = -KA*DEPOT;
dxdt_CENT = KA*DEPOT - (CL/V)*CENT;
$TABLE double CP  = CENT/V;
$CAPTURE CP
'

Compile and load the model with mcode

moda <- mcode("optim", code, atol=1e-8, rtol=1e-8,maxsteps=5000)
#> Building optim ... done.
ff_ode_mrg <- function(model_switch, xt, p, poped.db){
  times_xt <- drop(xt)  
  dose_times <- seq(from=0,to=max(times_xt),by=p[["TAU"]])
  time <- sort(unique(c(0,times_xt,dose_times)))
  is.dose <- time %in% dose_times
  
  data <- 
    tibble::tibble(ID = 1,
                      time = time,
                      amt = ifelse(is.dose,p[["DOSE"]], 0), 
                      cmt = ifelse(is.dose, 1, 0), 
                      evid = cmt,
                      CL = p[["CL"]], V = p[["V"]], KA = p[["KA"]])
  
  out <- mrgsim_q(moda, data=data)
  
  y <-  out$CP
  
  y <- y[match(times_xt,out$time)]
  
  return(list(y=matrix(y,ncol=1),poped.db=poped.db))
  
}

ODE solution using RxODE

The model written for RxODE:

modrx <- RxODE({
  d/dt(DEPOT) = -KA*DEPOT;
  d/dt(CENT) = KA*DEPOT - (CL/V)*CENT;
  CP=CENT/V;
})
ff_ode_rx <- function(model_switch, xt, p, poped.db){
  times_xt <- drop(xt)
  et(0,amt=p[["DOSE"]], ii=p[["TAU"]], until=max(times_xt)) %>%
    et(times_xt) -> data
  
  out <- rxSolve(modrx, p, data, atol=1e-8, rtol=1e-8,maxsteps=5000,
                 returnType="data.frame")
  
  y <-  out$CP[match(times_xt,out$time)]
  
  return(list(y=matrix(y,ncol=1),poped.db=poped.db))
  
}

Common model elements

Other functions are used to define BSV and RUV.


sfg <- function(x,a,bpop,b,bocc){
  parameters=c( 
    KA=bpop[1]*exp(b[1]),
    CL=bpop[2]*exp(b[2]),
    V=bpop[3]*exp(b[3]),
    DOSE=a[1],
    TAU=a[2])
  return( parameters ) 
}

feps <- function(model_switch,xt,parameters,epsi,poped.db){
  y <- do.call(poped.db$model$ff_pointer,list(model_switch,xt,parameters,poped.db))[[1]]
  y = y*(1+epsi[,1])+epsi[,2]
  return(list(y=y,poped.db=poped.db)) 
}

Create PopED databases

Next we define the model to use, the parameters of those models, the intial design design and design space for any design calculation. Here we create a number of databases that correspond to different model implementations.

The initial design is a 2 group design, with doses of 20 mg or 40 mg every 24 hours. Each group has the same sampling schedule, with 3 samples in the first day of the study and 2 on the 10th day of the study.

poped_db_analytic <- create.poped.database(ff_fun =ff,
                                  fg_fun =sfg,
                                  fError_fun=feps,
                                  bpop=c(KA=0.25,CL=3.75,V=72.8), 
                                  d=c(KA=0.09,CL=0.25^2,V=0.09), 
                                  sigma=c(prop=0.04,add=0.0025),
                                  m=2,
                                  groupsize=20,
                                  xt=c( 1,2,8,240,245),
                                  minxt=c(0,0,0,240,240),
                                  maxxt=c(10,10,10,248,248),
                                  bUseGrouped_xt=1,
                                  a=cbind(DOSE=c(20,40),TAU=c(24,24)),
                                  maxa=c(DOSE=200,TAU=24),
                                  mina=c(DOSE=0,TAU=24))


poped_db_ode_desolve <- create.poped.database(poped_db_analytic,ff_fun = ff_ode_desolve)
poped_db_ode_desolve_c <- create.poped.database(poped_db_analytic,ff_fun = ff_ode_desolve_c)
poped_db_ode_desolve_rcpp <- create.poped.database(poped_db_analytic,ff_fun = ff_ode_desolve_rcpp)
poped_db_ode_pkpdsim <- create.poped.database(poped_db_analytic,ff_fun = ff_ode_pkpdsim)
poped_db_ode_mrg <- create.poped.database(poped_db_analytic,ff_fun = ff_ode_mrg)
poped_db_ode_rx <- create.poped.database(poped_db_analytic,ff_fun = ff_ode_rx)

Model predictions

So are there difference in the model predictions between the different implementations?

Here is a visual representation of the model predictions for this study design, based on the analytic solution:

plot_model_prediction(poped_db_analytic,model_num_points = 500,PI=T,separate.groups = T) 

We can compare the different predictions in this plot accross model implementations. Here we see that the accuracy of the different methods are within machine precision (or very small).

pred_std <- model_prediction(poped_db_analytic,model_num_points = 500,include_sample_times = TRUE,PI = TRUE)

pred_ode_desolve <- model_prediction(poped_db_ode_desolve,
                                     model_num_points = 500,
                                     include_sample_times = TRUE,
                                     PI = TRUE)
all.equal(pred_std,pred_ode_desolve)
#> [1] TRUE

pred_ode_desolve_c <- model_prediction(poped_db_ode_desolve_c,
                                       model_num_points = 500,
                                       include_sample_times = TRUE,
                                       PI = TRUE)
all.equal(pred_std,pred_ode_desolve_c)
#> [1] TRUE

pred_ode_desolve_rcpp <- model_prediction(poped_db_ode_desolve_rcpp,
                                          model_num_points = 500,
                                          include_sample_times = TRUE,
                                          PI = TRUE)
all.equal(pred_std,pred_ode_desolve_rcpp)
#> [1] TRUE

pred_ode_pkpdsim <- model_prediction(poped_db_ode_pkpdsim,
                                     model_num_points = 500,
                                     include_sample_times = TRUE,
                                     PI = TRUE)
all.equal(pred_std,pred_ode_pkpdsim)
#> [1] "Component \"PI_l\": Mean relative difference: 3.593837e-08"

pred_ode_mrg <- model_prediction(poped_db_ode_mrg,
                                 model_num_points = 500,
                                 include_sample_times = TRUE,
                                 PI = TRUE)
all.equal(pred_std,pred_ode_mrg)
#> [1] TRUE

pred_ode_rx <- model_prediction(poped_db_ode_rx,
                                 model_num_points = 500,
                                 include_sample_times = TRUE,
                                 PI = TRUE)
all.equal(pred_std,pred_ode_rx)
#> [1] TRUE

Evaluate the design

Here we compare the computation of the Fisher Information Matrix (FIM). By comparing the \(ln(det(FIM))\) (the lnD-objective function value, or ofv).

(eval_std <- evaluate_design(poped_db_analytic))
#> $ofv
#> [1] 48.98804
#> 
#> $fim
#>                   KA           CL           V       d_KA        d_CL        d_V
#> KA       1695.742314 -11.73537527 -6.75450789    0.00000     0.00000    0.00000
#> CL        -11.735375  29.99735715 -0.03288331    0.00000     0.00000    0.00000
#> V          -6.754508  -0.03288331  0.04213359    0.00000     0.00000    0.00000
#> d_KA        0.000000   0.00000000  0.00000000  147.24270     1.52226  192.23403
#> d_CL        0.000000   0.00000000  0.00000000    1.52226  2254.55188    1.21987
#> d_V         0.000000   0.00000000  0.00000000  192.23403     1.21987  634.42055
#> sig_prop    0.000000   0.00000000  0.00000000  148.86724   844.57325  387.53816
#> sig_add     0.000000   0.00000000  0.00000000 6555.68433 14391.88132 8669.58391
#>             sig_prop     sig_add
#> KA            0.0000       0.000
#> CL            0.0000       0.000
#> V             0.0000       0.000
#> d_KA        148.8672    6555.684
#> d_CL        844.5733   14391.881
#> d_V         387.5382    8669.584
#> sig_prop   7759.5374  110702.705
#> sig_add  110702.7045 4436323.946
#> 
#> $rse
#>         KA         CL          V       d_KA       d_CL        d_V   sig_prop 
#>  16.285678   4.909749  11.209270 120.825798  34.448477  57.300408  36.104027 
#>    sig_add 
#>  24.339781

All the computations give very similar results:

eval_ode_desolve <- evaluate_design(poped_db_ode_desolve) 
all.equal(eval_std$ofv,eval_ode_desolve$ofv)
#> [1] "Mean relative difference: 2.493735e-08"

eval_ode_desolve_c <- evaluate_design(poped_db_ode_desolve_c) 
all.equal(eval_std$ofv,eval_ode_desolve_c$ofv)
#> [1] "Mean relative difference: 2.493735e-08"

eval_ode_desolve_rccp <- evaluate_design(poped_db_ode_desolve_rcpp) 
all.equal(eval_std$ofv,eval_ode_desolve_rccp$ofv)
#> [1] "Mean relative difference: 2.493735e-08"

eval_ode_pkpdsim <- evaluate_design(poped_db_ode_pkpdsim) 
all.equal(eval_std$ofv,eval_ode_pkpdsim$ofv)
#> [1] TRUE

eval_ode_mrg <- evaluate_design(poped_db_ode_mrg) 
all.equal(eval_std$ofv,eval_ode_mrg$ofv)
#> [1] "Mean relative difference: 2.361616e-08"

Speed of FIM computation

We can compare the speed of the computations. Analytic solutions are fast, as expected, in this case more than 20 times faster than any of the ODE methods. mrgsolve is the fastest of the ODE solvers in this example. Note that, in previous work (http://pkpdsim.ronkeizer.com/speed.html), much of the speed difference between mrgsolve, RxODE and PKPDsim has been found to be due to the overhead from pre- and post-processing of the simulation from ODE systems. Other ways of handling the pre- and post-processing may speed up these computations.

library(microbenchmark)
library(ggplot2)

compare <- microbenchmark(
  evaluate_design(poped_db_analytic),
  evaluate_design(poped_db_ode_desolve),
  evaluate_design(poped_db_ode_desolve_c),
  evaluate_design(poped_db_ode_desolve_rcpp),
  evaluate_design(poped_db_ode_pkpdsim),
  evaluate_design(poped_db_ode_mrg),
  evaluate_design(poped_db_ode_rx),
  times = 100L)

autoplot(compare)

Version information

devtools::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#>  setting  value                       
#>  version  R version 4.1.0 (2021-05-18)
#>  os       macOS Big Sur 10.16         
#>  system   x86_64, darwin17.0          
#>  ui       X11                         
#>  language (EN)                        
#>  collate  en_US.UTF-8                 
#>  ctype    en_US.UTF-8                 
#>  tz       Europe/Stockholm            
#>  date     2021-05-21                  
#> 
#> ─ Packages ───────────────────────────────────────────────────────────────────
#>  ! package       * version  date       lib source                            
#>    assertthat      0.2.1    2019-03-21 [2] CRAN (R 4.1.0)                    
#>    backports       1.2.1    2020-12-09 [2] CRAN (R 4.1.0)                    
#>    BH              1.75.0-0 2021-01-11 [2] CRAN (R 4.1.0)                    
#>    bookdown        0.22     2021-04-22 [2] CRAN (R 4.1.0)                    
#>    boot            1.3-28   2021-05-03 [2] CRAN (R 4.1.0)                    
#>    bslib           0.2.5.1  2021-05-18 [2] CRAN (R 4.1.0)                    
#>    cachem          1.0.5    2021-05-15 [2] CRAN (R 4.1.0)                    
#>    callr           3.7.0    2021-04-20 [2] CRAN (R 4.1.0)                    
#>    checkmate       2.0.0    2020-02-06 [2] CRAN (R 4.1.0)                    
#>    cli             2.5.0    2021-04-26 [2] CRAN (R 4.1.0)                    
#>    codetools       0.2-18   2020-11-04 [2] CRAN (R 4.1.0)                    
#>    colorspace      2.0-1    2021-05-04 [2] CRAN (R 4.1.0)                    
#>    crayon          1.4.1    2021-02-08 [2] CRAN (R 4.1.0)                    
#>    data.table      1.14.0   2021-02-21 [2] CRAN (R 4.1.0)                    
#>    DBI             1.1.1    2021-01-15 [2] CRAN (R 4.1.0)                    
#>    desc            1.3.0    2021-03-05 [2] CRAN (R 4.1.0)                    
#>    deSolve       * 1.28     2020-03-08 [2] CRAN (R 4.1.0)                    
#>    devtools        2.4.1    2021-05-05 [2] CRAN (R 4.1.0)                    
#>    digest          0.6.27   2020-10-24 [2] CRAN (R 4.1.0)                    
#>    dparser         1.3.1-4  2021-04-07 [2] CRAN (R 4.1.0)                    
#>    dplyr           1.0.6    2021-05-05 [2] CRAN (R 4.1.0)                    
#>    ellipsis        0.3.2    2021-04-29 [2] CRAN (R 4.1.0)                    
#>    evaluate        0.14     2019-05-28 [2] CRAN (R 4.1.0)                    
#>    fansi           0.4.2    2021-01-15 [2] CRAN (R 4.1.0)                    
#>    farver          2.1.0    2021-02-28 [2] CRAN (R 4.1.0)                    
#>    fastmap         1.1.0    2021-01-25 [2] CRAN (R 4.1.0)                    
#>    fs              1.5.0    2020-07-31 [2] CRAN (R 4.1.0)                    
#>    generics        0.1.0    2020-10-31 [2] CRAN (R 4.1.0)                    
#>    ggplot2         3.3.3    2020-12-30 [2] CRAN (R 4.1.0)                    
#>    glue            1.4.2    2020-08-27 [2] CRAN (R 4.1.0)                    
#>    gtable          0.3.0    2019-03-25 [2] CRAN (R 4.1.0)                    
#>    gtools          3.8.2    2020-03-31 [2] CRAN (R 4.1.0)                    
#>    highr           0.9      2021-04-16 [2] CRAN (R 4.1.0)                    
#>    htmltools       0.5.1.1  2021-01-22 [2] CRAN (R 4.1.0)                    
#>    httr            1.4.2    2020-07-20 [2] CRAN (R 4.1.0)                    
#>    jquerylib       0.1.4    2021-04-26 [2] CRAN (R 4.1.0)                    
#>    jsonlite        1.7.2    2020-12-09 [2] CRAN (R 4.1.0)                    
#>    kableExtra    * 1.3.4    2021-02-20 [2] CRAN (R 4.1.0)                    
#>    knitr         * 1.33     2021-04-24 [2] CRAN (R 4.1.0)                    
#>    labeling        0.4.2    2020-10-20 [2] CRAN (R 4.1.0)                    
#>    lifecycle       1.0.0    2021-02-15 [2] CRAN (R 4.1.0)                    
#>    lotri           0.3.1    2021-01-05 [2] CRAN (R 4.1.0)                    
#>    lubridate       1.7.10   2021-02-26 [2] CRAN (R 4.1.0)                    
#>    magrittr        2.0.1    2020-11-17 [2] CRAN (R 4.1.0)                    
#>    MASS            7.3-54   2021-05-03 [2] CRAN (R 4.1.0)                    
#>    memoise         2.0.0    2021-01-26 [2] CRAN (R 4.1.0)                    
#>    mrgsolve      * 0.11.1   2021-05-10 [2] CRAN (R 4.1.0)                    
#>    munsell         0.5.0    2018-06-12 [2] CRAN (R 4.1.0)                    
#>    mvtnorm         1.1-1    2020-06-09 [2] CRAN (R 4.1.0)                    
#>    pillar          1.6.1    2021-05-16 [2] CRAN (R 4.1.0)                    
#>    pkgbuild        1.2.0    2020-12-15 [2] CRAN (R 4.1.0)                    
#>    pkgconfig       2.0.3    2019-09-22 [2] CRAN (R 4.1.0)                    
#>    pkgdown         1.6.1    2020-09-12 [2] CRAN (R 4.1.0)                    
#>    pkgload         1.2.1    2021-04-06 [2] CRAN (R 4.1.0)                    
#>    PKPDsim       * 1.0.7    2021-05-21 [2] Github (InsightRX/PKPDsim@bd9b9e4)
#>  P PopED         * 0.6.0    2021-05-21 [?] local                             
#>    PreciseSums     0.4      2020-07-10 [2] CRAN (R 4.1.0)                    
#>    prettyunits     1.1.1    2020-01-24 [2] CRAN (R 4.1.0)                    
#>    processx        3.5.2    2021-04-30 [2] CRAN (R 4.1.0)                    
#>    ps              1.6.0    2021-02-28 [2] CRAN (R 4.1.0)                    
#>    purrr           0.3.4    2020-04-17 [2] CRAN (R 4.1.0)                    
#>    qs              0.24.1   2021-03-07 [2] CRAN (R 4.1.0)                    
#>    R6              2.5.0    2020-10-28 [2] CRAN (R 4.1.0)                    
#>    ragg            1.1.2    2021-03-17 [2] CRAN (R 4.1.0)                    
#>    RApiSerialize   0.1.0    2014-04-19 [2] CRAN (R 4.1.0)                    
#>    Rcpp          * 1.0.6    2021-01-15 [2] CRAN (R 4.1.0)                    
#>    RcppParallel    5.1.4    2021-05-04 [2] CRAN (R 4.1.0)                    
#>    remotes         2.3.0    2021-04-01 [2] CRAN (R 4.1.0)                    
#>    rlang           0.4.11   2021-04-30 [2] CRAN (R 4.1.0)                    
#>    rmarkdown       2.8      2021-05-07 [2] CRAN (R 4.1.0)                    
#>    rprojroot       2.0.2    2020-11-15 [2] CRAN (R 4.1.0)                    
#>    rstudioapi      0.13     2020-11-12 [2] CRAN (R 4.1.0)                    
#>    rvest           1.0.0    2021-03-09 [2] CRAN (R 4.1.0)                    
#>    RxODE         * 1.0.9    2021-05-11 [2] CRAN (R 4.1.0)                    
#>    sass            0.4.0    2021-05-12 [2] CRAN (R 4.1.0)                    
#>    scales          1.1.1    2020-05-11 [2] CRAN (R 4.1.0)                    
#>    sessioninfo     1.1.1    2018-11-05 [2] CRAN (R 4.1.0)                    
#>    stringfish      0.15.1   2021-03-16 [2] CRAN (R 4.1.0)                    
#>    stringi         1.6.2    2021-05-17 [2] CRAN (R 4.1.0)                    
#>    stringr         1.4.0    2019-02-10 [2] CRAN (R 4.1.0)                    
#>    svglite         2.0.0    2021-02-20 [2] CRAN (R 4.1.0)                    
#>    sys             3.4      2020-07-23 [2] CRAN (R 4.1.0)                    
#>    systemfonts     1.0.2    2021-05-11 [2] CRAN (R 4.1.0)                    
#>    testthat      * 3.0.2    2021-02-14 [2] CRAN (R 4.1.0)                    
#>    textshaping     0.3.4    2021-05-11 [2] CRAN (R 4.1.0)                    
#>    tibble          3.1.2    2021-05-16 [2] CRAN (R 4.1.0)                    
#>    tidyselect      1.1.1    2021-04-30 [2] CRAN (R 4.1.0)                    
#>    usethis         2.0.1    2021-02-10 [2] CRAN (R 4.1.0)                    
#>    utf8            1.2.1    2021-03-12 [2] CRAN (R 4.1.0)                    
#>    vctrs           0.3.8    2021-04-29 [2] CRAN (R 4.1.0)                    
#>    viridisLite     0.4.0    2021-04-13 [2] CRAN (R 4.1.0)                    
#>    webshot         0.5.2    2019-11-22 [2] CRAN (R 4.1.0)                    
#>    withr           2.4.2    2021-04-18 [2] CRAN (R 4.1.0)                    
#>    xfun            0.23     2021-05-15 [2] CRAN (R 4.1.0)                    
#>    xml2            1.3.2    2020-04-23 [2] CRAN (R 4.1.0)                    
#>    yaml            2.2.1    2020-02-01 [2] CRAN (R 4.1.0)                    
#> 
#> [1] /private/var/folders/g6/mk9zc6552f7g93c457xzgvfc0000gp/T/Rtmp2dRrIZ/temp_libpath7ebb50981b8f
#> [2] /Library/Frameworks/R.framework/Versions/4.1/Resources/library
#> 
#>  P ── Loaded and on-disk path mismatch.
#sessionInfo()