Specifying a PK Model using mrgsolve
Here we illustrate the approach with a two-compartment PK model defined with an ODE and covariates on the PK parameters Clearance and Volume.
codepkmodelcov <- '
$PARAM
CL = 4, V=10 , KA=0.5, Vp =50, Qp= 10,
CLALB = -0.8, CLSEX = 0.2, CLWT = 1,
VSEX = 0.07, VWT = 1,
WT=85, SEX=0, AGE=45, ALB =45
$CMT GUT CENT PER
$MAIN
double CLi = CL *
pow((ALB/45.0), CLALB)*
(SEX == 1.0 ? (1.0+CLSEX) : 1.0)*
pow((WT/85.0), CLWT)*exp(ETA(1));
double Vi = V *
(SEX == 1.0 ? (1.0+VSEX) : 1.0)*
pow((WT/85.0), VWT)*exp(ETA(2));
double KAi = KA;
double Vpi = Vp *pow((WT/70.0), 1);
double Qpi = Qp *pow((WT/70.0), 0.75);
double Keli = CLi/Vi ;
double Kpti = Qpi/Vi ;
double Ktpi = Qpi/Vpi ;
$OMEGA
0.3
0.01 0.3
$ODE
dxdt_GUT = -KAi*GUT;
dxdt_CENT = KAi*GUT-Keli*CENT-Kpti*CENT+Ktpi*PER;
dxdt_PER = Kpti*CENT-Ktpi*PER;
$TABLE
double CP = CENT/ Vi;
double CPER = PER/Vpi;
$CAPTURE CP CPER KAi CLi Vi Vpi Qpi
'
modcovsim <- mcode("codepkmodelcov", codepkmodelcov)
#> Building codepkmodelcov ... done.Simulate Reference Subjects with BSV
We simulate with between subject variability (BSV) and provide a plot of concentrations on linear and log linear scales.
idata <- tibble(
ID = 1:1000,
WT = 85,
ALB = 45,
SEX = 0
)
ev1 <- ev(time = 0, amt = 100, cmt = 1)
data.dose <- ev(ev1)
data.dose <- as.data.frame(data.dose)
data.all <- merge(idata, data.dose)
outputsim <- modcovsim %>%
data_set(data.all) %>%
carry.out(WT, AGE, SEX, ALB, CLi) %>%
mrgsim(end = 24, delta = 0.25)
outputsim <- as.data.frame(outputsim)
outputsim <- outputsim %>%
arrange(ID, time)
outputsim$SEX <- as.factor(outputsim$SEX)
set.seed(678549)
p1 <- ggplot(data = outputsim[outputsim$ID %in% sample(unique(outputsim$ID), 500), ],
aes(time, CP, group = ID)) +
geom_line(alpha = 0.2, size = 0.1) +
facet_grid(ALB ~ WT+SEX,
labeller = label_both) +
labs(y = "Plasma Concentrations", color = "Sex", x = "Time (h)")
set.seed(678549)
p2 <- ggplot(data = outputsim[outputsim$ID %in% sample(unique(outputsim$ID), 500), ],
aes(time, CP, group = ID)) +
geom_line(alpha = 0.2, size = 0.1) +
facet_grid(ALB ~ WT+SEX,
labeller = label_both) +
scale_y_log10() +
labs(y = expression(Log[10]~~Plasma~~Concentrations), color = "Sex", x = "Time (h)")
p1+p2
#> Warning: Transformation introduced infinite values in continuous y-axis
Compute NCA PK Parameters, Plot and Summarize BSV
We then compute simple NCA parameters, 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 NCA parameter.
NCATYPICAL <- outputsim %>%
group_by(ID, ALB, WT, SEX) %>%
summarise (
Cmax = max(CP, na.rm = TRUE),
Clast = CP[n()],
AUC = sum(diff(time) * na.omit(lead(CP) + CP)) / 2
) %>%
gather(paramname, paramvalue, Cmax, Clast, AUC)
p3 <- ggplot(NCATYPICAL,
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("#FF0000A0",
"#0000FFA0", "white","white",
"#0000FFA0","#FF0000A0"),
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()+
labs(x="NCA parameters",y="")+
scale_x_log10()+
coord_cartesian(expand = FALSE)
NCATYPICALREF <- NCATYPICAL%>%
group_by (paramname) %>%
mutate(medparam = median(paramvalue),
paramvalue = paramvalue / medparam)
p4 <- ggplot(NCATYPICALREF,
aes(x=paramvalue,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("#FF0000A0","#0000FFA0", "white","white", "#0000FFA0","#FF0000A0"),
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")+
labs(x="Standardized NCA parameters",y="")+
scale_x_log10()+
coord_cartesian(expand = FALSE)
p3+p4
BSVRANGES<- NCATYPICALREF %>%
summarize(
P05 = quantile(paramvalue, 0.05),
P25 = quantile(paramvalue, 0.25),
P50 = quantile(paramvalue, 0.5),
P75 = quantile(paramvalue, 0.75),
P95 = quantile(paramvalue, 0.95)
)
BSVRANGES| paramname | P05 | P25 | P50 | P75 | P95 |
|---|---|---|---|---|---|
| AUC | 0.5248993 | 0.7945283 | 1 | 1.215652 | 1.486990 |
| Clast | 0.2432740 | 0.6147814 | 1 | 1.477597 | 2.190552 |
| Cmax | 0.7851306 | 0.9229614 | 1 | 1.068844 | 1.150875 |
