Past

“In steady-state studies, wash-out of the last dose of the previous treatment can overlap with the build-up of the second treatment, provided the build-up period is sufficiently long (at least three times the dominating half-life).
“Whenever multiple dose studies are performed it should be demonstrated that steady state has been reached.
“In steady-state studies wash out of the previous treatment last dose can overlap with the build-up of the second treatment, provided the build-up period is sufficiently long (at least three times the terminal half-life).

A similar wording is used by the WHO,⁸ only that »at least five times the terminal half-life« are recommended.

In Europe regularly a linear regression of the last three pre-dose concentrations was performed, allowing to exclude subjects who are not in steady state. The idea behind is based on elementary pharmacokinetics. By definition in steady state the rate of inflow equals the rate of outflow. Hence, in all dosing intervals the pre-dose concentrations (\(\small{C_{\textrm{pd}}}\)) are expected to be similar. Therefore, a linear regression should give a slope close to zero. However, this approach was based on custom and not backed by a guideline.

Recall the accumulation index \(\small{R=1/(1-2^{-\tau\,/\,t_{1/2}}\,)}\). \(\small{R}\) is strictly valid for intravenous doses and approximate for extravascular administrations in a one-compartment model if \(\small{k_{01}\gg k_{10}}\).
An example with \(\small{\tau=t_{1/2}\;(R=2)}\), where \(\small{C_\textrm{pd}}\) is the pre-dose concentration and \(\small{C_\tau}\) is the concentration at the end of the dosing interval. Values ‘by definition’ are denoted by \(\small{\ast}\). \[\small{ \begin{array}{rrrrrr} \text{dose} & \text{time} & C_\textrm{pd}\phantom{00000} & \text{% change} & \text{moving slope} & C_\tau\phantom{000000}\\ \hline 1 & 0 & 0\phantom{.000000}\phantom{0} & - & - & 50\phantom{.00000000}\\ 2 & 24 & 50\phantom{.0000000} & 50\phantom{.0000000} & - & 75\phantom{.00000000}\\ 3 & 48 & 75\phantom{.0000000} & 16.6666667 & 0.78125 & 87.5\phantom{0000000}\\ 4 & 72 & 87.5\phantom{000000} & 7.1428571 & 0.39063 & 93.75\phantom{000000}\\ 5 & 96 & 93.75\phantom{00000} & 3.3333333 & 0.19531 & 96.875\phantom{00000}\\ 6 & 120 & 96.875\phantom{0000} & 1.6129032 & 0.09766 & 98.4375\phantom{0000}\\ 7 & 144 & 98.4375\phantom{000} & 0.7936508 & 0.04883 & 99.21875\phantom{000}\\ 8 & 168 & 99.21875\phantom{00} & 0.3937008 & 0.02441 & 99.609375\phantom{00}\\ 9 & 192 & 99.609375\phantom{0} & 0.1960784 & 0.01221 & 99.8046875\phantom{0}\\ 10 & 216 & 99.8046875 & 0.0978474 & 0.00610 & 99.90234375\\ \infty & \infty & 100\phantom{.0000000} & 0\ast\phantom{.00000} & 0\ast\phantom{.000} & 100\;(=C_\textrm{pd})\phantom{0}\\ \hline \end{array}}\]

Health Canada required 1992 – 2011 multiple dose studies for drugs which accumulated, i.e., if the extrapolated \(\small{AUC}\) after a single dose was > 20% of \(\small{AUC_{0-\infty}}\).

“7.6.2 Evidence of Steady State
For steady-state studies of drugs with uncomplicated characteristics, at least three consecutive pre-dose concentration levels (C_pd) are required to provide evidence of steady state. Generally, observations of C_pd for the test and reference products should be recorded at the same time of the day. One of these measurements could be taken based on the first sample of the study day in which a profile over the dosing interval is being established. Steady state is usually achieved when repeated doses of a formulation are administered over a period that exceeds five disposition half-lives of the modified-release form.

‘Evidence of steady state’ had to be assessed by a repeated measures analysis with an amazing number of 11 (eleven‼) factors. If you want to give it a try, download the data compiled from the guidance’s Table 12-A (Ran­do­mization), Table 12-G (Test), and Table 12-H (Reference).

data     <- read.csv("https://bebac.at/downloads/HC1996-TABLE-12.csv",
                     colClasses = c(rep("factor", 5), "integer"))
geo.mean <- function(x) {
  t  <- as.numeric(unlist(unique(x[1])))
  n  <- as.integer(unlist(unique(x[2])))
  gm <- NA_real_
  for (i in seq_along(t)) {
    y <- x[x[1] == t[i], ]
    if (sum(y$C <= 0) > 0) {
      y <- y[-which(y$C <= 0), ]
    }
    if (nrow(y) > 0) {
      gm[i] <- exp(mean(log(y$C), na.rm = TRUE))
    }
  }
  return(gm)
}
Subject <- unique(as.character(data$Subject))
Form    <- unique(as.character(data$Formulation))
slopes  <- data.frame(Formulation = rep(Form, each = length(Subject)),
                      Subject = Subject, slope = NA_real_,
                      lower.CL = NA_real_, upper.CL = NA_real_,
                      sig = FALSE)
t       <- data$Time
t       <- as.numeric(levels(t))[t]
ylim    <- range(data$C.pd)
xlim    <- range(t)
cols    <- c("blue", "red")
col     <- colorRampPalette(cols)(length(Subject))
col     <- paste0(col, "80")
op      <- par(no.readonly=TRUE)
par(mar = c(4.1, 4, 0.2, 0.2), cex.axis = 0.9)
split.screen(c(2, 1))
for (i in seq_along(Form)) {
  screen(i)
  ifelse (i == 1, xlab  <- "", xlab <- "Time (day)")
  plot(xlim, ylim, type = "n", las = 1, xlab = xlab,
       ylab = paste0("C.pd  (", Form[i], ")"))
  grid()
  x <- data.frame(t = t,
                  C = data[data$Formulation == Form[i], 6])
  lines(unique(t), geo.mean(x), lwd = 2)
  for (j in seq_along(Subject)) {
    t  <- data$Time[data$Formulation == Form[i] &
                    data$Subject == Subject[j]]
    t  <- as.numeric(levels(t))[t]
    C  <- data$C.pd[data$Formulation == Form[i] &
                    data$Subject == Subject[j]]
    points(t, C, pch = 16, cex = 1.1, col = col[j])
    x  <- data.frame(t = t, C = C)
    x  <- tail(x, 3)
    m  <- lm(C ~ t, data = x)
    lines(x$t, predict(m), col = substr(col[j], 1, 7))
    slopes$slope[slopes$Formulation == Form[i] &
                 slopes$Subject == Subject[j]] <- coef(m)[2]
    ci <- suppressWarnings(confint(m, "t"))
    slopes[slopes$Formulation == Form[i] &
           slopes$Subject == Subject[j], 4:5] <- sprintf("%+.4f", ci)
    if (sign(ci[1]) == sign(ci[2]))
      slopes$sig[slopes$Formulation == Form[i] &
                 slopes$Subject == Subject[j]] <- TRUE
  }
}
close.screen(all = TRUE)
par(op)
slopes$slope     <- sprintf("%+5.1f", slopes$slope)
names(slopes)[6] <- "sign # 0"
last.3 <- data[as.numeric(as.character(data$Time)) >= 6, ]
muddle <- lm(log(C.pd) ~ Sequence + Subject * Sequence + Period +
                         Formulation + Subject * Formulation * Sequence + 
                         Time + Time * Sequence + Time * Subject * Sequence +
                         Time * Period + Time * Period +
                         Time * Formulation * Subject * Sequence,
                         data = last.3)
print(slopes, row.names = FALSE); anova(muddle) # Oops!

“Of main importance in this analysis is the time and time*form interaction. The error term to test time is the time*ID(Seq). The error term to test time*form is time*form*ID(Seq). Should either the time or time*form effects be significant, the study may not have been at steady state for one or both formulations.

Back in the day I failed to reproduce the results given in Table 12-I.
In R I got …

Analysis of Variance Table
Response: log(C.pd)
                         Df  Sum Sq  Mean Sq F value Pr(>F)
Sequence                  1 0.01326 0.013260     NaN    NaN
Subject                  10 1.22059 0.122059     NaN    NaN
Period                    1 0.01825 0.018249     NaN    NaN
Formulation               1 0.22069 0.220686     NaN    NaN
Time                      2 0.11935 0.059675     NaN    NaN
Subject:Formulation      10 0.56173 0.056173     NaN    NaN
Sequence:Time             2 0.13041 0.065204     NaN    NaN
Subject:Time             20 0.71386 0.035693     NaN    NaN
Period:Time               2 0.06441 0.032204     NaN    NaN
Formulation:Time          2 0.05961 0.029804     NaN    NaN
Subject:Formulation:Time 20 0.60824 0.030412     NaN    NaN
Residuals                 0 0.00000      NaN               
Warning message:
In anova.lm(m) : ANOVA F-tests on an essentially perfect fit are unreliable

… and in Phoenix WinNonlin (all effects fixed):

ERROR 11068: There are no degrees of freedom for residual. No tests are possible.

When I specified all subject-related effects as random:

      Hypothesis    Numer_DF    Denom_DF      F_stat      P_value
-----------------------------------------------------------------
             int           1          60     4446.56       0.0000
        Sequence           1          60        0.0609162  0.8059
     Formulation           1          60        3.04158    0.0863
            Time           2          60        1.00778    0.3711
   Time*Sequence           2          60        1.19822    0.3088
     Time*Period           3          60        1.59434    0.2002
Time*Formulation           2          60        1.64308    0.2020

\[\small{\text{TABLE 12-I }\textsf{(in all of its doubtful beauty)}}\] \[\small{\begin{array}{lrcccc} \hline\text{Source} & \text{df} & \text{SS} & \text{MS} & \text{F} & \text{PR > F}\\\hline \text{Seq} & 1 & 0.01500416 & 0.01500416 & 0.12 & 0.7313\\ \text{ID(Seq)} & 10 & 1.20256303 & 1.20256303 & 2.11 & 0.1278\\ \text{Period} & 1 & 0.04556284 & 0.04556284 & 0.80 & 0.3926\\ \text{Form} & 1 & 0.22620312 & 0.22620312 & 3.96 & 0.0745\\ \text{ID*Form(Seq)} & 10 & 0.57066567 & 0.0570666\phantom{0} & - & -\\ \text{Time} & 2 & 0.11182757 & 0.05591378 & 1.60 & 0.2274\\ \text{Time*Seq} & 10 & 0.13927056 & 0.06963528 & 1.99 & 0.1631\\ \text{Time*ID(Seq)} & 20 & 0.70058566 & 0.03502928 & - & -\\ \text{Time*Period} & 2 & 0.07672565 & 0.03836282 & 1.24 & 0.3110\\ \text{Time*Form} & 2 & 0.06390201 & 0.03195101 & 1.03 & 0.3746\\ \text{Time*Form*ID(Seq)} & 20 & 0.61937447 & 0.03096872 & - & -\\ \text{Corrected Total} & 71 & 3.74093818 & - & - & -\\\hline \end{array}}\] Maybe you are equipped with Harry Potter’s magic wand. If you succeed in reproducing this table, drop me a note how you did it. I’m always eager to learn.

However, the outcome of such an analysis is dichotomous (pass | fail), risking to discard the entire [sic] study. The family-wise error rate of two simultaneous tests (here \(\small{\text{Time}}\) and \(\small{\text{Time*Form}}\)) at level \(\small{\alpha}\) is given by \[FWER=1-(1-\alpha)^2 \tag{1}\] or 9.75% for \(\small{\alpha=0.05}\). That means, following this method you may have to trash almost one of ten studies although steady state was achieved.
Given, that is only of historical interest.

Current

“Achievement of steady state can be evaluated by collecting pre-dose samples on the day before the PK assessment day and on the PK assessment day. A specific statistical method to assure that steady state has been reached is not considered necessary in bioequivalence studies. Descriptive data is sufficient.

The requirement of performing multiple dose studies was lifted in 2012 by Health Canada. If such studies are performed, the pre-dose concentrations have to be reported.¹¹

“Because single-dose studies are considered more sensitive in addressing the primary question of BE (e.g., release of the drug substance from the drug product into the systemic circulation), mul­tiple-dose studies are generally not recommended.

A statistical method is not given in any guidance of the FDA. A statement about the current practice:¹³

“Current FDA practice for general procedure for steady-state assessment is by comparing at least three pre-dose concentrations for each formulation using linear regression techniques to evaluate achievement of steady state and by comparing the BE outcome with or without the subjects whose terminal [sic] slope is significantly different from 0 (i.e., they did not attain steady state).

“Whenever multiple dose studies are performed it should be demonstrated that steady state has been reached. Achievement of steady-state is assessed by comparing at least three pre-dose concentrations for each formulation, unless otherwise justified.

Since the EMA’s guideline is adopted in Australia, the same story Down Under.¹⁵ Even less details are given by the ICH.¹⁶

“[…] applicants should document appropriate dosage administration and sampling to demonstrate the attainment of steady-state. For steady-state studies, the following PK parameters should be tabulated: 1) primary parameters for analysis: C_maxSS and AUC_(0-tauSS), and 2) additional parameters for analysis: C_tauSS, C_minSS, C_avSS, degree of fluctuation, swing, and T_max.

No details are given by the ANVISA.¹⁷

“Devendo ser assegurada a obtenção do estado de equilíbrio antes das coletas das amostras.
Em estudos com doses múltiplas, a amostra pré-dose deve ser coletada imediatamente antes da última administração da medicação.
[It must be ensured that equilibrium state is reached before the samples are taken.
In multiple dose studies, the pre-dose sample should be collected immediately prior to the last drug administration.]

Linear model

For every subject a linear model is set up \[Y=a+bX+\epsilon\small{\textsf{,}}\tag{2}\] where \(\small{Y}\) is the vector of responses (here at least three pre-dose concentrations), \(\small{X}\) the vector of predictors (the corresponding sampling time points), and \(\epsilon\) the unexplained error with \(\small{\mathcal{N}(0;\sigma^2)}\). Fitting the data gives estimates of the intercept \(\small{\hat{a}}\) and slope \(\small{\hat{b}}\). One can either perform a one-sample t-test of \(\small{\hat{b}}\) versus zero with \(\small{\alpha=0.05}\) (significantly different if \(\small{p(t)<0.05}\)) or assess whether zero is contained in the 95% CI. Since the latter is more ‘popular’, we employed it in the following.

The script is for a one-compartment model, first order absorption and elimination with an optional lag-time. Flip-flop PK (\(\small{k_{01}=k_{10}}\)) is implemented. Cave: 288 LOC.

sim.ss <- function(nsims, f, D, V, t12.a, t12.e, tlag, tau, fact,
                   CV.lo = 0.1, CV.hi = 0.4, delta = 0.2,
                   method = "linear", dw = FALSE, setseed = TRUE,
                   progress = FALSE, plots = FALSE) {
  suppressMessages(library(car)) # for Durbin-Watson
  # Cave: Long run time if nsims > 1e4
  #       Extremely long run time if dw = TRUE
  if (nsims < 1)      stop("nsims must be at least 1.")
  if (f <= 0 | f > 1) stop("f must be positive and <= 1.")
  if (D <= 0)         stop("D must be positive.")
  if (t12.a <= 0)     stop("t12.a (absorption half life) must be positive.")
  if (t12.e <= 0)     stop("t12.e (elimination half life) must be positive.")
  if (tlag < 0)       stop("tlag must >= 0.")
  if (tau < 0)        stop("tau must be positive.")
  if (CV.lo <= 0)     stop("CV.lo must be positive.")
  if (CV.hi <= 0)     stop("CV.hi must be positive.")
  if (delta < 0)      stop("delta must be positive.")
  if (fact < 3) {
    message("fact < 3 not reasonable. Forced to 3.")
    fact <- 3
  }
  if (!method %in% c("linear", "power")) stop("method not supported.")

  one.comp <- function(f, D, V, k01, k10, tlag, t) {
    # one-compartment model, first order absorption
    # and elimination
    if (!isTRUE(all.equal(k01, k10))) {# common: k01 != k10
      C    <- f * D * k01 / (V * (k01 - k10)) *
              (exp(-k10 * (t - tlag)) - exp(-k01 * (t - tlag)))
      tmax <- log(k01 / k10) / (k01 - k10) + tlag
      Cmax <- f * D * k01 / (V * (k01 - k10)) *
              (exp(-k10 * tmax) - exp(-k01 * tmax))
  }else {                           # flip-flop
      k    <- k10
      C    <- f * D / V * k * (t - tlag) * exp(-k * (t - tlag))
      tmax <- 1 / k
      Cmax <- f * D / V * k * tmax * exp(-k * tmax)
    }
    C[C <= 0] <- 0 # correct negatives caused by lag-time
    res <- list(C = C, Cmax = Cmax, tmax = tmax)
    return(res)
  }

  round.up <- function(x, y) {
    # round x up to the next multiple of y
    return(as.integer(y * (x %/% y + as.logical(x %% y))))
  }

  superposition <- function(f, D, V, k01, k10, tlag, tau, t, n.D) {
    # single dose
    x         <- data.frame(adm = 1L, dose = NA, t = t,
                            C = one.comp(f, D, V, k01, k10, tlag, t)$C)
    x$dose[1] <- D
    # multiple dose (nonparametric superposition)
    y <- x
    for (i in 1:nrow(y)) {
      if (isTRUE(all.equal(y$t[i] %% tau, 0))) y$dose[i] <- D
    }
    delta <- unique(diff(y$t))
    adm   <- seq(0, max(y$t), tau) # administration times
    adm.n <- 1:length(adm)         # number of administration
    if (isTRUE(all.equal(delta, tau))) {# for pre-dose concentrations
      y$adm <- adm.n
  }else {                            # for build-up
      t.tau <- sum(y$t < tau)
      y$adm <- head(rep(adm.n, each = t.tau), nrow(y))
    }
    z <- 2:(n.D + 1)
    for (i in seq_along(z)) {
      tmp.y              <- y$C[y$adm >= z[i]]
      tmp.x              <- x$C[1:length(tmp.y)]
      y$C[y$adm >= z[i]] <- tmp.y + tmp.x
    }
    res <- list(SD = x, MD = y)
    return(res)
  }

  CV2var <- function(CV) {
    # convert CV to variance
    return(log(CV^2 + 1))
  }

  mods <- function(x, method = "linear", dw = FALSE, delta = 0.2) {
    # last three pre-dose concentrations
    x <- tail(x, 3)
    # model and 95% CI of slope
    if (method == "linear") {
      m  <- lm(C ~ t, data = x)
      a  <- coef(m)[[1]]
      ci <- confint(m, "t", level = 0.95)
      if (dw) dw.p <- dwt(m, reps = 1e4)$p
    }
    if (method == "power") {
      m  <- lm(log(C) ~ log(t), data = x)
      a  <- exp(coef(m)[[1]])
      ci <- confint(m, "log(t)", level = 0.95)
      if (dw) dw.p <- dwt(m, reps = 1e4)$p
    }
    # assess whether the CI contains zero
    ifelse (sign(ci[1]) == sign(ci[2]), sig1 <- TRUE, sig1 <- FALSE)
    # assess whether the CI is within +/- delta
    ifelse ((ci[1] >= -delta & ci[2] <= +delta), sig2 <- FALSE, sig2 <- TRUE)
    res <- list(method = method, a = a, b = coef(m)[[2]],
                lwr = ci[[1]], upr = ci[[2]], sig1 = sig1, sig2 = sig2,
                t.start = x$t[1], t.end = x$t[3], dw.p = dw.p)
    return(res)
  }

  geo.mean <- function(x) {
    # geometric mean
    n <- as.integer(unlist(unique(x[1])))
    t <- as.numeric(unlist(unique(x[2])))
    gm <- NA_real_
    for (i in seq_along(t)) {
      y <- x[x[2] == t[i], ]
      if (sum(y$C <= 0) > 0) {
        y <- y[-which(y$C <= 0), ] # only positive values
      }
      if (nrow(y) > 0) {
        gm[i] <- exp(mean(log(y$C), na.rm = TRUE))
      }
    }
    return(gm)
  }

  if (nsims > 1e4) progress <- TRUE  # to keep the spirit high
  k01    <- log(2) / t12.a           # absorption rate constant
  k10    <- log(2) / t12.e           # elimination rate constant
  ifelse (t12.e >= t12.a, t.last <- fact * t12.e, t.last <- fact * t12.a)
  t.last <- round.up(t.last, tau)    # time of last dose
  t      <- seq(0, t.last, tau)      # time of administrations
  n.D    <- length(t)                # number of doses
  if (setseed) set.seed(1234567)     # for reproducibility
  # build-up (high resolution)
  t.bu     <- seq(0, t.last, length.out = (t.last * 10 + 1))
  bu       <- superposition(f, D, V, k01, k10, tlag, tau, t = t.bu, n.D)$MD
  # pre-dose concentrations
  MD       <- bu[!is.na(bu$dose), ]
  mod.a    <- mods(MD, method, dw)$a # intercept without error
  mod.b    <- mods(MD, method, dw)$b # slope/exponent without error
  C.ss     <- tail(MD$C, 1)          # last pre-dose of model
  # simulated pre-dose concentrations
  sim.lo   <- sim.hi <- data.frame(sim = rep(1:nsims, each = n.D),
                                   t = t, C = NA_real_)
  # results with low and high variability
  lo       <- hi <- data.frame(sim = 1:nsims, method = NA_character_,
                               a = NA_real_, b = NA_real_,
                               lwr = NA_real_, upr = NA_real_,
                               sig1 = NA, sig2 = NA,
                               t.start = NA_real_, t.end = NA_real_)
  if (progress) pb <- txtProgressBar(0, 1, 0, char = "\u2588",
                                     width = NA, style = 3)
  for (sim in 1:nsims) {
    # multiplicative error
    C.lo <- rlnorm(n       = n.D,
                   meanlog = log(MD$C) - 0.5 * CV2var(CV.lo),
                   sdlog   = sqrt(CV2var(CV.lo)))
    C.hi <- rlnorm(n       = n.D,
                   meanlog = log(MD$C) - 0.5 * CV2var(CV.hi),
                   sdlog   = sqrt(CV2var(CV.hi)))
    sim.lo$C[sim.lo$sim == sim] <- C.lo
    sim.hi$C[sim.hi$sim == sim] <- C.hi
    # get slope, 95% CI, and significance
    lo[sim, 2:11] <- mods(sim.lo[sim.lo$sim == sim, ], method, dw)
    hi[sim, 2:11] <- mods(sim.hi[sim.lo$sim == sim, ], method, dw)
    if (progress) setTxtProgressBar(pb, sim / nsims)
}# simulations running...
  if (progress) close(pb)
  if (method == "linear") {
    txt <- "slope"
    me  <- "linear model: C ~ a + b * t"
}else {
    txt <- "exponent"
    me  <- "power model: C ~ a * t^b"
  }
  me <- paste(me, "\n delta =", delta,
              "(acceptance limits {-delta, +delta})\n")
  # cosmetics
  fm <- paste0("CV = %4.1f%%, median of ", txt, "s = %+.5f")
  fm <- paste0(fm, "\n  %5.1f%% positive ", txt, "s")
  fm <- paste0(fm, "\n  %5.1f%% negative ", txt, "s")
  fm <- paste0(fm, "\n   significant ", txt, "s: %6.2f%% %s")
  fm <- paste0(fm, "\n   failed ", txt, "s     : %6.2f%% %s")
  # results
  txt <- paste(prettyNum(nsims, format = "d", big.mark =","), "simulations,",
               me, txt, "of model =", sprintf("%+.5f", mod.b),
               "\n last model’s pre-dose concentration =", signif(C.ss, 4), "\n",
               sprintf(fm, 100 * CV.lo, median(lo$b),
                       100 * sum(lo$b > 0) / nsims,
                       100 * sum(lo$b < 0) / nsims,
                       100 * sum(lo$sig1, na.rm = TRUE) / nsims,
                       "(95% CI does not contain zero)",
                       100 * sum(lo$sig2, na.rm = TRUE) / nsims,
                       "(95% CI not within acceptance limits)\n"))
  if (dw) {
    txt <- paste(txt, "  autocorrelation   :",
                 sprintf("%5.1f%%",
                 100 * sum(lo$dw.p < 0.05, na.rm = TRUE) / nsims), "\n")
  }
  txt <- paste(txt, sprintf(fm, 100 * CV.hi, median(hi$b),
                            100 * sum(hi$b > 0) / nsims,
                            100 * sum(hi$b < 0) / nsims,
                            100 * sum(hi$sig1, na.rm = TRUE) / nsims,
                            "(95% CI does not contain zero)",
                            100 * sum(hi$sig2, na.rm = TRUE) / nsims,
                            "(95% CI not within acceptance limits)\n"))
  if (dw) {
    txt <- paste(txt, "  autocorrelation   :",
                 sprintf("%5.1f%%",
                 100 * sum(hi$dw.p < 0.05, na.rm = TRUE) / nsims), "\n")
  }
  cat(txt)
  if (plots) {
    # use [Page][up] / [Page][down] keys to switch
    # between plots in the graphics device
    windows(width = 4.5, height = 4.5, record = TRUE)
    op <- par(no.readonly = TRUE)
    par(mar = c(4, 3.9, 0.2, 0.1), cex.axis = 0.9)
    ylim  <- range(bu$C)
    if (method == "power") {
      last.3 <- tail(MD$t, 3)
      t.tmp <- seq(last.3[1], last.3[3], length.out = 251)
    }
    for (sim in 1:nsims) {
      if (sim == 1) {
        plot(bu$t, bu$C, type = "l", ylim = ylim, axes = FALSE,
             xlab = "time (h)", ylab = "concentration  (unit)")
        grid(nx = NA, ny = NULL)
        abline(v = t, lty = 3, col = "lightgrey")
        axis(1, at = t)
        axis(2, las = 1)
        box()
        points(MD$t, MD$C, pch = 15, cex = 1.25)
      }
      points(sim.lo$t[sim.lo$sim == sim], sim.lo$C[sim.lo$sim == sim],
             pch = 21, cex = 1.1, col = "blue", bg = "#A0A0FF")
      if (lo$sig1[sim]) {
        col <- "red"
        lwd <- 2
    }else {
        col <- "#A0A0FF"
        lwd <- 1
      }
      if (method == "linear") {
        lines(x = lo[sim, 9:10],
              y = c(lo$a[sim] + lo$b[sim] * lo$t.start[sim],
                    lo$a[sim] + lo$b[sim] * lo$t.end[sim]), col = col, lwd = lwd)
      }
      if (method == "power") {
        lines(x = t.tmp, y = lo$a[sim] * t.tmp^lo$b[sim], col = col, lwd = lwd)
      }
    }
    points(t, geo.mean(sim.lo), pch = 24, cex = 1.5, bg = "#FF0808BF")
    for (sim in 1:nsims) {
      if (sim == 1) {
        plot(bu$t, bu$C, type = "l", ylim = ylim, axes = FALSE,
             xlab = "time (h)", ylab = "concentration  (unit)")
        grid(nx = NA, ny = NULL)
        abline(v = t, lty = 3, col = "lightgrey")
        axis(1, at = t)
        axis(2, las = 1)
        box()
        points(MD$t, MD$C, pch = 15, cex = 1.25)
      }
      points(sim.hi$t[sim.hi$sim == sim], sim.hi$C[sim.hi$sim == sim],
             pch = 21, cex = 1.1, col = "blue", bg = "#A0A0FF")
      if (hi$sig1[sim]) {
        col <- "red"
        lwd <- 2
    }else {
        col <- "#A0A0FF"
        lwd <- 1
      }
      if (method == "linear") {
        lines(x = hi[sim, 9:10],
              y = c(hi$a[sim] + hi$b[sim] * hi$t.start[sim],
                    hi$a[sim] + hi$b[sim] * hi$t.end[sim]), col = col, lwd = lwd)
      }
      if (method == "power") {
        lines(x = t.tmp, y = hi$a[sim] * t.tmp^hi$b[sim], col = col, lwd = lwd)
      }
    }
    points(t, geo.mean(sim.hi), pch = 24, cex = 1.5, bg = "#FF0808BF")
    par(op)
  }
}

Let us perform some simulations.
We explore low (\(\small{CV=10\%}\)) and high within-subject variability (\(\small{CV=40\%}\)). In a conservative manner we dose for seven half lives.
Additionally to testing whether the CI contains zero, in analogy to BE we test also whether it lies entirely within – arbitrary – acceptance limits of \(\small{\{-\Delta,+\Delta\}}\).

nsims  <- 1e4L # number of simulations
D      <- 200L # dose
f      <- 2/3  # fraction absorbed (BA)
V      <- 3    # volume of distribution
t12.a  <- 0.75 # absorption half life
t12.e  <- 12   # elimination half life
tlag   <- 0.25 # lag time
tau    <- 24L  # dosing interval
fact   <- 7L   # half life factor
CV.lo  <- 0.10 # drug with low variability
CV.hi  <- 0.40 # drug with high variability
delta  <- 0.2  # for acceptance limits = {1 - delta, 1 + delta}
dw     <- TRUE # Assess autocorrelation
sim.ss(nsims, f, D, V, t12.a, t12.e, tlag, tau, fact,
       CV.lo, CV.hi, delta, dw = dw)

# 10,000 simulations, linear model: C ~ a + b * t 
#  delta = 0.2 (acceptance limits {-delta, +delta})
#  slope of model = +0.01957 
#  last model’s pre-dose concentration = 15.97 
#  CV = 10.0%, median of slopes = +0.01977
#    66.6% positive slopes
#    33.4% negative slopes
#    significant slopes:   5.58% (95% CI does not contain zero)
#    failed slopes     :  78.27% (95% CI not within acceptance limits)
#    autocorrelation   :  31.9% 
#  CV = 40.0%, median of slopes = +0.02072
#    54.8% positive slopes
#    45.2% negative slopes
#    significant slopes:   4.18% (95% CI does not contain zero)
#    failed slopes     :  96.75% (95% CI not within acceptance limits)
#    autocorrelation   :  24.2%

With a half life of 12 hours and once daily dosing after five administrations we are close to the true steady state (≈16). Hence, the slope of the build-up based on the model is slightly positive. Only in true steady state the slope would be zero (see above). Even then we would expect 5% significant results by pure chance. That’s the false positive rate (FPR) of any test performed at level \(\small{\alpha=0.05}\).

With low variability we are above the FPR. That’s correct because the true slope is positive indeed.
On the other hand, with high variability the test is not powerful enough to detect the positive slope and we are close to the level of the test.

Note also that with low variability we see more positive slopes than negative ones – which is correct in principle. With high variability everything drowns in noise.

The confidence interval inclusion approach fails with flying colors. Definitely not a good idea.

Let’s explore a smaller subset, which reflects better what we might face in a particular study.

Fig. 1 First 24 of the simulated subjects, \(\small{CV=10\%}\).
Blue lines estimated regressions, red lines with significant slope, red triangles geometric means, black line model of build-up, black squares its pre-dose concentrations.

Fig. 2 First 24 of the simulated subjects, \(\small{CV=40\%}\).
Blue lines estimated regressions, red triangles geometric means, black line model of build-up, black squares its pre-dose concentrations.

Here with low variability three simulated subjects show a significant slope (FPR 12.5%) and with high variability none. When assessing the Durbin–Watson test, we find a significant autocorrelation in six simulated subjects with low variability and in five with high variability. Hence, it confirms that applying a regression is not justified.

Of course, in reality the level of the true steady state is unknown. The geometric means show that the pre-dose concentrations are stable, though with low variability they seem to rise a bit.

A shorter build-up (fewer doses) might be risky. Here planned with five half lives as recommended in the guide­lines (only four doses, last dosing time rounded up from 60 to 72 hours).

10,000 simulations, linear model: C ~ a + b * t 
 delta = 0.2 (acceptance limits {-delta, +delta})
 slope of model = +0.07828 
 last model’s pre-dose concentration = 15.78
 CV = 10.0%, median of slopes = +0.07712
   97.3% positive slopes
    2.7% negative slopes
   significant slopes:   9.46% (95% CI does not contain zero)
   failed slopes     :  85.79% (95% CI not within acceptance limits)
   autocorrelation   :   6.2%
 CV = 40.0%, median of slopes = +0.06995
   68.9% positive slopes
   31.1% negative slopes
   significant slopes:   4.98% (95% CI does not contain zero)
   failed slopes     :  96.26% (95% CI not within acceptance limits)
   autocorrelation   :   4.4%

The positive slopes are detected. With low variability 9.46% is clearly above the level of the test. With high variability we are close to the level of the test.

As before, we assess a small subset.

Fig. 3 First 24 of the simulated subjects, \(\small{CV=10\%}\).
Blue lines estimated regressions, red triangles geometric means, black line model of build-up, black squares its pre-dose concentrations.

Fig. 4 First 24 of the simulated subjects, \(\small{CV=40\%}\).
Blue lines estimated regressions, red line with significant slope, red triangles geometric means, black line model of build-up, black squares its pre-dose concentrations.

With low variability none of the simulated subjects show a significant slope and with high variability one. However, what might a regulatory assessor conclude? With low variability all slopes are positive and with high variability still 19.

An example from my files. A CNS-drug with polymorphic metabolism (\(\small{CV_\textrm{inter}\gg CV_\textrm{intra}}\)). Given the wide range of half lives reported in the literature and observed in our single dose studies, we dosed o.a.d. for 14 days, trans­lating into an average of 15.5 half lives (minimum 5.2, maximum 46). \[\small{\textrm{PK metric (ng/mL): }\bar{x}_\textrm{geom}\;(CV)}\]\[\small{ \begin{array}{lccc} \textrm{Treatment} & C_\textrm{pd} & C_\textrm{min} & C_\tau\\\hline \textrm{Test} & 28.47\;(80\%) & 26.99\;(79\%) & 33.41\;(59\%)\\ \textrm{Reference} & 29.80\;(77\%) & 28.19\;(75\%) & 32.84\;(61\%)\\\hline \end{array}}\]\(\small{C_\textrm{pd}>C_\textrm{min}}\), since in 26 of the 56 profiles we observed a lag time. Recall the accumulation example above. It is common that \(\small{C_\tau>C_\textrm{pd}}\) because we can never attain true steady state.

Although not assessed in the study, all of these PK metrics would pass conventional bioequivalence. \[\small{\begin{array}{lccccc} \textrm{PK metric} & PE & 90\%\;CL_\textrm{lo} & 90\%\;CL_\textrm{hi} & CV_\textrm{intra} & CV_\textrm{inter}\\\hline \phantom{000}C_\textrm{pd} & \phantom{1}95.54\% & 81.29\% & 112.29\% & 39.5\% & 133.3\%\\ \phantom{000}C_\textrm{min} & \phantom{1}95.74\% & 81.82\% & 112.03\% & 38.3\% & 128.2\%\\ \phantom{000}C_\tau & 100.48\% & 90.27\% & 111.85\% & 25.7\% & \phantom{1}94.3\%\\\hline \end{array}}\] Here the results of the first period.

linear model: C ~ a + b * t 
Δ = 0.5 (acceptance limits {–Δ, +Δ}) 
  42.9% positive slopes
  57.1% negative slopes
  significant slopes:  3.57% (95% CI does not contain zero)
  failed slopes     : 82.14% (95% CI not within acceptance limits)

Fig. 6 Spaghetti plot of a CNS-drug in 28 subjects. Red line regression with significant slope,
red triangles geometric means.

The one subject with a significant slope (with a 95% CI of 0.0039 – 1.038 it was a close shave) had low variability. How­ever, at the end of the profile day the concentration was lower than at 312 hours. Hence, it is reasonable to conclude that the subject was in steady state and the apparently rising concentrations were due to pure chance. There was yet another subject with a similar slope. Since its variability was higher, the slope was not sig­ni­ficant.
Note that ≈82% would fail even with an extremely liberal \(\small{\Delta=0.5}\).

According to the protocol we reported only individual data, descriptive statistics,¹⁰ individual and spaghetti plots. All European regulatory assessors were fine with that.

---

    

Interlude: Flip-flop Pharmacokinetics

\(\small{k_{01}}\): absorption rate constant
\(\small{k_{10}}\): elimination rate constant

Even in a simple one-compartment model \(\small{k_{01}}\) and \(\small{k_{10}}\) are not identifiable. Where for IR formulations it is reasonable to assume \(\small{k_{10}\ll k_{01}}\), i.e., that absorption is much faster than elimination, for prolonged release formulations this must not be – and very likely is not – the case. If you are a dinosaur of biopharmaceutics like me, perhaps you remember the slogan of the 1980s »the flatter is better«. A main target of product development is to decrease \(\small{k_{01}}\) (increase the absorption half life). Never plan the study based on the apparent elimination half life of an IR formulation.

Here an example of flip-flop PK with \(\small{k_{01}=k_{10}=\textsf{8 h}}\), designed with the naïve dosing of seven elimination half lives (56 rounded up to 72 hours).

Fig. 5 24 simulated subjects, \(\small{CV=10\%}\). Blue lines estimated regressions, red triangles geometric means, black line model of build-up, black squares its pre-dose concentrations.

That’s a recipe for disaster. All slopes are positive. Though none is significantly different from zero, it would need a blind assessor to accept your claim that steady state was achieved. And they would be absolutely right. It would need dosing for ten half lives to end up with something pretty similar to Fig. 1.
As usual: Know the drug and know your formulation!

---

Power model

Given the behavior of the build-up, which is concave, one might be tempted set up a power model \[Y=a\cdot X^{\,b}+\epsilon\small{\textsf{.}}\tag{3}\] Like in dose-proportionality²⁰ for convenience its linearized form \[\log_{e}Y=\log_{e}a+b\cdot\log_{e}X+\epsilon \tag{4}\] can be fitted. The tests are the same like in the linear model \(\small{(2)}\).

nsims  <- 1e4L # number of simulations
D      <- 200L # dose
f      <- 2/3  # fraction absorbed (BA)
V      <- 3    # volume of distribution
t12.a  <- 0.75 # absorption half life
t12.e  <- 12   # elimination half life
tlag   <- 0.25 # lag time
tau    <- 24L  # dosing interval
fact   <- 5L   # half life factor
CV.lo  <- 0.10 # drug with low variability
CV.hi  <- 0.40 # drug with high variability
delta  <- 0.2  # for acceptance limits = {1 - delta, 1 + delta}
dw     <- TRUE # Assess autocorrelation
method <- "power"
sim.ss(nsims, f, D, V, t12.a, t12.e, tlag, tau, fact,
       CV.lo, CV.hi, dw = dw, method = method)

# 10,000 simulations, power model: C ~ a * t^b 
#  delta = 0.2 (acceptance limits {-delta, +delta})
#  exponent of model = +0.25554 
#  last model’s pre-dose concentration = 15.78 
#  CV = 10.0%, median of exponents = +0.25334
#    97.8% positive exponents
#     2.2% negative exponents
#    significant exponents:  11.55% (95% CI does not contain zero)
#    failed exponents     :  98.88% (95% CI not within acceptance limits)
#    autocorrelation   :  12.2% 
#  CV = 40.0%, median of exponents = +0.25455
#    69.9% positive exponents
#    30.1% negative exponents
#    significant exponents:   5.23% (95% CI does not contain zero)
#    failed exponents     :  99.59% (95% CI not within acceptance limits)
#    autocorrelation   :  10.8%

Do we gain something? Not really.

As above the subset of 24 subjects.

Fig. 7 First 24 of the simulated subjects, \(\small{CV=10\%}\).
Blue lines estimated regressions, red triangles geometric means, black line model of build-up, black squares its pre-dose concentrations.

Fig. 8 First 24 of the simulated subjects, \(\small{CV=40\%}\).
Blue lines estimated regressions, red line with significant slope, red triangles geometric means, black line model of build-up, black squares its pre-dose concentrations.

We would arrive at the same conclusion than with the linear model. When we run both models we see that in some cases the power model fits better (based on an ANOVA or the AIC). However, we see also that some fits are convex and that is from a theoretical perspective simply wrong. That’s a dead end.

Alternative?

So why are we not fitting the underlying function \[C=C_\textrm{ss}\big{(}1-\exp(-\lambda \times t)\big{)} \tag{5}\] of the build-up directly? There are two problems, a technical and a theoretical one.

We cannot linearize \(\small{(5)}\) like we did with the power model \(\small{(3)}\). Hence, we require nonlinear modeling as provided in R’s function nls(), which requires suitable start values. In many cases \(\small{\max (C) \Rightarrow C_\textrm{ss}}\) and \(\small{\log_{e}2\,/\,t_{1/2}\Rightarrow\lambda}\) work. Alternatively we can set up the model in Phoenix Win­Nonlin. However, in most cases calculating a CI fails.

Even if we obtain a CI, which parameter should we test and against what?
The true steady state is unknown and therefore, testing \(\small{\widehat{C_\textrm{ss}}}\) is not an option. We face a similar problem with testing \(\small{\widehat{\lambda}}\). We have only the half life used in study planning and no clue what the subjects’ true predominant rate constants are.