Building on the first simple example, we can be more ambitious, and try a full PBPK model. This one was published for mavoglurant (Wendling et al. 2016).

# nlmixr model

library(nlmixr)

library(xpose)
library(xpose.nlmixr)
library(ggplot2)

pbpk <- function(){
ini({
##theta=exp(c(1.1, .3, 2, 7.6, .003, .3))
lKbBR = 1.1
lKbMU = 0.3
lCLint = 7.6
lKbBO = 0.03
lKbRB = 0.3
eta.LClint ~ 4
prop.err <- 10
})
model({
KbBR = exp(lKbBR)
KbMU = exp(lKbMU)
CLint= exp(lCLint + eta.LClint)
KbBO = exp(lKbBO)
KbRB = exp(lKbRB)

## Regional blood flows
CO  = (187.00*WT^0.81)*60/1000;         # Cardiac output (L/h) from White et al (1968)
QHT = 4.0 *CO/100;
QBR = 12.0*CO/100;
QMU = 17.0*CO/100;
QSK = 5.0 *CO/100;
QSP = 3.0 *CO/100;
QPA = 1.0 *CO/100;
QLI = 25.5*CO/100;
QST = 1.0 *CO/100;
QGU = 14.0*CO/100;
QHA = QLI - (QSP + QPA + QST + QGU); # Hepatic artery blood flow
QBO = 5.0 *CO/100;
QKI = 19.0*CO/100;
QRB = CO - (QHT + QBR + QMU + QAD + QSK + QLI + QBO + QKI);
QLU = QHT + QBR + QMU + QAD + QSK + QLI + QBO + QKI + QRB;

## Organs' volumes = organs' weights / organs' density
VLU = (0.76 *WT/100)/1.051;
VHT = (0.47 *WT/100)/1.030;
VBR = (2.00 *WT/100)/1.036;
VMU = (40.00*WT/100)/1.041;
VSK = (3.71 *WT/100)/1.116;
VSP = (0.26 *WT/100)/1.054;
VPA = (0.14 *WT/100)/1.045;
VLI = (2.57 *WT/100)/1.040;
VST = (0.21 *WT/100)/1.050;
VGU = (1.44 *WT/100)/1.043;
VBO = (14.29*WT/100)/1.990;
VKI = (0.44 *WT/100)/1.050;
VAB = (2.81 *WT/100)/1.040;
VVB = (5.62 *WT/100)/1.040;
VRB = (3.86 *WT/100)/1.040;

## Fixed parameters
BP = 0.61;      # Blood:plasma partition coefficient
fup = 0.028;    # Fraction unbound in plasma
fub = fup/BP;   # Fraction unbound in blood

KbLU = exp(0.8334);
KbHT = exp(1.1205);
KbSK = exp(-.5238);
KbSP = exp(0.3224);
KbPA = exp(0.3224);
KbLI = exp(1.7604);
KbST = exp(0.3224);
KbGU = exp(1.2026);
KbKI = exp(1.3171);

##-----------------------------------------
S15 = VVB*BP/1000;
C15 = Venous_Blood/S15

##-----------------------------------------
d/dt(Lungs) = QLU*(Venous_Blood/VVB - Lungs/KbLU/VLU);
d/dt(Heart) = QHT*(Arterial_Blood/VAB - Heart/KbHT/VHT);
d/dt(Brain) = QBR*(Arterial_Blood/VAB - Brain/KbBR/VBR);
d/dt(Muscles) = QMU*(Arterial_Blood/VAB - Muscles/KbMU/VMU);
d/dt(Skin) = QSK*(Arterial_Blood/VAB - Skin/KbSK/VSK);
d/dt(Spleen) = QSP*(Arterial_Blood/VAB - Spleen/KbSP/VSP);
d/dt(Pancreas) = QPA*(Arterial_Blood/VAB - Pancreas/KbPA/VPA);
d/dt(Liver) = QHA*Arterial_Blood/VAB + QSP*Spleen/KbSP/VSP + QPA*Pancreas/KbPA/VPA + QST*Stomach/KbST/VST + QGU*Gut/KbGU/VGU - CLint*fub*Liver/KbLI/VLI - QLI*Liver/KbLI/VLI;
d/dt(Stomach) = QST*(Arterial_Blood/VAB - Stomach/KbST/VST);
d/dt(Gut) = QGU*(Arterial_Blood/VAB - Gut/KbGU/VGU);
d/dt(Bones) = QBO*(Arterial_Blood/VAB - Bones/KbBO/VBO);
d/dt(Kidneys) = QKI*(Arterial_Blood/VAB - Kidneys/KbKI/VKI);
d/dt(Arterial_Blood) = QLU*(Lungs/KbLU/VLU - Arterial_Blood/VAB);
d/dt(Rest_of_Body) = QRB*(Arterial_Blood/VAB - Rest_of_Body/KbRB/VRB);

})
}

dat$occ = unlist(with(dat, tapply(EVID, ID, function(x) cumsum(x>0)))) dat = subset(dat, occ==1) dat = subset(dat, ID<812) ## First 20 dat = subset(dat, EVID>0 | DV>0) dat$CMT[dat$CMT == 0] <- 1; dat$CMT[dat$EVID == 1] <- "Venous_Blood" ## Compartment dosed to is Venous Blood dat$CMT[dat$EVID != 1] <- "C15" ## Observing C15 gofs <- function(fit){ ################################################################################ ## Standard plots ################################################################################ plot(fit); xpdb <- xpose_data_nlmixr(fit) ## Convert to nlmixr object print(dv_vs_pred(xpdb) + ylab("Observed Mavoglurant Concentrations (ng/mL)") + xlab("Population Predicted Mavoglurant Concentrations (ng/mL)")); print(dv_vs_ipred(xpdb) + ylab("Observed Mavoglurant Concentrations (ng/mL)") + xlab("Individual Predicted Mavoglurant Concentrations (ng/mL)")); print(res_vs_pred(xpdb) + ylab("Conditional Weighted Residuals") + xlab("Population Predicted Mavoglurant Concentrations (ng/mL)")); print(res_vs_idv(xpdb) + ylab("Conditional Weighted Residuals") + xlab("Time (h)")); if (!is.null(fit$saem)){
print(prm_vs_iteration(xpdb));
}

print(absval_res_vs_idv(xpdb, res = 'IWRES') +
ylab("Individual Weighted Residuals") +
xlab("Time (h)"))

print(absval_res_vs_pred(xpdb, res = 'IWRES') +
ylab("Individual Weighted Residuals") +
xlab("Population Predicted Mavoglurant Concentrations (ng/mL)"))

print(ind_plots(xpdb, nrow=3, ncol=4) +
ylab("Predicted and Observed Mavoglurant Concentrations (ng/mL)") +
xlab("Time (h)"))

print(res_distrib(xpdb) +
ylab("Density") +
xlab("Conditional Weighted Residuals"));
##Visual Predictive Checks
f1 <- vpc.ui(fit,n=500,stratify=c("dose"), show=list(obs_dv=T), log_y=TRUE,
bins = c(0, 2, 4, 6, 8, 10, 20, 30, 40, 50),
ylab = "Mavoglurant Concentrations (ng/mL)",
xlab = "Time (hours)")
f2 <- vpc.ui(fit,n=500, show=list(obs_dv=T), bins = c(0, 2, 4, 6, 8, 10, 20, 30, 40, 50), log_y=TRUE,
ylab = "Mavoglurant Concentrations (ng/mL)", xlab = "Time (hours)")
plot(f1)
plot(f2)
}

# SAEM

fit.addProp.S <- nlmixr(pbpk, dat, est="saem", control=list(print=0),
table=list(cwres=TRUE, npde=TRUE))
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#>
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#>
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#>
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#>
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#>
#> [1] "WT"

gofs(fit.addProp.S);

## Change error to lognormal

fit.lnorm.S <- pbpk %>%
model({C15 ~ lnorm(lnorm.err)}) %>% # Change C15 to be log-normally distributed
## Requires data since piping from pbpk model
nlmixr(dat,est="saem", control=list(print=0),
table=list(cwres=TRUE, npde=TRUE))
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#>
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#>
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#>
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#>
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#>
#> [1] "WT"

NOTE: lognormal distribution AIC/loglik/etc is on normal scale. Therefore, you can compare the AICs between fit.lnorm and fit.addProp since they are calculated on the same scale.

In this case you can see that the AIC for the log-normal model is better than the AIC for the addProp model.

gofs(fit.lnorm.S);

# Piping to FOCEi

You can pipe models from different estimation methods to new estimation methods.

fit.addProp.F <- fit.addProp.S %>%
nlmixr(est="focei",
control=list(print=0),
table=list(cwres=TRUE, npde=TRUE));
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#>
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#>
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#>
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#>
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#>
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#>
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#>
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#>
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#>
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#>
#> [1] "WT"
#> calculating covariance matrix
#> [====|====|====|====|====|====|====|====|====|====] 0:01:57
#> done
## Since this was a model pipline, the data
## remains the same as the last fit.

gofs(fit.addProp.F);

## Lognoral

fit.lnorm.F <- fit.addProp.F %>%
model({C15 ~ lnorm(lnorm.err)}) %>%
nlmixr(est="focei",
control=list(print=0),
table=list(cwres=TRUE, npde=TRUE));
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#>
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#>
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#>
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#>
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#>
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#>
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#>
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#>
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#>
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#>
#> [1] "WT"
#> calculating covariance matrix
#> [====|====|====|====|====|====|====|====|====|====] 0:00:30
#> done
## In this model pipline we are changing the fit method to focei.

gofs(fit.lnorm.F);

# Traditional lognormal estimates are identical

datL  <- dat
datL$DV <- log(datL$DV);

pbpkL <- function() {
ini({
##theta=exp(c(1.1, .3, 2, 7.6, .003, .3))
lKbBR = 1.1
lKbMU = 0.3
lCLint = 7.6
lKbBO = 0.03
lKbRB = 0.3
eta.LClint ~ 4
})
model({
KbBR = exp(lKbBR)
KbMU = exp(lKbMU)
CLint= exp(lCLint + eta.LClint)
KbBO = exp(lKbBO)
KbRB = exp(lKbRB)

## Regional blood flows
CO  = (187.00*WT^0.81)*60/1000;         # Cardiac output (L/h) from White et al (1968)
QHT = 4.0 *CO/100;
QBR = 12.0*CO/100;
QMU = 17.0*CO/100;
QSK = 5.0 *CO/100;
QSP = 3.0 *CO/100;
QPA = 1.0 *CO/100;
QLI = 25.5*CO/100;
QST = 1.0 *CO/100;
QGU = 14.0*CO/100;
QHA = QLI - (QSP + QPA + QST + QGU); # Hepatic artery blood flow
QBO = 5.0 *CO/100;
QKI = 19.0*CO/100;
QRB = CO - (QHT + QBR + QMU + QAD + QSK + QLI + QBO + QKI);
QLU = QHT + QBR + QMU + QAD + QSK + QLI + QBO + QKI + QRB;

## Organs' volumes = organs' weights / organs' density
VLU = (0.76 *WT/100)/1.051;
VHT = (0.47 *WT/100)/1.030;
VBR = (2.00 *WT/100)/1.036;
VMU = (40.00*WT/100)/1.041;
VSK = (3.71 *WT/100)/1.116;
VSP = (0.26 *WT/100)/1.054;
VPA = (0.14 *WT/100)/1.045;
VLI = (2.57 *WT/100)/1.040;
VST = (0.21 *WT/100)/1.050;
VGU = (1.44 *WT/100)/1.043;
VBO = (14.29*WT/100)/1.990;
VKI = (0.44 *WT/100)/1.050;
VAB = (2.81 *WT/100)/1.040;
VVB = (5.62 *WT/100)/1.040;
VRB = (3.86 *WT/100)/1.040;

## Fixed parameters
BP = 0.61;      # Blood:plasma partition coefficient
fup = 0.028;    # Fraction unbound in plasma
fub = fup/BP;   # Fraction unbound in blood

KbLU = exp(0.8334);
KbHT = exp(1.1205);
KbSK = exp(-.5238);
KbSP = exp(0.3224);
KbPA = exp(0.3224);
KbLI = exp(1.7604);
KbST = exp(0.3224);
KbGU = exp(1.2026);
KbKI = exp(1.3171);

##-----------------------------------------
S15 = VVB*BP/1000;
C15 = Venous_Blood/S15
lnC15 = log(C15);

##-----------------------------------------
d/dt(Lungs) = QLU*(Venous_Blood/VVB - Lungs/KbLU/VLU);
d/dt(Heart) = QHT*(Arterial_Blood/VAB - Heart/KbHT/VHT);
d/dt(Brain) = QBR*(Arterial_Blood/VAB - Brain/KbBR/VBR);
d/dt(Muscles) = QMU*(Arterial_Blood/VAB - Muscles/KbMU/VMU);
d/dt(Skin) = QSK*(Arterial_Blood/VAB - Skin/KbSK/VSK);
d/dt(Spleen) = QSP*(Arterial_Blood/VAB - Spleen/KbSP/VSP);
d/dt(Pancreas) = QPA*(Arterial_Blood/VAB - Pancreas/KbPA/VPA);
d/dt(Liver) = QHA*Arterial_Blood/VAB + QSP*Spleen/KbSP/VSP + QPA*Pancreas/KbPA/VPA + QST*Stomach/KbST/VST + QGU*Gut/KbGU/VGU - CLint*fub*Liver/KbLI/VLI - QLI*Liver/KbLI/VLI;
d/dt(Stomach) = QST*(Arterial_Blood/VAB - Stomach/KbST/VST);
d/dt(Gut) = QGU*(Arterial_Blood/VAB - Gut/KbGU/VGU);
d/dt(Bones) = QBO*(Arterial_Blood/VAB - Bones/KbBO/VBO);
d/dt(Kidneys) = QKI*(Arterial_Blood/VAB - Kidneys/KbKI/VKI);
d/dt(Arterial_Blood) = QLU*(Lungs/KbLU/VLU - Arterial_Blood/VAB);
d/dt(Venous_Blood) = QHT*Heart/KbHT/VHT + QBR*Brain/KbBR/VBR +
QSK*Skin/KbSK/VSK + QLI*Liver/KbLI/VLI + QBO*Bones/KbBO/VBO +
QKI*Kidneys/KbKI/VKI + QRB*Rest_of_Body/KbRB/VRB - QLU*Venous_Blood/VVB;
d/dt(Rest_of_Body) = QRB*(Arterial_Blood/VAB - Rest_of_Body/KbRB/VRB);

})
}

fit.lnorm.trans <- pbpkL %>%
nlmixr(datL,est="saem",
control=list(print=0),
table=list(npde=TRUE, cwres=TRUE))
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#>
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#>
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#>
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#>
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#>
#> [1] "WT"

NOTE: the estimates are the same but the AIC is different since it is calculated on the log scale.