Skip to contents

Flurbiprofen (Zhang 2018)

Source: vignettes/articles/Zhang_2018_flurbiprofen.Rmd
Zhang_2018_flurbiprofen.Rmd

Model and source

  • Citation: Zhang J, Zhang H, Zhao L, Gu J, Feng Y, An H. Population pharmacokinetic modeling of flurbiprofen, the active metabolite of flurbiprofen axetil, in Chinese patients with postoperative pain. J Pain Res. 2018;11:3061-3070. doi:10.2147/JPR.S176475. EXPRESSION OF CONCERN: Population Pharmacokinetic Modeling Of Flurbiprofen, The Active Metabolite Of Flurbiprofen Axetil, In Chinese Patients With Postoperative Pain [Expression Of Concern]. J Pain Res. 2019;12:2799. doi:10.2147/JPR.S231531. The Editor-in-Chief and Publisher raise concerns about the validity of the data: a reader alleged that authors Zhang H, Feng Y and Gu J had previously published the same results in Chin J Anaesthesiol. 2011;31(04):432-434 (doi:10.3760/cma.j.issn.0254-1416.2011.04.012), and that the clinical trial registration number ChiCTR-TRC-11001791 cited in the article was taken from an unrelated study. Journal review confirmed the registered trial differs from the published study in age range (60-120 vs 18-72 years), indication (cognitive dysfunction in the elderly vs postoperative pain) and size (100 vs 72 participants), and that the number had been wrongly ascribed in the Methods as an ethics approval number. Two editorial board members compared the two articles but could not conclude with certainty that the data were from the same study. The authors and their institution did not respond to a request for an explanation and the raw data. The notice is not a retraction; the Editor-in-Chief alerts readers to reach their own conclusions. CORRIGENDUM: Population Pharmacokinetic Modeling of Flurbiprofen, the Active Metabolite of Flurbiprofen Axetil, in Chinese Patients with Postoperative Pain [Corrigendum]. J Pain Res. 2020;13:553. doi:10.2147/JPR.S251590. The journal’s investigation determined that the study was not a clinical trial and that the CTR number ChiCTR-TRC-11001791 was used in error; the Methods sentence on ethics approval is corrected to read Ethics Approval No. 2010-30 (Medical Ethics Committee of Peking University People’s Hospital). All relevant study data and ethical approval documents were reviewed and verified by the authors’ institution, the data-integrity and prior-publication issues were resolved, and the journal notes that all concerns raised in the Expression of Concern have been resolved.
  • Description: One-compartment IV population PK plus Holford-Sheiner effect-compartment for cerebrospinal fluid (CSF) disposition of flurbiprofen, the active metabolite of flurbiprofen axetil, in Chinese adults with postoperative pain receiving 1 mg/kg IV flurbiprofen axetil (Zhang 2018, Tables 1-2, Eq. 3 covariate form). Final-model typical values CL = 1.55 L/h, Vd = 7.91 L, plasma-CSF equilibration rate Ke = 0.0015/h; linear-multiplicative covariate effects of weight and height on Ke centered on the population medians (68.5 kg, 165 cm). NOTE: the source article was the subject of a Journal of Pain Research Expression of Concern (2019;12:2799) over data integrity, overlap with an earlier Chinese-language report, and a clinical trial registration number wrongly given as the ethics approval number. The journal’s 2020 Corrigendum (J Pain Res. 2020;13:553) closed that investigation: the study was not a clinical trial, the correct ethics approval is Peking University People’s Hospital No. 2010-30, the study data and ethics documents were reviewed and verified by the authors’ institution, and all concerns raised in the Expression of Concern were stated to be resolved. No parameter, equation, table or unit was changed; see the reference field and the vignette.
  • Article: https://doi.org/10.2147/JPR.S176475
  • Expression of Concern: https://doi.org/10.2147/JPR.S231531
  • Corrigendum (resolves the Expression of Concern): https://doi.org/10.2147/JPR.S251590

Expression of Concern (2019), resolved by Corrigendum (2020). The source article was the subject of a Journal of Pain Research Expression of Concern (2019;12:2799). A reader alleged that three of the authors (Zhang H, Feng Y, Gu J) had already published the same results in Chin J Anaesthesiol. 2011;31(04):432-434, and that the clinical trial registration number cited in the article (ChiCTR-TRC-11001791) was taken from an unrelated study (cognitive dysfunction in the elderly, ages 60-120, n = 100) and wrongly given in the Methods as an ethics approval number. At that point the authors had not supplied the raw data.

The journal’s Corrigendum (J Pain Res. 2020;13:553) closed the investigation. It found that the study was not a clinical trial and that the CTR number was used in error. It corrects the Methods to give the ethics approval as Peking University People’s Hospital Medical Ethics Committee, Ethics Approval No. 2010-30. It also states that all relevant study data and ethical approval documents were reviewed and verified by the authors’ institution, that the data-integrity and prior-publication issues were resolved, and that all concerns raised in the Expression of Concern have been resolved. The Corrigendum changes no parameter, equation, table or unit, so the model below is unaffected. Both notices are recorded in the model file’s reference metadata so they travel with the model. See “Assumptions and deviations” below.

Population

The model was developed at Peking University People’s Hospital from 72 Chinese adults undergoing surgery under subarachnoid anesthesia who received a single 1 mg/kg intravenous injection of flurbiprofen axetil (Zhang 2018 Table 1). Mean age was 50.5 years (SD 13.6, range 18-72), mean body weight 67.5 kg (SD 11.5, range 45-96, median 68.5), and mean height 164.5 cm (SD 9.1, range 149-185, median 165). The cohort was 45 female and 27 male (62.5% female). Each subject contributed one paired plasma + cerebrospinal-fluid sample drawn simultaneously at one of nine post-dose times (5, 10, 15, 20, 25, 30, 35, 40, or 45 min) under a stratified random-allocation design with 8 subjects per time-point group. Observed plasma flurbiprofen concentrations were 3.48-14.56 mg/L (mean 8.30, SD 2.85) and observed CSF concentrations were 0-20.80 ng/mL (mean 5.58, SD 3.66) (Zhang 2018 Table 1).

The same information is available programmatically via readModelDb("Zhang_2018_flurbiprofen")$population.

Source trace

Per-parameter origin is recorded as an in-file comment next to each ini() entry in inst/modeldb/specificDrugs/Zhang_2018_flurbiprofen.R. The table below collects them for review.

Equation / parameter Value Source location
Structural PK model 1-compartment IV with first-order elimination Zhang 2018 Results “Population pharmacokinetic modeling” paragraph; Methods Eq. 1
CSF disposition Holford-Sheiner effect compartment linked to plasma Zhang 2018 Methods, end of “Population pharmacokinetic model” subsection; Results final-model equations p. 3064
Covariate model form Linear-multiplicative [1 + theta2 * (cov - median)] on Ke Zhang 2018 Methods Eq. 3 + Results final-model Ke equation, p. 3064
lcl (CL) log(1.55) L/h Table 2 final model: theta_CL = 1.55 L/h
lvc (Vd) log(7.91) L Table 2 final model: theta_Vd = 7.91 L
lkeo (Ke / ke0) log(0.0015) 1/h Table 2 final model: theta_Ke = 0.0015 /h
e_wt_keo -0.0080 (1/kg) Table 2 final model: theta_weight = -0.0080
e_ht_keo -0.0162 (1/cm) Table 2 final model: theta_height = -0.0162
etalcl IIV variance 0.5217 omega_CL = 82.75% (Table 2) -> omega^2 = log(1 + 0.8275^2)
etalvc IIV variance 0.00789 omega_Vd = 8.90% (Table 2) -> omega^2 = log(1 + 0.089^2)
etalkeo IIV variance 0.000484 omega_Ke = 2.20% (Table 2) -> omega^2 = log(1 + 0.022^2)
addSd (plasma residual) 0.0023 mg/L Table 2 final model: sigma = 0.0023 mg/L
addSd_Ccsf (CSF residual) 0.0023 mg/L Same source sigma; replicated for the per-endpoint nlmixr2 convention
Reference median weight 68.5 kg Table 1 median body weight
Reference median height 165 cm Table 1 median height
d/dt(central) = -kel * central n/a One-compartment IV with kel = CL/Vd
d/dt(effect) = keo * (Cc - effect) n/a Holford-Sheiner effect-compartment link, paper Methods

The paper also reports a derived population time-to-maximum CSF concentration of 10 minutes (Zhang 2018 Discussion p. 3066); this is reproduced by the simulation below rather than encoded as an independent parameter.

Virtual cohort

Individual-level data from the 72-patient cohort are not publicly available. The virtual cohort below approximates the demographic spread reported in Zhang 2018 Table 1: body weight and height are sampled from independent normal distributions truncated to the paper’s reported ranges, and sex is sampled to match the 62.5% female fraction (sex is not a model covariate and is carried as a descriptive label). Each subject receives a single 1 mg/kg IV bolus of flurbiprofen axetil at time 0; observations are collected densely from 0 to 1 hour so the simulated plasma + CSF profiles overlay the paper’s 5-45 min observation window.

set.seed(20180925)
n_subj <- 200L

cohort <- tibble(
  id        = seq_len(n_subj),
  WT        = pmin(pmax(rnorm(n_subj, mean = 67.5, sd = 11.5), 45), 96),
  HT        = pmin(pmax(rnorm(n_subj, mean = 164.5, sd = 9.1), 149), 185),
  SEXF      = rbinom(n_subj, 1, prob = 0.625),
  treatment = factor("1 mg/kg IV flurbiprofen axetil")
)

Simulation

Each subject receives a single IV bolus of 1 mg/kg flurbiprofen axetil into the central compartment at time 0 (the paper treats flurbiprofen axetil as the dosed species and reports plasma flurbiprofen concentrations as the analytical target; see Assumptions and deviations for the molecular-weight caveat). Observations span 0 to 1 h with 1-min resolution so both the 5-45 min observation window of the paper and the broader plasma decay profile are visible.

sim_hours <- 1.0
obs_grid  <- sort(unique(c(seq(0, sim_hours, by = 1/60),
                           c(5, 10, 15, 20, 25, 30, 35, 40, 45) / 60)))

dose_rows <- cohort |>
  dplyr::mutate(time = 0,
                amt  = 1 * WT,        # 1 mg/kg * body weight (kg) -> mg
                cmt  = "central",
                evid = 1L)

obs_cc <- cohort |>
  tidyr::crossing(time = obs_grid) |>
  dplyr::mutate(amt = NA_real_, cmt = "Cc", evid = 0L)

obs_csf <- cohort |>
  tidyr::crossing(time = obs_grid) |>
  dplyr::mutate(amt = NA_real_, cmt = "Ccsf", evid = 0L)

events <- dplyr::bind_rows(dose_rows, obs_cc, obs_csf) |>
  dplyr::select(id, time, amt, cmt, evid, WT, HT, SEXF, treatment) |>
  dplyr::arrange(id, time, dplyr::desc(evid))

stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid", "cmt")])))
mod <- rxode2::rxode2(readModelDb("Zhang_2018_flurbiprofen"))
#> ℹ parameter labels from comments will be replaced by 'label()'
conc_unit <- mod$units[["concentration"]]
sim <- rxode2::rxSolve(
  mod, events = events,
  keep = c("WT", "HT", "SEXF", "treatment"),
  returnType = "data.frame"
)

sim_cc   <- sim |> dplyr::filter(CMT == 3)   # Cc   = plasma observation endpoint (DVID 1)
sim_ccsf <- sim |> dplyr::filter(CMT == 4)   # Ccsf = CSF effect-compartment endpoint (DVID 2)

Replicate published figures

Figure 4: plasma + CSF VPC over the 5-45 min observation window

Zhang 2018 Figure 4 shows the prediction-corrected VPC for plasma (Panel A) and CSF (Panel B) flurbiprofen across the 5-45 min observation window. The plasma envelope spans roughly 4-13 mg/L, and the CSF envelope rises from near 0 at 5 min to roughly 0-20 ng/mL by 45 min. The block below produces the analogous 5-50-95 percentile envelopes from the 200-subject virtual cohort.

sim_cc |>
  dplyr::filter(time >= 5/60, time <= 45/60) |>
  dplyr::mutate(time_min = time * 60) |>
  dplyr::group_by(time_min, treatment) |>
  dplyr::summarise(
    Q05 = quantile(Cc, 0.05, na.rm = TRUE),
    Q50 = quantile(Cc, 0.50, na.rm = TRUE),
    Q95 = quantile(Cc, 0.95, na.rm = TRUE),
    .groups = "drop"
  ) |>
  ggplot(aes(time_min, Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.25) +
  geom_line(linewidth = 0.8) +
  labs(x = "Time after dose (min)",
       y = paste0("Plasma flurbiprofen (", conc_unit, ")"),
       title = "Replicates Figure 4 Panel A: plasma VPC",
       caption = "5-50-95 percentile envelope across 200 simulated subjects.")

sim_ccsf |>
  dplyr::filter(time >= 5/60, time <= 45/60) |>
  dplyr::mutate(time_min = time * 60,
                Ccsf_ngml = Ccsf * 1000) |>   # mg/L -> ng/mL
  dplyr::group_by(time_min, treatment) |>
  dplyr::summarise(
    Q05 = quantile(Ccsf_ngml, 0.05, na.rm = TRUE),
    Q50 = quantile(Ccsf_ngml, 0.50, na.rm = TRUE),
    Q95 = quantile(Ccsf_ngml, 0.95, na.rm = TRUE),
    .groups = "drop"
  ) |>
  ggplot(aes(time_min, Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.25) +
  geom_line(linewidth = 0.8) +
  labs(x = "Time after dose (min)",
       y = "CSF flurbiprofen (ng/mL)",
       title = "Replicates Figure 4 Panel B: CSF VPC",
       caption = "5-50-95 percentile envelope across 200 simulated subjects.")

Typical-value plasma + CSF profile

Reproducing the typical-value (zeroed random effects) profile gives a clean view of the structural model behaviour. The plasma curve is the expected exponential decay with kel = CL/Vd = 1.55/7.91 = 0.196 /h (plasma t1/2 = 3.54 h); the CSF curve rises near-linearly because the slow ke0 = 0.0015 /h keeps the effect compartment far from its eventual quasi-steady-state.

mod_typical <- mod |> rxode2::zeroRe()

typical_cohort <- tibble(
  id        = 1L,
  WT        = 68.5,
  HT        = 165,
  SEXF      = 1L,
  treatment = factor("Typical subject (68.5 kg, 165 cm), 1 mg/kg IV")
)
typical_dose <- typical_cohort |>
  dplyr::mutate(time = 0, amt = 1 * WT, cmt = "central", evid = 1L)
typical_obs_cc <- typical_cohort |>
  tidyr::crossing(time = obs_grid) |>
  dplyr::mutate(amt = NA_real_, cmt = "Cc", evid = 0L)
typical_obs_csf <- typical_cohort |>
  tidyr::crossing(time = obs_grid) |>
  dplyr::mutate(amt = NA_real_, cmt = "Ccsf", evid = 0L)
typical_events <- dplyr::bind_rows(typical_dose, typical_obs_cc, typical_obs_csf) |>
  dplyr::select(id, time, amt, cmt, evid, WT, HT, SEXF, treatment) |>
  dplyr::arrange(id, time, dplyr::desc(evid))

sim_typical <- rxode2::rxSolve(
  mod_typical, events = typical_events,
  keep = c("WT", "HT", "SEXF", "treatment"),
  returnType = "data.frame"
)
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalkeo'
sim_typ_cc  <- sim_typical |> dplyr::filter(CMT == 3)
sim_typ_csf <- sim_typical |> dplyr::filter(CMT == 4)

bind_rows(
  sim_typ_cc  |> dplyr::transmute(time_min = time * 60,
                                  value = Cc,
                                  panel = "Plasma (mg/L)"),
  sim_typ_csf |> dplyr::transmute(time_min = time * 60,
                                  value = Ccsf * 1000,
                                  panel = "CSF (ng/mL)")
) |>
  ggplot(aes(time_min, value)) +
  geom_line(linewidth = 1) +
  facet_wrap(~panel, scales = "free_y") +
  labs(x = "Time after dose (min)",
       y = NULL,
       title = "Typical-value plasma vs CSF profile",
       caption = "68.5 kg, 165 cm subject (population medians); 1 mg/kg IV.")

peak_plasma <- sim_typ_cc  |>
  dplyr::summarise(Cmax = max(Cc, na.rm = TRUE),
                   Tmax = time[which.max(Cc)] * 60)
peak_csf <- sim_typ_csf |>
  dplyr::summarise(Cmax = max(Ccsf, na.rm = TRUE) * 1000,
                   Tmax = time[which.max(Ccsf)] * 60)

knitr::kable(
  tibble::tibble(
    Quantity = c("Plasma Cmax (mg/L)",
                 "Plasma Tmax (min)",
                 "CSF Cmax (ng/mL)",
                 "CSF Tmax (min)"),
    `Zhang 2018 reference` = c(
      "~8.5 (dose / Vd = 68.5/7.91)",
      "0 (IV bolus)",
      "0-20.8 observed across 5-45 min",
      "~10 (paper Discussion p. 3066: 'time to maximum concentration of 10 minutes')"),
    Simulation = c(sprintf("%.2f", peak_plasma$Cmax),
                   sprintf("%.2f", peak_plasma$Tmax),
                   sprintf("%.2f", peak_csf$Cmax),
                   sprintf("%.2f", peak_csf$Tmax))
  ),
  caption = "Typical-value peak metrics: simulation vs Zhang 2018 reported values."
)
Typical-value peak metrics: simulation vs Zhang 2018 reported values.
Quantity Zhang 2018 reference Simulation
Plasma Cmax (mg/L) ~8.5 (dose / Vd = 68.5/7.91) 8.66
Plasma Tmax (min) 0 (IV bolus) 0.00
CSF Cmax (ng/mL) 0-20.8 observed across 5-45 min 11.79
CSF Tmax (min) ~10 (paper Discussion p. 3066: ‘time to maximum concentration of 10 minutes’) 60.00

The simulated CSF Tmax over the 0-60 min window is at the end of the window because, with ke0 = 0.0015 /h, the CSF is still in the rising linear-buildup phase across the entire observation window. The paper’s Discussion statement of ‘time to maximum concentration of 10 minutes’ appears to refer to plasma Tmax (a 1 mg/kg IV bolus reaches Cmax at t = 0 in a one-compartment model; the ‘peak observed’ would be the first sampling time, which is 5-10 min) rather than CSF Tmax. The CSF Cmax is not reached within the paper’s observation window; the simulation shows this is consistent with the structural model.

PKNCA validation

PKNCA is run on the simulated plasma profile across the 0-1 h window (a single-dose IV analysis). Reported NCA values are compared against the literature values cited by Zhang 2018 (Suri 2002, Qayyum 2008, Galasko 2003).

sim_nca <- sim_cc |>
  dplyr::filter(time <= sim_hours) |>
  dplyr::transmute(id, time, Cc, treatment)

dose_df <- cohort |>
  dplyr::transmute(id, time = 0, amt = 1 * WT, treatment)

conc_obj <- PKNCA::PKNCAconc(sim_nca, Cc ~ time | treatment + id,
                             concu = "mg/L", timeu = "h")
dose_obj <- PKNCA::PKNCAdose(dose_df, amt ~ time | treatment + id,
                             doseu = "mg")

intervals <- data.frame(
  start      = 0,
  end        = 1,
  cmax       = TRUE,
  tmax       = TRUE,
  auclast    = TRUE,
  aucinf.obs = TRUE,
  half.life  = TRUE
)

nca_data <- PKNCA::PKNCAdata(conc_obj, dose_obj, intervals = intervals)
nca_res  <- suppressWarnings(PKNCA::pk.nca(nca_data))
knitr::kable(summary(nca_res),
             caption = "Simulated single-dose plasma NCA over 0-1 h window (n = 200).")
Simulated single-dose plasma NCA over 0-1 h window (n = 200).
Interval Start Interval End treatment N AUClast (h*mg/L) Cmax (mg/L) Tmax (h) Half-life (h) AUCinf,obs (h*mg/L)
0 1 1 mg/kg IV flurbiprofen axetil 200 7.40 [21.2] 8.38 [19.2] 0.000 [0.000, 0.000] 4.61 [3.77] 42.7 [87.2]

Comparison against literature CL and Vd

Zhang 2018 Discussion (p. 3066) compares the final-model CL = 1.55 L/h and Vd = 7.91 L against three published flurbiprofen studies in healthy volunteers: Suri et al. 2002 (CL = 1.52 L/h, Vd = 6.00 L; US subjects, oral 100 mg S-flurbiprofen), Qayyum et al. 2008 (CL = 1.89 L/h, Vd = 12.07 L; Pakistani subjects), and Galasko et al. 2003 (CL = 1.465 L/h after multidose R-flurbiprofen). The simulated CL and Vd derived from PKNCA should track the Zhang 2018 typical values to within sampling noise.

nca_tbl <- as.data.frame(nca_res$result)

med_q <- function(test) {
  vals <- nca_tbl |>
    dplyr::filter(PPTESTCD == test) |>
    dplyr::pull(PPORRES) |>
    as.numeric()
  c(median = median(vals, na.rm = TRUE),
    q05    = quantile(vals, 0.05, na.rm = TRUE, names = FALSE),
    q95    = quantile(vals, 0.95, na.rm = TRUE, names = FALSE))
}

sim_cmax  <- med_q("cmax")
sim_aucl  <- med_q("auclast")
sim_auci  <- med_q("aucinf.obs")
sim_thalf <- med_q("half.life")

# Derive CL and Vd from AUCinf and Cmax (per-subject; uses each subject's WT)
nca_wide <- nca_tbl |>
  dplyr::filter(PPTESTCD %in% c("cmax", "aucinf.obs", "half.life")) |>
  dplyr::select(treatment, id, PPTESTCD, PPORRES) |>
  tidyr::pivot_wider(names_from = PPTESTCD, values_from = PPORRES) |>
  dplyr::mutate(dplyr::across(c(cmax, aucinf.obs, half.life), as.numeric))

nca_wide <- nca_wide |>
  dplyr::inner_join(cohort |> dplyr::select(id, WT), by = "id") |>
  dplyr::mutate(
    CL_Lph = (1 * WT) / aucinf.obs,                  # CL = dose / AUCinf
    Vd_L   = CL_Lph * half.life / log(2)             # Vd = CL * t1/2 / ln2 (one-cmt)
  )

cl_med <- median(nca_wide$CL_Lph, na.rm = TRUE)
vd_med <- median(nca_wide$Vd_L,   na.rm = TRUE)
cl_lo  <- quantile(nca_wide$CL_Lph, 0.05, na.rm = TRUE, names = FALSE)
cl_hi  <- quantile(nca_wide$CL_Lph, 0.95, na.rm = TRUE, names = FALSE)
vd_lo  <- quantile(nca_wide$Vd_L, 0.05, na.rm = TRUE, names = FALSE)
vd_hi  <- quantile(nca_wide$Vd_L, 0.95, na.rm = TRUE, names = FALSE)

knitr::kable(
  tibble::tibble(
    Quantity = c("Plasma Cmax (mg/L)",
                 "AUCinf (mg*h/L)",
                 "Plasma t1/2 (h)",
                 "Derived CL (L/h)",
                 "Derived Vd (L)"),
    `Zhang 2018 / literature` = c(
      "~8.5 (Zhang Table 1: 8.30 +/- 2.85 observed)",
      "~34 (= 67.47 mg / 1.55 L/h * 0.78 over 0-1 h is preliminary)",
      "3.54 (= ln2 * 7.91 / 1.55)",
      "1.55 (Zhang); 1.52 Suri; 1.89 Qayyum; 1.47 Galasko",
      "7.91 (Zhang); 6.00 Suri; 12.07 Qayyum"
    ),
    `Simulated median (5-95%)` = c(
      sprintf("%.2f (%.2f-%.2f)", sim_cmax["median"], sim_cmax["q05"], sim_cmax["q95"]),
      sprintf("%.2f (%.2f-%.2f)", sim_auci["median"], sim_auci["q05"], sim_auci["q95"]),
      sprintf("%.2f (%.2f-%.2f)", sim_thalf["median"], sim_thalf["q05"], sim_thalf["q95"]),
      sprintf("%.2f (%.2f-%.2f)", cl_med, cl_lo, cl_hi),
      sprintf("%.2f (%.2f-%.2f)", vd_med, vd_lo, vd_hi)
    )
  ),
  caption = "Simulated NCA-derived CL and Vd compared with Zhang 2018 final-model values and literature."
)
Simulated NCA-derived CL and Vd compared with Zhang 2018 final-model values and literature.
Quantity Zhang 2018 / literature Simulated median (5-95%)
Plasma Cmax (mg/L) ~8.5 (Zhang Table 1: 8.30 +/- 2.85 observed) 8.46 (6.14-11.21)
AUCinf (mg*h/L) ~34 (= 67.47 mg / 1.55 L/h * 0.78 over 0-1 h is preliminary) 43.11 (12.74-170.31)
Plasma t1/2 (h) 3.54 (= ln2 * 7.91 / 1.55) 3.53 (1.06-12.20)
Derived CL (L/h) 1.55 (Zhang); 1.52 Suri; 1.89 Qayyum; 1.47 Galasko 1.53 (0.43-5.10)
Derived Vd (L) 7.91 (Zhang); 6.00 Suri; 12.07 Qayyum 7.98 (6.97-9.07)

The PKNCA window is 0-1 h, which is shorter than ideal for an AUCinf estimate when the plasma t1/2 is 3.54 h; the half-life extrapolation within PKNCA still yields a CL estimate close to the Zhang 2018 typical value because the one-compartment exponential decay is well-defined from the available samples.

Assumptions and deviations

  • Expression of Concern and its resolution. The source article carried a Journal of Pain Research Expression of Concern (2019;12:2799, doi:10.2147/JPR.S231531) about data validity, which the 2020 Corrigendum (2020;13:553, doi:10.2147/JPR.S251590) resolved – see the quoted notices under “Model and source” above. The Corrigendum revises only the ethics approval sentence in the Methods; every value used in this model comes from the original 2018 article and is unchanged. The NCA checks below recover the paper’s own typical-value CL and Vd from the simulated profiles. That is a self-consistency check on the published estimates, not independent evidence about the data they were fitted to. The journal did not state whether the cohort overlaps with the Chin J Anaesthesiol. 2011 report; it said only that the issue was resolved. As a precaution, do not pool the two reports as independent studies. The two other flurbiprofen entries in this library – Kumpulainen_2010_flurbiprofen and the Yao_2025_flurbiprofen_r / Yao_2025_flurbiprofen_s enantiomer pair, the latter also with paired plasma/CSF sampling – are useful comparators.
  • Ethics approval number corrected by the 2020 Corrigendum. Zhang 2018 Methods (p. 3062) gives “ChiCTR-TRC-11001791” as the Peking University People’s Hospital Medical Ethics Committee approval number. The Expression of Concern established that this is a clinical trial registration number belonging to an unrelated study. The Corrigendum states the study was not a clinical trial and corrects the sentence to “Ethics Approval No. 2010-30”. population$disease_state records the corrected identifier.
  • Covariate model form: paper text contradicts the equation, the equation wins. Zhang 2018 Methods (p. 3063) introduces continuous covariates with both ‘a power function after normalization to the median value’ and (in the next paragraph) ‘A linear model was used for the continuous covariate candidates,’ followed by Eq. 3 P_i = theta * [1 + theta2 * (cov - median)] * exp(eta_i). The final- model equations on p. 3064 use the linear-multiplicative form explicitly: K_e = theta_Ke * [1 + theta_weight * (weight - median)] * [1 + theta_height * (height - median)] * exp(eta_Ke). The library follows Eq. 3 and the final-model equation, which is the only form consistent with the reported parameter magnitudes (a true power function with exponents -0.008 and -0.016 would change Ke by < 1% across the cohort’s WT and HT ranges and could not explain the IIV reduction from omega_Ke = 0.0871 in the basic model to 0.0220 in the final model).
  • Ke is the effect-compartment equilibration rate constant, not plasma elimination. Zhang 2018 Methods describes theta_Ke as ‘the elimination rate constant’ (p. 3062), but the final model has three independent typical-value parameters (Vd, CL, Ke) and parameterises the effect-compartment ODE with Ke (p. 3064). The plasma elimination rate constant is kel = CL/Vd = 0.196 /h, giving a plasma t1/2 of 3.54 h. The reported Ke = 0.0015 /h (t1/2 ~462 h) is far too slow to be plasma elimination and is consistent with the slow plasma-CSF equilibration the paper aims to capture (CSF concentrations rise from near-zero at 5 min to roughly 9 ng/mL at 45 min, i.e. CSF stays in the linear-buildup phase across the entire observation window). The library names the parameter lkeo per the Park 2001 ketoprofen precedent (Park_2001_ketoprofen.R).
  • Single residual error split across two endpoints. Zhang 2018 Eq. 2 defines a single additive sigma = 0.0023 mg/L on the ‘serum flurbiprofen concentration,’ but the joint Phoenix NLME fit applies the same sigma to both plasma and CSF observations (Figure 1A and 1B IPRED-vs-DV alignment confirms this). The library replicates the paper’s single sigma into two named parameters addSd (plasma) and addSd_Ccsf (CSF) because nlmixr2 requires distinct residual-error parameters per endpoint; both carry the same initial value 0.0023 to honour the source fit.
  • IIV CV-to-omega2 conversion. Zhang 2018 Table 2 reports omega values under the ‘Estimate’ column with row labels including ‘(%)’; the values are interpreted as CV (fraction) per the Phoenix NLME default and converted to internal log-scale variance via omega^2 = log(1 + CV^2). The basic-model omega_Ke = 0.0871 (CV 8.71%) drops to 0.0220 (CV 2.20%) after adding WT and HT covariates, a 75% reduction consistent with the linear-multiplicative covariate form (a power model with the reported coefficients would produce a negligible reduction).
  • No molecular-weight correction for axetil-to-flurbiprofen. The paper dosed flurbiprofen axetil (MW 330.35) and analysed flurbiprofen (MW 244.26); on a molar basis 1 mg/kg axetil yields ~0.74 mg/kg flurbiprofen. The fit absorbs this conversion into the apparent Vd and CL of flurbiprofen, so the library simulation dosed at 1 mg/kg reproduces the paper’s plasma concentrations directly without applying the 0.74 conversion explicitly. Users wishing to dose flurbiprofen directly (e.g. for a hypothetical IV flurbiprofen formulation) should multiply the desired flurbiprofen mg/kg dose by 1/0.74 = 1.35 before passing to the model.
  • Virtual cohort body weight and height distributions. The cohort of 200 subjects samples WT from a normal distribution truncated to 45-96 kg (mean 67.5, SD 11.5) and HT from a normal distribution truncated to 149-185 cm (mean 164.5, SD 9.1), matching Zhang 2018 Table 1 summary statistics. WT and HT are sampled independently even though physiological correlation exists; the simulation does not rely on this correlation because the covariate effects on Ke enter multiplicatively without an interaction term.
  • Sex is not a model covariate. Zhang 2018 stepwise selection did not retain sex as a covariate on any PK parameter (Discussion p. 3066). The virtual cohort carries SEXF as a label only.
  • CYP2C9 genotype not encoded. Zhang 2018 discusses CYP2C9 genotype as a known driver of inter-individual variability in flurbiprofen clearance but did not include it as a covariate due to small sample size (p. 3068). The library follows the paper and does not encode CYP2C9 status.