Juvena PK study pipeline
Bill Chen
2/15/2021
PK study
In order to utilize most of the PK library, we need to have a three steps for get best of the results.
- Using linear regression script to find initial values for model search
- Using NlmixR to solve the results.
- Uisng mrgsolve to simuulate
1. Linear regression script
library(readr)
library(tidyverse)## -- Attaching packages --------------------------------------- tidyverse 1.3.0 --## v ggplot2 3.3.3 v dplyr 1.0.4
## v tibble 3.0.6 v stringr 1.4.0
## v tidyr 1.1.2 v forcats 0.5.0
## v purrr 0.3.4## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(ggplot2)
library(nlmixr)
pk_dat <- read_csv("data/011921_nlmixr.csv")##
## -- Column specification --------------------------------------------------------
## cols(
## ID = col_double(),
## TIME = col_double(),
## DV = col_double(),
## AMT = col_double(),
## EVID = col_double(),
## CMT = col_double(),
## WT = col_double()
## )pk_dat## # A tibble: 32 x 7
## ID TIME DV AMT EVID CMT WT
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 0 4.96 500 1 2 22.3
## 2 1 3 11.2 0 0 0 22.3
## 3 1 10 11.5 0 0 0 22.3
## 4 1 15 9.67 0 0 0 22.3
## 5 1 30 8.51 0 0 0 22.3
## 6 1 60 10.4 0 0 0 22.3
## 7 1 1440 2.71 0 0 0 22.3
## 8 1 2880 0.642 0 0 0 22.3
## 9 2 0 0.436 500 1 2 25.4
## 10 2 3 9.32 0 0 0 25.4
## # ... with 22 more rowsHere time is in minutes, most of the PK study using hours, convert to hour
pk_dat$TIME <- pk_dat$TIME / 60
pk_dat## # A tibble: 32 x 7
## ID TIME DV AMT EVID CMT WT
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 0 4.96 500 1 2 22.3
## 2 1 0.05 11.2 0 0 0 22.3
## 3 1 0.167 11.5 0 0 0 22.3
## 4 1 0.25 9.67 0 0 0 22.3
## 5 1 0.5 8.51 0 0 0 22.3
## 6 1 1 10.4 0 0 0 22.3
## 7 1 24 2.71 0 0 0 22.3
## 8 1 48 0.642 0 0 0 22.3
## 9 2 0 0.436 500 1 2 25.4
## 10 2 0.05 9.32 0 0 0 25.4
## # ... with 22 more rowsNow using the pre-define model to find the initial value: - One compartment 
Two compartment
Three compartment
Using pk script to find the initial value
source("previous_implementation.R")
fitting_func <- function(dat_to_fit, num_of_compartments=1, approach='bolus') {
formula_df <- formula_func()
experiment <- formula_df[formula_df[["num_of_compartments"]]==num_of_compartments &
formula_df[["delivery_approach"]]==approach,]
nls_formula <- as.formula(experiment[['formula_str']])
paras <- init_para(dat_to_fit, num_of_compartments)
nls_out <- NULL
predicted_t <- seq(from = dat_to_fit$t[1], to = dat_to_fit$t[nrow(dat_to_fit)], length.out = 100)
predicted_c <- NULL
if (num_of_compartments==1 & approach=='bolus'){
start_point <- c(paras[['init_para']][['init_val']])
start_point <- list(a = start_point[1], b = start_point[2])
}
else if (num_of_compartments==2 & approach=='bolus') {
start_point <- c(paras[['init_para']][['init_val']])
start_point <- list(a = start_point[1],
b = start_point[2],
c = start_point[3],
d = start_point[4])
}
else if (num_of_compartments==3 & approach=='bolus') {
start_point <- c(paras[['init_para']][['init_val']])
start_point <- list(a = start_point[1],
b = start_point[2],
c = start_point[3],
d = start_point[4],
e = start_point[5],
g = start_point[6])
}
nls_out <- try(nls(nls_formula, data = dat_to_fit, lower = paras[['init_para']][['min_val']], upper = paras[['init_para']][['max_val']],
start = start_point, algorithm = "default"), silent = TRUE)
if (class(nls_out)=="try-error") {
fn <- nls_bolus_error_solution(dat_to_fit, num_of_compartments, approach)
# print(fn)
nls_out <- try(nloptr::bobyqa(as.vector(unlist(start_point)), fn, lower = paras[['init_para']][['min_val']], upper = paras[['init_para']][['max_val']]),
silent = TRUE)
}
return(nls_out)
}dat_init <- pk_dat[c("TIME", "DV", "AMT", "ID")]
colnames(dat_init) <- c('t',"con",'dose',"ID")
temp_output <- fitting_func(dat_init)## Warning in nls(nls_formula, data = dat_to_fit, lower = paras[["init_para"]]
## [["min_val"]], : upper and lower bounds ignored unless algorithm = "port"print(temp_output)## Nonlinear regression model
## model: con ~ a * exp(-b * t)
## data: dat_to_fit
## a b
## 11.66041 0.04068
## residual sum-of-squares: 937
##
## Number of iterations to convergence: 5
## Achieved convergence tolerance: 1.58e-06temp_output <- coef(temp_output)2. Using NlmixR to solve the results.
cl / v = ke = k_10
pk_dat$EVID<- ifelse(pk_dat$EVID==1,101,pk_dat$EVID)
one.cmt <- function() {
ini({
tka <- 0.45 # Log Ka
tcl <- 1 # Log Cl
tv <- 3.45 # Log V
eta.ka ~ 0.6
eta.cl ~ 0.3
eta.v ~ 0.1
add.sd <- 0.7
})
model({
ka <- exp(tka + eta.ka)
cl <- exp(tcl + eta.cl)
v <- exp(tv + eta.v)
d/dt(depot) = -ka * depot
d/dt(center) = ka * depot - cl / v * center
cp = center / v
cp ~ add(add.sd)
})
}
f <- nlmixr(one.cmt)## i parameter labels from comments will be replaced by 'label()'fit2 <- nlmixr(one.cmt, pk_dat, list(print=0), est="saem")## i parameter labels from comments will be replaced by 'label()'## > generate SAEM model## v done## qs v0.23.5.## RxODE 1.0.3 using 2 threads (see ?getRxThreads)
## no cache: create with `rxCreateCache()`## Calculating covariance matrix## Calculating residuals/tables## doneprint(fit2)## -- nlmixr SAEM(ODE); OBJF not calculated fit -----------------------------------
##
## Gaussian/Laplacian Likelihoods: AIC() or $objf etc.
## FOCEi CWRES & Likelihoods: addCwres()
##
## -- Time (sec $time): -----------------------------------------------------------
##
## saem setup table covariance other
## elapsed 5.99 1.576 0.03 0.01 0.584
##
## -- Population Parameters ($parFixed or $parFixedDf): ---------------------------
##
## Parameter Est. SE %RSE Back-transformed(95%CI) BSV(CV%)
## tka Log Ka 6.2 14.9 240 495 (1.02e-010, 2.4e+015) 0.772
## tcl Log Cl 0.555 0.106 19.1 1.74 (1.41, 2.14) 0.374
## tv Log V 3.63 0.154 4.23 37.8 (28, 51.1) 30.9
## add.sd 1.5 1.5
## Shrink(SD)%
## tka 91.5%
## tcl 96.1%
## tv -13.8%
## add.sd
##
## Covariance Type ($covMethod): linFim
## Fixed parameter correlations in $cor
## No correlations in between subject variability (BSV) matrix
## Full BSV covariance ($omega) or correlation ($omegaR; diagonals=SDs)
## Distribution stats (mean/skewness/kurtosis/p-value) available in $shrink
##
## -- Fit Data (object is a modified tibble): -------------------------------------
## # A tibble: 28 x 19
## ID TIME DV PRED RES IPRED IRES IWRES eta.ka eta.cl eta.v cp
## <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 0.05 11.2 13.2 -2.01 10.5 0.713 0.477 5.36e-4 1.36e-6 0.231 10.5
## 2 1 0.167 11.5 13.1 -1.58 10.4 1.12 0.746 5.36e-4 1.36e-6 0.231 10.4
## 3 1 0.25 9.67 13.1 -3.40 10.4 -0.721 -0.482 5.36e-4 1.36e-6 0.231 10.4
## # ... with 25 more rows, and 7 more variables: depot <dbl>, center <dbl>,
## # ka <dbl>, cl <dbl>, v <dbl>, tad <dbl>, dosenum <dbl>plot(fit2)
iter <- fit2$par.hist.stacked
iter$Parameter[iter$par=="add.sd"] <- "Additive error"
iter$Parameter[iter$par=="eta.cl"] <- "IIV CL/F"
iter$Parameter[iter$par=="eta.v"] <- "IIV V/F"
iter$Parameter[iter$par=="eta.ka"] <- "IIV ka"
iter$Parameter[iter$par=="tcl"] <- "log(CL/F)"
iter$Parameter[iter$par=="tv"] <- "log(V/F)"
iter$Parameter[iter$par=="tka"] <- "log(ka)"
iter$Parameter <- ordered(iter$Parameter, c("log(CL/F)", "log(V/F)", "log(ka)",
"IIV CL/F", "IIV V/F", "IIV ka",
"Additive error"))
ggplot(iter, aes(iter, val)) +
geom_line(col="red") +
scale_x_continuous("Iteration") +
scale_y_continuous("Value") +
facet_wrap(~ Parameter, scales="free_y") +
labs(title="Theophylline single-dose", subtitle="Parameter estimation iterations")
Diagnosis
library(xpose)##
## Attaching package: 'xpose'## The following object is masked from 'package:nlmixr':
##
## vpc## The following object is masked from 'package:stats':
##
## filterlibrary(xpose.nlmixr)##
## Attaching package: 'xpose.nlmixr'## The following object is masked from 'package:xpose':
##
## vpcxp <- xpose_data_nlmixr(fit2)## Calculating residuals/tables## done## Warning in xpose_data_nlmixr(fit2): Added CWRES to fit (using nlmixr::addCwres)## Warning: `as.tibble()` is deprecated as of tibble 2.0.0.
## Please use `as_tibble()` instead.
## The signature and semantics have changed, see `?as_tibble`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.## Warning: `arrange_()` is deprecated as of dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.dv_vs_pred(xp)## `geom_smooth()` using formula 'y ~ x'## `geom_smooth()` using formula 'y ~ x'
dv_vs_ipred(xp)## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
absval_res_vs_pred(xp, res="IWRES")## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'