Introduction

If this article is perceived as overly focused on statistics, I apologize. This is due to my professional background, which has led me to be less skilled at crafting engaging narratives.

‘Bioavailability’ (a portmanteau of ‘biologic availability’) in its current meaning was coined in 1973¹ and ‘Bio­equi­va­lence’ saw the light of day in 1975.²

The site of action (i.e., a receptor) is inaccessible. There should be no space for believes in science. The best definition of bioequivalence is given by the ICH.³

“Two drug products containing the same drug substance(s) are considered bioequivalent if their relative bio­availability (BA) (rate and extent of drug absorption) after administration in the same molar dose lies with­in acceptable predefined limits. These limits are set to ensure com­par­able in vivo performance, i.e., si­mi­la­ri­ty in terms of safety and efficacy.

We will use a simple example in the following. A two-treatment two-sequence two-period (2×2×2) crossover design, subjects 1–6 were in sequence \(\small{\text{TR}}\) and subjects 7–12 in sequence \(\small{\text{RT}}\). \[\small{\begin{array}{ccc} \textsf{Table I}\phantom{0}\\ \text{subject} & \text{T} & \text{R}\\\hline \phantom{1}1 & 71 & 81\\ \phantom{1}2 & 61 & 65\\ \phantom{1}3 & 80 & 94\\ \phantom{1}4 & 66 & 74\\ \phantom{1}5 & 94 & 54\\ \phantom{1}6 & 97 & 63\\ \phantom{1}7 & 70 & 85\\ \phantom{1}8 & 76 & 90\\ \phantom{1}9 & 54 & 53\\ 10 & 99 & 56\\ 11 & 83 & 90\\ 12 & 51 & 68\\\hline \end{array}}\]

top of section ↩︎

The 1970s

Problems were reported with formulations of Narrow Therapeutic Index Drugs (NTIDs) like phenytoin,^{4 5 6 7} digoxin,^{1 8 9} warfarin,¹⁰ theophylline,¹¹ primidone.¹² Some show nonlinear pharmacokinetics (phenytoin) or are auto-inducers (warfarin).

• Excipient changed from CaS0₄ to lactose^{5 6}
• The API was altered (e.g., particle size,^{7 9} amorphous to crystalline¹⁰)
• Variable disintegration time
• Dissolution testing not mandatory
• No in vivo studies were performed comparing the new to the approved formulation
• Breakthrough-seizures⁴ and intoxications^{5 6} (phenytoin) and variable or poor effect (digoxin, theophylline)

Generic drugs in the current sense did not yet exist at that time; only the content had to meet the USP requirements.

“Although in 1969 Professor John Wagner demonstrated to the Bureau of Medicine, methods for comparing areas under the serum versus time curve (AUC) to estimate bioequivalence, his approach was ignored inasmuch as the FDA hierarchy did not believe a problem existed, and there­fore such studies would not be ne­cessary. For their part the Offices of Pharmaceutical Re­search and Compliance in the Bureau of Medicine and the Com­mis­sio­ner’s Office believed that the “Bioavailability Problem” as some called it was a “Content Uni­formity Problem”.¹³ In 1971 for example, when notified of a “Bioavailability Problem” with a generic di­goxin product, FDA in­vestigated and ascertained that one manufacturer first added all the excipients into a 55-gal drum, then added di­gox­in, closed the lid, and mixed it by rolling the drum across the floor a few times. The content uniformity of those tablets varied from 10% to 156%.

Following a ‘Conference on Bioavailability of Drugs’ held at the National Academy of Sciences of the United States in 1971, a guideline was published the following year.¹⁵

“[…] the mean of AUC of the generic had to be within 20% of the mean AUC of the approved product. At first this was determined by using serum versus time plots on specially weighted paper, cutting the plot out and then weighing each separately.

top of section ↩︎ previous section ↩︎

example          <- data.frame(subject   = rep(1:12, each = 2),
                               sequence  = c(rep("TR", 12), rep("RT", 12)),
                               treatment = c(rep(c("T", "R"), 6),
                                             rep(c("R", "T"), 6)),
                               period    = rep(1:2, 12),
                               Y         = c(71, 81, 61, 65, 80, 94,
                                             66, 74, 94, 54, 97, 63,
                                             85, 70, 90, 76, 54, 53,
                                             56, 99, 90, 83, 68, 51))
factors          <- c("subject", "period", "treatment")
example[factors] <- lapply(example[factors], factor) # factorize the data
# additive model (untransformed data, differences); sequence not in the model!
muddle           <- lm(Y ~ subject + period + treatment, data = example)
CI               <- as.numeric(confint(muddle, level = 0.95)["treatmentT", ])
PE               <- coef(muddle)[["treatmentT"]]
# Percentages (flawed!)
mean.T           <- mean(example$Y[example$treatment == "T"])
mean.R           <- mean(example$Y[example$treatment == "R"])
PE.pct           <- 100 * mean.T / mean.R
CI.pct           <- 100 * (CI + mean.R) / mean.R
result           <- data.frame(method = c("differences", "percentages"),
                               PE = c(sprintf("%+.3f", PE),
                                      sprintf("%6.2f%%",  PE.pct)),
                               lower = c(sprintf("%+.3f", CI[1]),
                                         sprintf("%.2f%%",  CI.pct[1])),
                               upper = c(sprintf("%+.3f", CI[2]),
                                         sprintf("%6.2f%%",  CI.pct[2])),
                               BE = c("", "fail"))
if (CI.pct[1] >= 80 & CI.pct[2] <= 120) result$BE[2] <- "pass"
names(result)[3:4] <- c("lower CL", "upper CL")
print(result, row.names = FALSE)

#       method      PE lower CL upper CL   BE
#  differences  +2.250  -12.807  +17.307     
#  percentages 103.09%   82.42%  123.76% fail

The Roaring 1980s

The generic boom started 1984 in the U.S. with the ‘Drug Price Competition and Patent Term Restoration Act’ (informally known as ‘Hatch-Waxman Act’).²³

The approval process was different for innovator (originator) and generic companies.

Innovators:

• Preclinical data
• Documentation of pharmaceutical quality
• In clinical phase I documentation of pharmacokinetics (PK) in healthy subjects, dose finding, safety / tolerability, food effect
• In phase II efficacy & safety in a small groups of patients
• In phase III demonstration of efficacy & safety versus placebo in well-powered studies:
  Efficacy: Non-Inferiority/Superiority
  Safety: Non-Superiority

Generic companies:

• Documentation of pharmaceutical quality
• Not required:
  • Any in vivo study
  • Sometimes comparison of disintegration, rarely comparison of dissolution was performed

Regulatory concerns about generic substitution arose, leading to extensive discussions which method could be used to compare formulations.

• Pharmaceutical equivalence
• Bioequivalence (BE)
• Therapeutic equivalence

There was an early agreement that pharmaceutical equivalence is too permissive and therapeutic equivalence would require extremely large studies in patients.²⁴ Hence, comparing the bioavailability (BA) in healthy volunteers seemed to be a reasonable compromise.¹⁷

“What is the justification for studying bioequivalence in healthy volunteers?
“Variability is the enemy of therapeutics” and is also the enemy of bioequivalence. We are trying to determine if two dosage forms of the same drug behave similarly. Therefore we want to keep any other variability not due to the dosage forms at a minimum. We choose the least vari­able “test tube”, that is, a healthy vo­lun­teer.
Disease states can definitely change bioavailability, but we are testing for bioequivalence, not bio­avail­ability.

Whereas in pharmacokinetics (PK) by bioavailability exclusively the Area under Curve extrapolated to infinite time (\(\small{AUC_{0-\infty}}\)) is meant, the FDA introduced two new terms, namely

1. the ‘rate of bioavailability’ – measured by the maximum concentration (\(\small{C_\text{max}}\)) and
2. the ‘extent of bioavailability’ – measured by the \(\small{AUC}\).

The former is understood as a surrogate for the absorption rate \(\small{k\,_\text{a}}\) in a PK model. I pre­fer – like the ICH³ and the FDA since 2003²⁶ – rate and extent of absorption, in order not to contaminate the original meaning of BA in PK.

Let us consider the basic equation of pharmacokinetics \[\frac{f\cdot D}{CL}=\frac{f\cdot D}{V\cdot k_\text{ el}}=AUC_{0-\infty}=\int_{0}^{\infty}C(t)\,dt,\tag{3}\] where \(\small{f}\) is the fraction absorbed (we are interested in the comparison of formulations), \(\small{D}\) is the dose, \(\small{CL}\) is the clear­ance, \(\small{V}\) is the apparent volume of distribution, \(\small{k\,_\text{el}}\) is the elimination rate constant, and \(\small{C(t)}\) is the plasma concentration with time. We see immediately that for identical²⁷ doses and invariate²⁸ \(\small{CL}\), \(\small{V}\), \(\small{k\,_\text{el}}\) (which are drug-spe­ci­fic), com­paring the \(\small{AUC}\text{s}\) allows to compare the frac­tions absorbed.

“Pharmacokinetics: one of the magic arts of divination whereby needles are stuck into dummies in an attempt to predict profits.

It must be mentioned that \(\small{C_\text{max}}\) is not sensitive to even substantial changes in the rate of absorption \(\small{k\,_\text{a}}\), since it is a composite metric.²⁹ In a one compartment model it depends on \(\small{k\,_\text{a}}\), \(\small{f}\) and both the elimination rate con­stant \(\small{k\,_\text{el}}\) and \(\small{V}\) (or \(\small{CL}\) if you belong to the other church).³⁰ Whereas \(\small{k\,_\text{a}}\) and \(\small{f}\) are properties of the formulation – we are interested in – the others are properties of the drug. \[\eqalign{ t_\textrm{max}&=\frac{\log_{e}(k\,_\text{a}/k\,_\text{el})}{k\,_\text{a}-k\,_\text{el}}\\ C_\textrm{max}&=\frac{f\cdot D\cdot k\,_\text{a}}{V\cdot (k\,_\text{a}-k\,_\text{el})}\large(\small\exp(-k\,_\text{el}\cdot t_\textrm{max})-\exp(-k\,_\text{a}\cdot t_\textrm{max})\large)\tag{4}}\] Therefore, when using it as a surrogate for the absorption rate one must keep in mind that formulations with different fractions absorbed and \(\small{t_\text{max}}\) might show the same \(\small{C_\text{max}}\).
It took ten years before the alternative metric \(\small{C_\text{max}/AUC}\) (based on theo­re­tical considerations and simulations) was proposed.^{31 32 33} Apart from being independent from \(\small{f}\), it is substantially less variable than \(\small{C_\text{max}}\). Regrett­ably, it was never implemented in any guideline.

In the early 1980s originators failed in trying to falsify the concept (i.e., comparing BE in healthy volunteers to large the­ra­peu­tic tqui­va­lence (TE) studies in patients): If BE passed, TE passed as well and vice versa. If they would have succeeded (BE passed while TE failed), generic companies would have to demonstrate TE in order to get pro­ducts approved. Such studies would have to be much larger than the originators’ phase III studies, making them economically infeasible.²⁴ Essentially, that would have meant an early end of the young generic industry.

However, comparative BA is also used by originators in scale-up of formulations used in phase III to the to-be-mar­keted formulation, supporting post-approval changes, in line extensions of approved products, and for testing of drug-drug interactions or food effects. Hence, a substantial part of BE trials are performed by originators. If they had been successful to refute the concept, they would have shot into their own foot.

In the mid 1980s a consensus was reached, i.e., that generic approval should only be acceptable after suitable in vivo equivalence.

The main assumption in BE was (and still is) that ‘similar’ plasma concentrations in healthy volunteers will lead to similar concentrations at the target site (i.e., a receptor) and thus, to similar effects in patients. It was still an open issue whether BE should be interpreted as a surrogate of clinical efficacy/safety or a measure of pharmaceutical quality. Where­as in the 1980s the former was prevalent, since the 1990s the latter is mainstream.
A somewhat naïve interpretation of the PK metrics is that \(\small{AUC}\) directly translates to efficacy and \(\small{C_\text{max}}\) to safety. Especially the latter is not correct because any difference in \(\small{C_\text{max}}\) leads to a relatively smaller difference in the ma­xi­mum effect \(\small{E_\text{max}}\).

There was no consensus about the definition of ‘similarity’ and the statistical methodology to compare plasma profiles. Two early methods are outlined in the following.

top of section ↩︎ previous section ↩︎

example       <- data.frame(subject   = rep(1:12, each = 2),
                            sequence  = c(rep("TR", 12), rep("RT", 12)),
                            treatment = c(rep(c("T", "R"), 6),
                                          rep(c("R", "T"), 6)),
                            period    = rep(1:2, 12),
                            Y         = c(71, 81, 61, 65, 80, 94,
                                          66, 74, 94, 54, 97, 63,
                                          85, 70, 90, 76, 54, 53,
                                          56, 99, 90, 83, 68, 51))
rule.75.75    <- reshape(example, idvar = "subject", timevar = "treatment",
                         drop = c("sequence", "period"), direction = "wide")
names(rule.75.75)[2:3] <- c("T", "R")
rule.75.75$T.R <- 100 * (rule.75.75$T / rule.75.75$R)
for (i in 1:nrow(rule.75.75)) {
  if (rule.75.75$T.R[i] >= 75 & rule.75.75$T.R[i] <= 125) {
    rule.75.75$BE[i]     <- TRUE
    rule.75.75$within[i] <- "yes"
  } else {
    rule.75.75$BE[i]     <- FALSE
    rule.75.75$within[i] <- "no"
  }
}
names(rule.75.75)[c(4, 6)] <- c("T/R (%)", "±25%")
BE            <- "Failed BE by the"
if (sum(rule.75.75$BE) / nrow(rule.75.75) >= 0.75) BE <- "Passed BE by the"
print(rule.75.75[, c(1:4, 6)], row.names = FALSE); cat(BE, "75/75 Rule.\n")

#  subject  T  R   T/R (%) ±25%
#        1 71 81  87.65432  yes
#        2 61 65  93.84615  yes
#        3 80 94  85.10638  yes
#        4 66 74  89.18919  yes
#        5 94 54 174.07407   no
#        6 97 63 153.96825   no
#        7 70 85  82.35294  yes
#        8 76 90  84.44444  yes
#        9 53 54  98.14815  yes
#       10 99 56 176.78571   no
#       11 83 90  92.22222  yes
#       12 51 68  75.00000  yes
# Passed BE by the 75/75 Rule.

t-test

Another suggestion was testing for a statistically significant difference at level \(\small{\alpha=0.05}\). The null hypothesis was that formulations are equal (\(\small{\mu_\text{T}-\mu_\text{R}=0}\)).

Let’s assess our example in R again:

example        <- data.frame(subject   = rep(1:12, each = 2),
                             sequence  = c(rep("TR", 12), rep("RT", 12)),
                             treatment = c(rep(c("T", "R"), 6),
                                           rep(c("R", "T"), 6)),
                             period    = rep(1:2, 12),
                             Y         = c(71, 81, 61, 65, 80, 94,
                                           66, 74, 94, 54, 97, 63,
                                           85, 70, 90, 76, 54, 53,
                                           56, 99, 90, 83, 68, 51))
tt             <- reshape(example, idvar = "subject", timevar = "treatment",
                          drop = c("sequence", "period"), direction = "wide")
tt$T.R         <- tt[, 2] - tt[, 3]
names(tt)[2:4] <- c("T", "R", "T–R")
p              <- t.test(x = tt$T, y = tt$R, paired = TRUE)$p.value
BE             <- "Failed BE"
if (p >= 0.05) BE <- "Passed BE"
print(tt, row.names = FALSE); cat(sprintf("%s by a paired t-test (p = %.4f).\n", BE, p))

#  subject  T  R T–R
#        1 71 81 -10
#        2 61 65  -4
#        3 80 94 -14
#        4 66 74  -8
#        5 94 54  40
#        6 97 63  34
#        7 70 85 -15
#        8 76 90 -14
#        9 53 54  -1
#       10 99 56  43
#       11 83 90  -7
#       12 51 68 -17
# Passed BE by a paired t-test (p = 0.7381).

We calculate a \(\small{p}\)-value of 0.7381, which is statistically not significant (\(\small{\geq\alpha}\)) and the study would pass again.

However, we face a similar problem like with the 75/75 Rule. If the differences show high variability, the study would pass. On the other hand, if there is low variability in the differences, the study would fail. This is counterintuitive and actually the opposite of what regulators want.

    

Interlude 1

One of my early sins³⁷ – it was not the last…
After phenytoin intoxications in Austria³⁸ we compared three ge­ne­rics (containing the free acid like the ori­gi­nator, Na-, or Ca-salt) to the reference in a cross­over design. All formulations have been approved and were marketed in Austria. Although at that time I already calculated a 95% CI, the reviewers of our manuscript insisted in testing for a significant difference ‘because it is state of the art’.

Fig. 1 Phenytoin 3 × 100 mg equivalent, single dose fasting.

Two generics were statistically significant different from the reference (\(\small{\text{T}_1}\) containing the free acid like the originator and \(\small{\text{T}_3}\) containing the Ca-salt). \(\small{\text{T}_2}\) containing the Na-salt was statistically not significant different and, thus, considered equi­va­lent – despite its high \(\small{\text{T}/\text{R}\textsf{-}}\)ratio (Table II). \[\small{ \begin{array}{ccccc} \textsf{Table II}\phantom{0000}\\ \text{formulation} & \text{T}/\text{R (%)} & p & & \text{BE}\\\hline \text{T}_1 & 146 & 0.0195\phantom{6} & \text{*} & \text{fail}\\ \text{T}_2 & 134 & 0.151\phantom{96} & \text{n.s.} & \text{pass}\\ \text{T}_3 & \phantom{1}28 & 0.00596 & \text{**} & \text{fail}\\\hline \end{array}}\] If we would evaluate the study according to current standards (i.e., by the 90% CI inclusion approach based on \(\small{\log_{e}\textsf{-}}\)trans­formed data and acceptance limits of 80.00–125.00%), all generics would fail. \(\small{\text{T}_3}\) would even be bio­in­equi­valent because its upper CL is way below 80% (Table III).
If we would adjust for multiplicity (\(\small{\alpha_\text{adj}=0.05/3=0.1\dot{6}\mapsto 96.6\dot{6}\text{% CI}}\)) – although not required in an explo­ra­tory study – the outcome would be even worse (Table IV). \[\small{\begin{array}{ccccc} \textsf{Table III}\phantom{0000}\\ \text{formulation} & \text{PE (%)} & \text{CL}_\text{lower}\text{(%)} & \text{CL}_\text{upper}\text{ (%)} & \text{BE}\\\hline \text{T}_1 & 151.12 & 118.75 & 192.32 & \text{fail (inconclusive)}\\ \text{T}_2 & 139.39 & \phantom{1}95.91 & 202.60 & \text{fail (inconclusive)}\\ \text{T}_3 & \phantom{1}21.67 & \phantom{1}10.25 & \phantom{2}45.81 & \text{fail (inequivalent)}\\\hline \end{array}}\] \[\small{\begin{array}{ccccc} \textsf{Table IV}\phantom{0000}\\ \text{formulation} & \text{PE (%)} & \text{CL}_\text{lower}\text{(%)} & \text{CL}_\text{upper}\text{ (%)} & \text{BE}\\\hline \text{T}_1 & 151.12 & 106.67 & 214.09 & \text{fail (inconclusive)}\\ \text{T}_2 & 139.39 & \phantom{1}81.20 & 239.28 & \text{fail (inconclusive)}\\ \text{T}_3 & \phantom{1}21.67 & \phantom{10}7.34 & \phantom{2}63.93 & \text{fail (inequivalent)}\\\hline \end{array}}\] Given the nonlinear PK of phenytoin,^{39 40} switching a patient from the originator to the generics with high \(\small{\text{T}/\text{R}\textsf{-}}\)ratios would be problematic – potentially leading to toxicity after multiple doses. Even worse would be switching from the ge­ne­ric \(\small{\text{T}_3}\) with its low \(\small{\text{T}/\text{R}\textsf{-}}\)ratio to any of the other formulations.

top of section ↩︎ previous section ↩︎

    

ANOVA and beyond

An Analysis of Variance (ANOVA) instead of a t-test allows to take period-effects into account.^{41 42 43} This decade was also the heyday of Bayesian methods.^{44 45 46 47} Nomograms for sample size estimation were also Bayesian⁴⁸ but happily mis­used by frequentists. New parametric^{49 50} as well as nonparametric methods entered the stage.^{50 51} Metrics to compare controlled release formulations in steady state were proposed.^{52 53 54} The first software to evaluate 2×2×2 crossover studies was released in the public domain.⁵⁵

    

The acceptance range in bioequivalence is based on a ‘clinically relevant difference’ \(\small{\Delta}\), i.e., for data following a lognormal dis­tri­bu­tion \[\left\{\theta_1,\theta_2\right\}=\left\{100\,(1-\Delta),100\,(1-\Delta)^{-1}\right\}\tag{5}\] It must be mentioned that the commonly assumed \(\small{\Delta=20\%}\) leading to \(\small\left\{80.00\%,125.00\%\right\}\)⁵⁶ is arbitrary (as is any other).

    

An important leap forward was the Two One-Sided Tests Procedure (TOST)¹⁶ – al­though it was never implemented in its original form \(\small{(6)}\) in regulatory practice. In­stead, the confidence interval inclusion approach \(\small{(7)}\) made it to the guidelines. Al­though these approaches are operationally identical (i.e., their outcomes [pass | fail] are the same), these are statistically different methods:

1. The TOST Procedure gives two \(\small{p}\)-values, namely \(\small{p(\theta_0\geq\theta_1)}\) and \(\small{p(\theta_0\leq\theta_2)}\). BE is concluded if both \(\small{p}\)-values are \(\small{\leq\alpha}\).

\[\begin{matrix}\tag{6} H_\textrm{0L}:\frac{\mu_\textrm{T}}{\mu_\textrm{R}}\leq\theta_1\:vs\:H_\textrm{1L}:\frac{\mu_\textrm{T}}{\mu_\textrm{R}}>\theta_1\\ H_\textrm{0U}:\frac{\mu_\textrm{T}}{\mu_\textrm{R}}\geq\theta_2\:vs\:H_\textrm{1U}:\frac{\mu_\textrm{T}}{\mu_\textrm{R}}<\theta_2 \end{matrix}\]

2. In the CI inclusion approach BE is concluded if the two-sided \(\small{1-2\,\alpha}\) CI lies entirely within the acceptance range \(\small{\left\{\theta_1,\theta_2\right\}}\). For an explanation why a 90% CI (and not a 95% CI like in phase III) is used, see another article. \[H_0:\frac{\mu_\textrm{T}}{\mu_\textrm{R}}\ni\left\{\theta_1,\theta_2\right\}\:vs\:H_1:\theta_1<\frac{\mu_\textrm{T}}{\mu_\textrm{R}}<\theta_2\tag{7}\]

When we evaluate our example by \(\small{(6)}\), we get \(\small{p(\theta_0\geq\theta_1)=0.0155}\) and \(\small{p(\theta_0\leq\theta_2)=0.0515}\). Since one of the \(\small{p\textsf{-}}\)values is \(\small{>\alpha}\), the study would fail.

    

Interlude 2

It is a misconception that a certain CI of a sample (i.e., a particular study) contains the – true but unknown – population mean \(\small{\mu}\) with \(\small{1-\alpha}\) probabilty. Let’s simulate some studies and evaluate them by \(\small{(7)}\):

invisible(library(PowerTOST))
set.seed(123) # for reproducibility of simulations
mue      <- 1 # true population mean
CV       <- 0.25
studies  <- 100
x        <- sampleN.TOST(CV = CV, theta0 = mue, targetpower = 0.8, print = FALSE)
subjects <- x[["Sample size"]]
power    <- x[["Achieved power"]]
# simulate subjects within studies, lognormal distribution
samples  <- data.frame(study     = rep(1:studies, each = subjects * 2),
                       subject   = rep(rep(1:subjects, studies), each = 2),
                       period    = rep(rep(1:2, studies), 2),
                       sequence  = rep(c(rep(c("TR"), subjects),
                                         rep(c("RT"), subjects)), studies),
                       treatment = c(rep(c("T", "R"), subjects / 2),
                                     rep(c("R", "T"), subjects / 2)),
                       Y         = rlnorm(n = subjects * studies * 2,
                                          meanlog = log(mue) - 0.5 * log(CV^2 + 1),
                                          sdlog = sqrt(log(CV^2 + 1))))
facs     <- c("subject", "period", "treatment")
samples[facs] <- lapply(samples[facs], factor) # factorize the data
result   <- data.frame(study = 1:studies, PE = NA_real_,
                       lower = NA_real_, upper = NA_real_,
                       BE = FALSE, contain = TRUE)
grand.PE <- numeric(studies)
for (i in 1:studies) {
  temp           <- samples[samples$study == i, ]
  heretic        <- lm(log(Y) ~ period + subject + treatment, data = temp)
  result$PE[i]   <- 100 * exp(coef(heretic)[["treatmentT"]])
  result[i, 3:4] <- 100 * exp(confint(heretic, level = 0.90)["treatmentT", ])
  if (round(result[i, 3], 2) >= 80 & round(result[i, 4], 2) <= 125)
    result$BE[i] <- TRUE
  if (result$lower[i] > 100 * mue | result$upper[i] < 100 * mue) result$contain[i] <- FALSE
  grand.PE[i]    <- mean(result$PE[1:i]) # (cumulative) grand means
}
dev.new(width = 4.5, height = 4.5)
op       <- par(no.readonly = TRUE)
par(mar = c(3.05, 2.9, 1.4, 0.75), cex.axis = 0.9, mgp = c(2, 0.5, 0))
xlim     <- range(c(min(result$lower), 1e4 / min(result$lower),
                    max(result$upper), 1e4 / max(result$upper)))
plot(1:2, 100 * rep(mue, 2), type = "n", log = "x", xlab = "PE [90% CI]",
     ylab = "study  #", axes = FALSE,
     xlim = xlim, ylim = range(result$study))
abline(v = 100 * c(0.8, mue, 1.25), lty = c(2, 1, 2))
axis(1, at = c(125, pretty(xlim)),
     labels = sprintf("%.0f%%", c(125, pretty(xlim))))
axis(2, at = c(1, pretty(1:studies)[-1]), las = 1)
axis(3, at = 100 * mue, label = expression(mu))
box()
lines(grand.PE, 1:studies, lwd = 2)
for (i in 1:studies) {
  if (result$BE[i]) {       # pass
    clr <- "blue"
  } else {                  # fail
    if (result$contain[i]) {# mue within CI
      clr <- "magenta"
    } else {                # mue not in CI
      clr <- "red"
    }
  }
  lines(c(result$lower[i], result$upper[i]), rep(i, 2), col = clr)
  points(result$PE[i], i, pch = 16, cex = 0.6, col = clr)
}
par(op)

Fig. 2 2×2×2 crossover studies (\(\small{\mu}\) = 100%, \(\small{CV}\) = 25%: \(\small{n}\) = 24 for ≥80% power).

In 7% of studies the population mean \(\small{\mu}\) is not contained in the 90% CI (red lines). In other words, given the result of a single study we can never know where \(\small{\mu}\) lies. Only the grand mean (mean of sample means \(\small{\frac{1}{n}\sum_{i=1}^{i=n}\overline{x_i}}\)) approaches \(\small{\mu}\) for a large num­ber of samples. After the 100^th study it is with 99.44% pretty close to \(\small{\mu}\) (for geeks: The convergence is poor; when simulating 25,000 studies, it is 100.23%). How­ever, nobody would repeat a – passing – study (blue lines) for such a rather un­inter­esting information, right?
This explains also why a particular study might fail by pure chance even if a formulation is equivalent (here 15% of studies; red or magenta lines). Such cases are related to the producer’s risk (Type II Error = 1 – power), which is for the given conditions 16.3%. On the other hand, it is also possible that a formulation which is not equivalent might pass. These cases are related to the patient’s risk (Type I Error).
For details see the articles about hypotheses, treatment effects, post hoc power, and sample size estimation. Science is a cruel mistress.

    

At a hearing in 1986 the FDA confirmed that \(\small{(6)}\) or \(\small{(7)}\) of untransformed data should be used with \(\small{\Delta=20\%}\). If clinically relevant, tighter limits (\(\small{\Delta=10\%}\)) might be needed.⁵⁷

The first German guideline was drafted by the Working Group for Pharmaceutical Pro­cess Engineering (Ar­beits­ge­mein­schaft für Phar­ma­zeu­tische Ver­fah­rens­tech­nik) in 1985.⁵⁸ It was presented and discussed in 1987.^{59 60 61}

In 1988 wider acceptance limits of 70 – 130% were proposed for \(\small{C_\text{max}}\) due to its inherent high variability⁶² (as a one-point metric practically always larger than the one of the integrated metric \(\small{AUC}\)).

The Australian draft guideline was published in 1988.⁶³ It was the first covering not only the design and evaluation but also validation of bioanalytical methods. The model with effects period, subject, treatment^{20 43} was rec­om­mend­ed and a test for se­quence-ef­fects was not considered necessary. The problematic conversion of differences to percentages was acknowledged and Fieller’s CI^{21 22} discussed. Kudos to both!

In 1989 a series of loose-leaf binders was started.⁶⁴ It contained raw-data of generic drugs marketed in Germany, the evaluation provided by companies, as well as results recalculated by the ZL (Central Laboratory of German Phar­ma­cists). Including the 6^th supplement of 1996 it contained more than 2,000 pages… It was an indispensible resource for planning new studies and also showed the ‘journey’ of dossiers (i.e., the same study being used by different companies).

The BioInternational conference series set milestones in the development of testing for bioequivalence. The first in Toronto 1989 dealt with the \(\small{\log_{e}\textsf{-}}\)transformation of data and the definition of highly variable drugs (HVDs).⁶⁵ There was a poll among the participants about the \(\small{\log_{e}\textsf{-}}\)transformation. Out­come: ⅓ never, ⅓ always, ⅓ case by case (i.e., perform both analyses and report the one with narrower CI ‘because it fits the data better’). Let’s be silent about the last team.⁶⁶ HVDs were defined as drugs with intra-subject variabilities of more than 30% but problems might be evident already at 25%.

top of section ↩︎ previous section ↩︎

The Boring (?) 1990s

The original acceptance range was symmetrical around 100%. In \(\small{\log_{e}\textsf{-}}\)scale it should be symmetrical around \(\small{0}\) (because \(\small{\log_{e}1=0}\)). What happens to our \(\small{\Delta}\), which should still be 20%? Due to the positive skewness of the lognormal distribution a lively discussion started after early publications proposing 80 – 125%.^{19 41} Keeping 80 – 120% would have been flawed because the maximum power should be obtained at \(\small{\mu_\text{T}/\mu_\text{R}=1}\) for \[\exp\left((\log_{e}\theta_1+\log_{e}\theta_2)/2\right),\tag{8}\] which works only if \(\small{\theta_2=\theta_1^{-1}}\) or \(\small{\theta_1=\theta_2^{-1}}\). Keeping the original limits, maximum power would be obtained at \(\small{\mu_\text{T}/\mu_\text{R}=\exp((\log_{e}0.8+\log_{e}1.2)/2)\approx0.979796}\).

Fig. 3 Power for a 2×2×2 design and limits 0.80 – 1.20.
Note that the \(\small{\theta}\)-axis is in log-scale.

There were three parties (all agreed that the acceptance range should be symmetrical in \(\small{\log_{e}\textsf{-}}\)scale and consequently asymmetrical when back-transformed). These were their arguments and suggestions:

\[\left\{\theta_1,\theta_2\right\}=\left\{100\,(1-\Delta),100/(1-\Delta)\right\}=80-125\%\tag{11}\]

The 90% CI inclusion approach \(\small{(7)}\) based on \(\small{\log_{e}\textsf{-}}\)transformed data with acceptance limits of 80.00 – 125.00% \(\small{(5)}\) was the winner.

Fig. 4 Power for a 2×2×2 design and limits 0.80 – 1.25.
Note the symmetry: power for any \(\small{1/\theta=\theta}\).

First sample size tables for the multiplicative model with the acceptance range 80 – 125% were published⁶⁷ and ex­tended for narrower (90 – 111%) and wider (70 – 143%) acceptance ranges.⁶⁸ The nonparametric method was improved taking period-effects into account.^{69 70} Drug-drug and food-in­ter­action studies should be assessed for equi­va­lence.⁷¹ The general applicability of average BE was challenged and the concept of individual and population bioequivalence outlined.^{72 73 74} The first textbook dealing exclusively with BA/BE was published.⁷⁵

This was also the decade of updated and new guidelines. A European draft guidance was published in 1990;⁷⁶ the final guideline was published in December 1991 and came into force in June 1992.⁷⁷ The 90% CI inclusion approach of \(\small{\log_{e}\textsf{-}}\)transformed data with an acceptance range of 80 – 125% was recommended and for NTIDs the acceptance range may need to be tightened. Due to its inherent higher variability a wider acceptance range may be acceptable for \(\small{C_\text{max}}\). If inevitable and clinically acceptable, a wider acceptance range may also be used for \(\small{AUC}\). Only if clinically relevant, a nonparametric analysis of \(\small{t_\text{max}}\) was re­comm­end­ed.
An in vivo stuy was not required if the new formulation is

1. to be parenterally administered as a solution and contains the same API(s) and excipients in the same concentrations as the reference or
2. is a liquid oral form in solution (elixir, syrup, etc.) containing the API(s) in the same concentration and form as the reference, not containing excipients that may significantly affect gastric passage or absorption of the active substance.

Similar statements about solutions were given in all later guidelines. The second lead to application of the Bio­phar­ma­ceu­tic Classi­fi­cation System (BCS).⁷⁸ More about that later.

In July 1992 the first guidance of the FDA was published.⁷⁹ An ANOVA of \(\small{\log_{e}\textsf{-}}\)transformed data was re­com­mend­ed and the nested subject(sequence) term in the statistical model entered the scene. It must be mentioned that in com­pa­rative BA studies subjects are usually uniquely coded. Hence, the term subject(sequence) is a bogus one⁸⁰ and could be replaced by the simple subject as well (see below for an example). Regrettably this model was implemented in all global guidelines ever since.

In the same year the Canadian guidance for Immediate Release (IR) formulations was published.⁸¹ To that time is was the most extensive one because it gave not only the method of evaluation, but information about the study design, sample size, ethics, bioanalytics, etc. It differed from the others in the relaxed requirement for \(\small{C_\text{max}}\), where only the \(\small{\text{T}/\text{R}\textsf{-}}\)ratio has to lie within 80 – 125% (instead of its CI).

In 1998 the World Health Organization published its first guideline,⁸² which was similar to the European one.

Table V shows the result of the example evaluated by the various methods. \[\small{\begin{array}{lcccc} \textsf{Table V}\phantom{0}\\ \phantom{0}\text{Method} & \text{Model} & \text{PE} & \text{power},p,\text{CI} & \text{BE?}\\\hline \text{80/20 Rule} & \text{additive} & - & 47.22\% & \text{fail}\\ \text{TOST} & \text{additive} & +2.250\;(103.09\%) & 0.0155,\,0.0515 & \text{fail}\\ \text{95% CI} & \text{additive} & +2.250\;(103.09\%) & -12.807\,,+17.307\;(82.61-123.76\%) & \text{fail}\\ \text{Westlake} & \text{additive} & \pm0.000\;(100.00\%) & \pm16.143\;(\pm21.48\%) & \text{fail}\\\hline \text{80/20 Rule} & \text{multiplicative} & - & 73.57\% & \text{fail}\\ \text{TOST} & \text{multiplicative} & 102.82\% & 0.0099,\,0.0283 & \text{pass}\\ \text{90% CI} & \text{multiplicative} & 102.82\% & \phantom{1}87.25-121.17\% & \text{pass}\\ \text{Westlake} & \text{multiplicative} & 100.00\% & \pm17.72\% & \text{pass}\\ \text{75/75 Rule} & \text{multiplicative} & - & - & \text{pass}\\\hline \end{array}}\] In the additive model the acceptance range was 80 – 120%, whereas in the multiplicative model it is 80 – 125%. Since in the former differences are assessed – wrong – percentages are given in brackets.

As of today only the 90% CI inclusion approach is globally accepted. Our example in R again:

example       <- data.frame(subject   = rep(1:12, each = 2),
                            sequence  = c(rep("TR", 12), rep("RT", 12)),
                            treatment = c(rep(c("T", "R"), 6),
                                          rep(c("R", "T"), 6)),
                            period    = rep(1:2, 12),
                            Y         = c(71, 81, 61, 65, 80, 94,
                                          66, 74, 94, 54, 97, 63,
                                          85, 70, 90, 76, 54, 53,
                                          56, 99, 90, 83, 68, 51))
facs          <- c("subject", "sequence", "treatment", "period")
example[facs] <- lapply(example[facs], factor) # factorize the data
txt           <- paste("nested model : period, subject(sequence), treatment",
                       "\nsimple model : period, subject, sequence, treatment",
                       "\nheretic model: period, subject, treatment\n\n")
result        <- data.frame(model = c("nested", "simple", "heretic"),
                            PE = NA, lower = NA, upper = NA, BE = "fail", na = 0)
for (i in 1:3) {
  if (result$model[i] == "nested") { # bogus nested model (guidelines)
    nested         <- lm(log(Y) ~ period +
                                  subject %in% sequence +
                                  treatment, data = example)
    result$PE[i]   <- 100 * exp(coef(nested)[["treatmentT"]])
    result[i, 3:4] <- 100 * exp(confint(nested, level = 0.90)["treatmentT", ])
    result[i, 6]   <- sum(is.na(coef(nested)))
  }
  if (result$model[i] == "simple") { # simple model (subjects are uniquely coded)
    simple         <- lm(log(Y) ~ period +
                                  subject +
                                  sequence +
                                  treatment, data = example)
    result$PE[i]   <- 100 * exp(coef(simple)[["treatmentT"]])
    result[i, 3:4] <- 100 * exp(confint(simple, level = 0.90)["treatmentT", ])
    result[i, 6]   <- sum(is.na(coef(simple)))
  }
  if (result$model[i] == "heretic") { # heretic model (without sequence)
    heretic        <- lm(log(Y) ~ period +
                                  subject +
                                  treatment, data = example)
    result$PE[i]   <- 100 * exp(coef(heretic)[["treatmentT"]])
    result[i, 3:4] <- 100 * exp(confint(heretic, level = 0.90)["treatmentT", ])
    result[i, 6]   <- sum(is.na(coef(heretic)))
  }
  # rounding acc. to guidelines
  if (round(result[i, 3], 2) >= 80 & round(result[i, 4], 2) <= 125)
    result$BE[i] <- "pass"
}
# cosmetics
result$PE     <- sprintf("%6.2f%%", result$PE)
result$lower  <- sprintf("%6.2f%%", result$lower)
result$upper  <- sprintf("%6.2f%%", result$upper)
names(result)[c(3:4, 6)] <- c("lower CL", "upper CL", "NE")
cat(txt); print(result, row.names = FALSE)

# nested model : period, subject(sequence), treatment 
# simple model : period, subject, sequence, treatment 
# heretic model: period, subject, treatment
# 
#    model      PE lower CL upper CL   BE NE
#   nested 102.82%   87.25%  121.17% pass 13
#   simple 102.82%   87.25%  121.17% pass  1
#  heretic 102.82%   87.25%  121.17% pass  0

As already outlined above, the nested model recommended in all [sic] guidelines is over-specified because subjects are uniquely coded. In the example we get 13 not estimable (aliased) effects (in the output of R lines with NA, in SAS ., and in Win­Non­lin not estimable). Correct, because we asking for something the data cannot provide.⁸⁰ In the simple mod­el only one effect cannot be estimated. Even sequence can be removed from the model. I call it he­re­tic because regulators will grill you if you are using it. It was the model proposed by Westlake^{20 43} and I used it in hundreds (‼) of my stud­ies. Note that the results of all models are exactly the same; if you don’t believe me, try it with one of your stud­ies.

A ‘Positive List’ was published by the German regulatory authority, i.e., for 90 drugs BE was not required.⁸³ In order to comply with the European Note for Guidance of 2001⁸⁴ it had to be removed by the BfArM.

Two (of five) sessions of the BioInternational ’92 conference in Bad Homburg dealt with BE of Highly Variable Drugs.^{85 86} Vari­ous approaches have been discussed: Multiple dose instead of single dose studies, metabolite instead of the parent compound, stable isotope tech­niques,⁸⁷ add-on designs, and – for the first time – replicate designs.

Although the BioInternational 2 in Munich 1994 was with over 600 participants the largest in the series, no sub­stan­ti­al progress for HVD(P)s was achieved.⁸⁸ Following a suggestion⁸⁹ at a joint AAPS/FDA workshop in 1995 widening the conventional acceptance limits of 80.00 – 125.00% was considered.⁹⁰

“For some highly variable drugs and drug products, the bioequivalence standard should be modified by changing the BE limits while maintaining the current confidence interval at 90%. […] the bioequivalence limits should be determined based in part upon the intrasubject varia­bility for the reference product.

A hot topic ever since… Why are we discussing it for 35 (‼) years (since the first Bio­Inter­national conference)? Is it really that com­pli­cated⁹¹ or are we too stupid?

Studies in steady-state were proposed as an option for HVD(P)s in a European draft guideline⁹² but was removed from the final version of 2001.⁸⁴

Validation of bioanalytical methods^{93 94 95 96} was partly covered in Australia and Canada. However, no specific guideline existed. A series of conferences (informally known as ‘Crys­tal City’) was initiated in 1990.⁹⁷ Procedures stated in the conference report⁹⁸ were discussed at the Bio­In­ter­na­tio­nal 2 in Munich 1994 and quickly adopted by bioanalytical sites. Updates were subsequently published.^{99 100}

TODO: SUPAC (FDA)

top of section ↩︎ previous section ↩︎

21^st century

After a wealth of – controversal – publications in the 1990s,^{72 73 74 101 102 103 104 105 106 107 108 109} the FDA introduced two new concepts as alternatives to average bio­equi­va­lence (ABE), namely population bioequivalence (PBE) and individual bio­equi­va­lence (IBE).¹¹⁰ ABE focuses only on the comparison of po­pu­lation averages of the PK metrics and not their variances of formulations. It does also not assess a sub­ject-by-for­mu­lat­ion interaction variance, that is, the variation in the average \(\small{\text{T}}\) and \(\small{\text{R}}\) difference among individuals. In contrast, PBE and IBE include com­pa­ri­sons of both averages and variances of PK metrics. The PBE approach assesses total variability of the PK metrics in the population. The IBE approach assesses within-subject variability for the \(\small{\text{T}}\) and \(\small{\text{R}}\) formulations, as well as the sub­ject-by-formulation interaction.
Demonstrated PBE would support ‘Prescribability’ (i.e., a drug naïve patient could start treatment with a generic), whereas IBE support ‘Switchability’ (i.e., a patient could switch formulations during treatment).¹⁰⁹ Contrary to ABE, both PBE and IBE require studies in a full replicate design, which means that both \(\small{\text{T}}\) and \(\small{\text{R}}\) are administered twice. The acceptance limits for ABE were kept at 80.00–125.00% but for the others scaling to the variability of the reference was possible. That would mean an incentive for test formulations with lower variability than the reference but a penalty for ones with a higher variability.

However, the underlying statistical concepts were not trivial and the result practically incomprehensible for non-statisticians. Furthermore, both approaches had a discontinuity (when moving from constant- to reference-scaling), which lead to an inflated type I error (patient’s risk) of approximately 6.5% if CV_wR 18.1–20.2%.^{110 111 112}
The PBE/IBE faced criticism, e.g., »responses [to the guidance] were still doubt-filled as to whether the new bioequivalence criteria really provided added value compared to average bioequivalence«¹¹³ and was regarded a »‘theoretical’ solution to a ‘thoretical’ problem«¹¹⁴ leading to its omission from a subsequent guid­ance,¹¹⁵ and a return to conventional ABE.¹¹⁶

“[ABE should suffice based upon grounds of] ‘practicality, plausibility, historical adequacy, and purpose’ and ‘because we have better things to do.’ […] ‘Statisticians have a bad track record in bioequivalence, […] the literature is full of ludicrous recommendations from statisticians, […] regulatory recommendations (of dubious validity) have been hastily implemented, and practical realities have been ignored’.

I remember a Dutch regulator standing up in the BioInternational conference (London 2003) saying: »I’m glad that PBE and IBE are dead. I never understood them.«

Poland happily adopted Germany’s ‘Positive List’⁸³ only when it wanted to join the European Union to learn that in the mean­time Germany abandoned it. Until 2015 a similar (but shorter) list existed in The Netherlands for national market authori­sa­tions only. Must have been a schizophrenic situation for assessors of the MEB: In the morning a dossier for national MA with­out any in vivo comparison → ☑. In the afternoon another dossier of the same product in the course of a European submission. BE performed, but lower 90% CI 79.99% → ☒. Bizarre.
Until 2012 Denmark required for NTIDs that the 90% CI had to include 100% (i.e., that there is no significant treatment effect). Bizarre as well. For details see Example 3 in this article.

The first bioanalytical method validation guidance was published by the FDA in 2001 and revised in 2018.^{118 119} Before the European draft guideline was published in 2009,¹²⁰ some inspectors raised an eyebrow if sites worked according to the FDA’s guidance.

“The validation of bioanalytical methods and the analysis of study samples should be per­form­ed in accordance with the principles of Good Laboratory Practice (GLP). However, as human bio­ana­ly­ti­cal studies fall outside of the scope of GLP, as defined in Directive 2004/10/EC, the sites con­duct­ing the human studies are not required to be monitored as part of a national GLP compliance programme.

Well roared, lions! My CRO was GLP-certified since 1991, although we performed only phase I studies. In other countries (e.g., Spain), this was not possible. In Germany GLP is subject to state law. Hence, it was possible to get certified in one federal state but not in another… However, this ‘issue’ was resolved with the final guideline published in 2011¹²¹ and the ICH M10 guideline of 2022,¹²² superseding all local guidelines.

TODO: BCS-based biowaivers, reference-scaling, two-stage designs, NTIDs, current guidelines in various jurisdictions…

Still unresolved or not harmonized issues:

1. Scaled ABE for HVD(P)s (RSABE¹²³ or ABEL^{124 125 126});
   control of the type I error,¹²⁷ agreement on which of the metrics can be scaled, outliers^{124 125}
2. Method for NTIDs (fixed narrower acceptance limits¹²⁴ or reference-scaling^{123 128})
3. Comparison of ‘early exposure’¹²⁹ if clinically relevant? (\(\small{t_\text{max}}\) by a nonparametric method or first partial \(\small{AUC}\)); see also this article
4. Cut-off times of partial \(\small{AUC}\textsf{s}\) (based on PD – like the FDA or PK – like the EMA?)
5. Alternative surrogate for the rate of absorption (\(\small{C_\text{max}/AUC}\)^{31 32 33})?
6. Reduce variability of \(\small{AUC}\)¹³⁰ of HVDs by using \(\small{AUC/\hat{\lambda}_z}\)?
7. Studies in fed state mandatory?
8. Multiple dose studies of modified release products really¹³¹ necessary?
9. Adaptive sequential two-stage designs (only exact or simulation-based as well?)
10. Potency-correction if measured contents differ by more than 5% (arbitrary)

A word of warning: The textbooks dealing with statistics (marked with ★ in the references) are rather tough cookies and not recommended for beginners.

top of section ↩︎ previous section ↩︎