Bioequivalence and Bioavailability Forum

Dear all,

I tried to implement PK model I (one-compartment open, lag-time):^*

parameter                  distribution  mean    CV(%)  truncation
Volume of distribution     log-normal    1        10
elimination rate constant  log-normal    0.17375  20
absorption rate constant   log-normal    1.39     20
bioavailability            uniform       0.5      —     0.4–0.6
lag time                   normal        1        50    0–2

LLOQ 1% of theoretical C_max, analytical error: multiplicative component (log-normal, assay error 10%), additive component (log-normal, LLOQ × assay error).

My clumsy code (requires package truncnorm for the truncated normal distribution):
#set.seed(29081957)# uncomment this line to compare results only
require(truncnorm)
n      <- 16       # no of sims
t      <- c(0,0.25,0.5,0.75,1,1.25,1.5,2,2.5,3,3.5,4,5,6,9,12,16,24)
D      <- 500
# Index ".d" denotes default values, ".c" CV
F.d    <- 0.5      # fraction absorbed (BA)
F.c    <- 0.1      # ±0.1, uniform
V.d    <- 1        # volume of distribution
V.c    <- 0.1      # CV 10%, lognormal
k01.d  <- 1.39     # absorption rate constant (Wagner: ka)
k01.c  <- 0.2      # CV 20%, lognormal
k10.d  <- 0.17375  # elimination rate constant (Wagner: K)
k10.c  <- 0.2      # CV 20%, lognormal
tlag.d <- 1
tlag.c <- 0.5      # CV 50%, truncated normal
tlag.l <- 0        # lower truncation 0
tlag.u <- 2        # upper truncation 2
AErr   <- 0.1      # analytical error CV 10%, lognormal
curve(F.d*D/V.d*k01.d/(k01.d-k10.d)*(exp(-k10.d*(x-tlag.d))
  -exp(-k01.d*(x-tlag.d))), from=tlag.d, to=24, n=(24-tlag.d)*5+1,
  type="l", lwd=2, col="red", xlim=c(0,24), ylim=c(0,350),
  xlab="time", ylab="concentration", add=FALSE)
tmax.th<- log(k01.d/k10.d)/(k01.d-k10.d)+tlag.d
Cmax.th<- F.d*D*k01.d/(V.d*(k01.d-k10.d))*(exp(-k10.d*(tmax.th-tlag.d))
  -exp(-k01.d*(tmax.th-tlag.d)))
LLOQ   <- 0.01*Cmax.th # 1% of theoretical Cmax
abline(v=0)
abline(h=0)
abline(h=LLOQ, lty=3)
for (i in 1:n){
  F      <- runif(n=1, min=F.d-F.c, max=F.d+F.c)
  V      <- rlnorm(n=1, meanlog=log(V.d)-0.5*log(V.c^2+1),
              sdlog=sqrt(log(V.c^2+1)))
  k01    <- rlnorm(n=1, meanlog=log(k01.d)-0.5*log(k01.c^2+1),
              sdlog=sqrt(log(k01.c^2+1)))
  k10    <- rlnorm(n=1, meanlog=log(k10.d)-0.5*log(k10.c^2+1),
              sdlog=sqrt(log(k10.c^2+1)))
  tlag   <- rtruncnorm(n=1, a=tlag.l, b=tlag.u, mean=tlag.d, sd=tlag.c)
  # C    <- F*D/V*k01/(k01-k10)*(exp(-k10*(t-tlag))-exp(-k01*(t-tlag)))
  # explicit (compute macro-constants, Wagner’s notation)
  # John G. Wagner, Pharmacokinetics for the Pharmaceutical Scientist, pp 4–5, 1993)
  ka     <- k01
  K      <- k10
  C0     <- F*D/V  # fictive iv concentration at t=0
  B1     <- C0*ka/(ka-K)*exp(K*tlag)
  B2     <- C0*ka/(ka-K)*exp(ka*tlag)
  C      <- B1*exp(-K*t)-B2*exp(-ka*t)
  for (j in 1:length(t)){ # analytical error (multiplicative 10% C)
    if (j==1)
      AErr1 <- rlnorm(n=1, meanlog=0, sdlog=sqrt(log(C[j]*AErr^2+1)))
    else
      AErr1 <- c(AErr1, rlnorm(n=1, meanlog=0, sdlog=sqrt(log(C[j]*AErr^2+1))))
  }
  AErr2  <- rep(AErr*LLOQ, length(t)) # analytical error (constant 10% LLOQ)
  C      <- C+AErr1+AErr2 # add analytical error
  C[C<LLOQ] <- NA  # assign NAs to Cs below LLOQ
  if (i==1) P <- log(C) else P <- P+log(C) # Prod. for geom. means
  if (i==1) S <- C else S <- S+C           # Sum  for arithm. means
  lines(t, C, type="l") # spaghetti plots
  }
GM       <- exp(P/n)  # geometric means
AM       <- S/n       # arithmetic means
points(t, GM, pch=16, col="blue")
lines(t, GM, type="l", col="blue", lwd=2)
points(t, AM, pch=16, col="#12DD00")
lines(t, AM, type="l", col="#12DD00", lwd=2)
legend("topright",
  c("theoretical", "geometric means", "arithmetic means"),
  pch=c(3,16,16),
  col=c("red","blue","#12DD00"))

A run of the set random seed below:



Is this correct? I think that I screwed up the analytical error. Original text:

Analytical assay errors were generated from log-normal distributions with no bias, a CV of 10%, plus a constant term equal to the product of the assay CV and the limit of quantification, LQ.

Shouldn’t I rather use a normal distribution instead (AErr1 <- rnorm(n=1, mean=0, sd=abs(C[j]*AErr)))? …which would give:

---

• Csizmadia F, Endrényi L. Model-Independent Estimation of Lag Times with First-Order Absorption and Disposition. J Pharm Sci. 1998;87(5):608–12. doi:10.1021/js9703333