Specifying a PK Model using mrgsolve
codepkmodelcov <- '
$PARAM @annotated
KA : 0.5 : Absorption rate constant Ka (1/h)
CL : 4 : Clearance CL (L/h)
V : 10 : Central volume Vc (L)
Vp : 50 : Peripheral volume Vp (L)
Qp : 10 : Intercompartmental clearance Q (L/h)
CLALB : -0.8 : Ablumin on CL (ref. 45 g/L)
CLSEX : 0.2 : Sex on CL (ref. Female)
CLWT : 1 : Weight on CL (ref. 85 kg)
VSEX : 0.07 : Sex on Vc (ref. Female)
VWT : 1 : Weight on Vc (ref. 85 kg)
$PARAM @annotated // reference values for covariate
WT : 85 : Weight (kg)
SEX : 0 : Sex (0=Female, 1=Male)
ALB : 45 : Albumin (g/L)
$PKMODEL cmt="GUT CENT PER", depot=TRUE, trans=11
$MAIN
double CLi = CL *
pow((ALB/45.0), CLALB)*
(SEX == 1.0 ? (1.0+CLSEX) : 1.0)*
pow((WT/85.0), CLWT)*exp(nCL);
double V2i = V *
(SEX == 1.0 ? (1.0+VSEX) : 1.0)*
pow((WT/85.0), VWT)*exp(nVC);
double KAi = KA;
double V3i = Vp *pow((WT/85.0), 1);
double Qi = Qp *pow((WT/85.0), 0.75);
$OMEGA @annotated @block
nCL : 0.09 : ETA on CL
nVC : 0.01 0.09 : ETA on Vc
$TABLE
double CP = CENT/V2i;
$CAPTURE CP KAi CLi V2i V3i Qi WT SEX ALB
'
modcovsim <- mcode("codepkmodelcov", codepkmodelcov)
partab <- setDT(modcovsim@annot$data)[block=="PARAM", .(name, descr, unit)]
partab <- merge(partab, melt(setDT(modcovsim@param@data), meas=patterns("*"), var="name"))
knitr::kable(partab)| name | descr | unit | value |
|---|---|---|---|
| ALB | Albumin | g/L | 45.00 |
| CL | Clearance CL | L/h | 4.00 |
| CLALB | Ablumin on CL | ref. 45 g/L | -0.80 |
| CLSEX | Sex on CL | ref. Female | 0.20 |
| CLWT | Weight on CL | ref. 85 kg | 1.00 |
| KA | Absorption rate constant Ka | 1/h | 0.50 |
| Qp | Intercompartmental clearance Q | L/h | 10.00 |
| SEX | Sex | 0=Female, 1=Male | 0.00 |
| V | Central volume Vc | L | 10.00 |
| VSEX | Sex on Vc | ref. Female | 0.07 |
| VWT | Weight on Vc | ref. 85 kg | 1.00 |
| Vp | Peripheral volume Vp | L | 50.00 |
| WT | Weight | kg | 85.00 |
Simulate a Reference Subjects with BSV
We simulate the reference subject having the reference covariate values defined in the model which are:
Weight = 85 kg, Sex = Female and Albumin = 45 g/L. We also keep the between subject variability (BSV) to illustrate its effects on the concentration-time profiles on linear and log linear scales.
idata <- data.table(ID=1:nbsvsubjects, WT=85, SEX=0, ALB=45)
ev1 <- ev(time = 0, amt = 100, cmt = 1)
data.dose <- ev(ev1)
data.dose <- setDT(as.data.frame(data.dose))
data.all <- data.table(idata, data.dose)
outputsim <- modcovsim %>%
data_set(data.all) %>%
mrgsim(end = 24, delta = 0.25) %>%
as.data.frame %>%
as.data.table
outputsim$SEX <- factor(outputsim$SEX, labels="Female")
# Only plot a random sample of N=500
set.seed(678549)
plotdata <- outputsim[ID %in% sample(unique(ID), 500)]
p1 <- ggplot(plotdata, aes(time, CP, group = ID)) +
geom_line(alpha = 0.2, size = 0.1) +
facet_grid(~ WT + SEX + ALB, labeller = label_both) +
labs(y = "Plasma Concentrations", x = "Time (h)")
p2 <- ggplot(plotdata, aes(time, CP, group = ID)) +
geom_line(alpha = 0.2, size = 0.1) +
facet_grid(~ WT + SEX + ALB, labeller = label_both) +
scale_y_log10() +
labs(y = expression(Log[10]~Plasma~Concentrations), x = "Time (h)")
egg::ggarrange(p1, p2, ncol = 2)
Compute PK Parameters, Plot and Summarize BSV
In this section we compute the PK parameters of interest, provide a plot of the parameters as well as of the standardized ones. We also summarize and report the BSV as ranges of 50 and 90% of patients for each PK parameter.
derive.exposure <- function(time, CP) {
n <- length(time)
x <- c(
Cmax = max(CP),
Clast = CP[n],
AUC = sum(diff(time) * (CP[-1] + CP[-n])) / 2
)
data.table(paramname=names(x), paramvalue=x)
}
refbsv <- outputsim[, derive.exposure(time, CP), by=.(ID, WT, SEX, ALB)]
p3 <- ggplot(refbsv, aes(
x = paramvalue,
y = paramname,
fill = factor(..quantile..),
height = ..ndensity..)) +
facet_wrap(~ paramname, scales="free", ncol=1) +
stat_density_ridges(
geom="density_ridges_gradient", calc_ecdf=TRUE,
quantile_lines=TRUE, rel_min_height=0.001, scale=0.9,
quantiles=c(0.05, 0.25, 0.5, 0.75, 0.95)) +
scale_fill_manual(
name = "Probability",
values = c("white", "#FF000050", "#FF0000A0", "#FF0000A0", "#FF000050", "white"),
labels = c("(0, 0.05]", "(0.05, 0.25]",
"(0.25, 0.5]", "(0.5, 0.75]",
"(0.75, 0.95]", "(0.95, 1]")) +
theme_bw() +
theme(
legend.position = "none",
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.title.y = element_blank()) +
labs(x="PK Parameters", y="") +
scale_x_log10() +
coord_cartesian(expand=FALSE)
# Obtain the standardized parameter value by dividing by the median.
refbsv[, stdparamvalue := paramvalue/median(paramvalue), by=paramname]
p4 <- ggplot(refbsv, aes(
x = stdparamvalue,
y = paramname,
fill = factor(..quantile..),
height = ..ndensity..)) +
facet_wrap(~ paramname, scales="free_y", ncol=1) +
stat_density_ridges(
geom="density_ridges_gradient", calc_ecdf=TRUE,
quantile_lines=TRUE, rel_min_height=0.001, scale=0.9,
quantiles=c(0.05, 0.25, 0.5, 0.75, 0.95)) +
scale_fill_manual(
name="Probability",
values=c("white", "#FF000050", "#FF0000A0", "#FF0000A0", "#FF000050", "white"),
labels = c("(0, 0.05]", "(0.05, 0.25]",
"(0.25, 0.5]", "(0.5, 0.75]",
"(0.75, 0.95]", "(0.95, 1]")) +
theme_bw() +
theme(
legend.position = "none",
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.title.y = element_blank()) +
labs(x="Standardized PK Parameters", y="") +
scale_x_log10() +
coord_cartesian(expand=FALSE)
p3+p4
Ranges of BSV for each PK Parameter:
bsvranges <- refbsv[,list(
P05 = quantile(stdparamvalue, 0.05),
P25 = quantile(stdparamvalue, 0.25),
P50 = quantile(stdparamvalue, 0.5),
P75 = quantile(stdparamvalue, 0.75),
P95 = quantile(stdparamvalue, 0.95)), by = paramname]
bsvranges| paramname | P05 | P25 | P50 | P75 | P95 |
|---|---|---|---|---|---|
| Cmax | 0.7918743 | 0.9184431 | 1 | 1.077698 | 1.185392 |
| Clast | 0.4752010 | 0.7899361 | 1 | 1.286067 | 1.735681 |
| AUC | 0.7063883 | 0.8917846 | 1 | 1.123854 | 1.313933 |
