Bioequivalence and Bioavailability Forum

Dear all!

Having thought a long time about handling endogenous compounds I would like to share my approach how to model such data assuming a two-compartmental model following an IV bolus. I think the model below has an advantage compared to the frequently used turnover model based on differential equations (e.g. Gabrielsson and Weiner, 2000) as the presented model takes also the measured pre-dose concentration(s) straightforwardly into account.

I would be grateful if you could have a look at the code and share your opinion on this PK model.

best regards

martin

Gabrielsson J. and Weiner D. (2000). Pharmacokinetic and Pharmacodynamic Data Analysis: Concepts and Applications. 4th Edition. Swedish Pharmaceutical Press, Stockholm.

#### modified data from Gabrielsson and Weiner (2000, page 743)
#### endogenous concentration is assumed to be constant over time 
dose <- 36630 
time <- c(-1, 0.167E-01, 0.1167, 0.1670, 0.25, 0.583, 0.8330, 1.083, 1.583, 2.083, 4.083, 8.083, 12, 23.5, 24.25, 26.75, 32)
conc <- c(20.34, 3683, 884.7, 481.1, 215.6, 114, 95.8, 87.89, 60.19, 60.17, 34.89, 20.99, 20.54, 19.28, 18.18, 19.39, 22.72)
data <- data.frame(conc,time)

## specify indicator variable enabling inclusion of pre-dose concentration for fitting 
data$i1 <- ifelse(data$time <0, 1, 0)
data$i2 <- ifelse(data$time <0, 0, 1)

## assuming constant absolute error: ordinary least squares
mod.ols <- nls(conc ~ i1*base + i2*(base + a1*exp(-k1*time) + a2*exp(-k2*time)),
               start=c(base=20.34, a1=4500, k1=15, a2=125, k2=0.6), data=data, trace=TRUE)
 
## assuming constant relative error: weighted least squares 
mod.wls <- nls(conc ~ i1*base + i2*(base + a1*exp(-k1*time) + a2*exp(-k2*time)),
               start=c(base=20.34, a1=4500, k1=15, a2=125, k2=0.6), data=data,
               weight=1/predict(mod.ols)^2, trace=TRUE)

## assuming constant relative error: iteratively re-weighted least squares
mod.irwls <- mod.wls
for(i in 1:10){
   mod.irwls <- nls(conc ~ i1*base + i2*(base + a1*exp(-k1*time) + a2*exp(-k2*time)), 
                     start=c(base=20.34, a1=4500, k1=15, a2=125, k2=0.6),
                     data=data, weight=1/predict(mod.irwls)^2)
   print(as.vector(coef(mod.irwls)))
}
 
summary(mod.ols)
summary(mod.wls)
summary(mod.irwls)

newdata <- data.frame(time=seq(0,32,0.01))
newdata$i1 <- ifelse(newdata$time <0, 1, 0)
newdata$i2 <- ifelse(newdata$time <0, 0, 1)
plot(conc ~ time, data=data, ylim=c(10,1E4), log='y', yaxt='n',
     xlab='Time (hours)', ylab='Log of concentration (pmol/L)')
axis(side=2, at=c(10, 100, 1000, 10000), las=1)
axis(side=2, at=seq(1E1,1E2,1E1), tcl=-0.25, labels=FALSE)
axis(side=2, at=seq(1E2,1E3,1E2), tcl=-0.25, labels=FALSE)
axis(side=2, at=seq(1E3,1E4,1E3), tcl=-0.25, labels=FALSE)
points(x=newdata$time, y=predict(mod.irwls, newdata), type='l')