Skip to contents

Bepirovirsen (Youssef 2024)

Source: vignettes/articles/Youssef_2024_bepirovirsen.Rmd
Youssef_2024_bepirovirsen.Rmd

Model and source

  • Citation: Youssef AS, Ismail M, Han K, Magee M, Nader A. (2024). Population pharmacokinetics of bepirovirsen in healthy participants and participants with chronic hepatitis B virus infection: results from phase 1, 2a, and 2b studies. Infect Dis Ther 13:1515-1530. doi:10.1007/s40121-024-00980-9
  • Description: Three-compartment population PK model for bepirovirsen (a novel antisense oligonucleotide in development for chronic hepatitis B virus (HBV) infection) in healthy participants and participants with chronic HBV infection (Youssef 2024; N = 479 subjects, 11021 plasma concentrations pooled from three studies – a phase 1 single- and multiple-ascending-dose study in healthy volunteers, a phase 2a study and the phase 2b B-Clear study in chronic HBV infection). Disposition is linear and apparent throughout (CL/F, V/F): first-order subcutaneous absorption from a depot with an absorption lag time, distribution to a shallow and a deep peripheral compartment, and first-order elimination from the central compartment. Typical values for a 73 kg healthy participant are CL/F = 3.10 L/h, V2/F = 11.7 L, ka = 0.237 1/h, ALAG1 = 0.232 h, Q3/F = 0.107 L/h, V3/F = 33.2 L (shallow), Q4/F = 0.0428 L/h and V4/F = 63.7 L (deep). Two covariates were retained. Body weight acts on CL/F and V2/F as a power model referenced to 73 kg with ESTIMATED exponents of 0.494 and 1.01 – the clearance exponent is well below the theoretical 0.75. Chronic HBV infection acts as a fractional shift on the absorption lag time (-46.6%) and on the shallow peripheral volume (-33.7%), so the tabulated lag time and V3/F are the HEALTHY-participant values. Age, albumin, race, nucleos(t)ide-analogue treatment status and baseline HBsAg category were screened and not retained; despite the abstract’s claim to the contrary, Asian versus non-Asian race is NOT in the final model. Inter-individual variability is estimated on all eight disposition parameters and is very large on the deep-compartment terms (146% on Q4/F, 198% on V4/F); residual error is proportional (21.2%).
  • Article: https://doi.org/10.1007/s40121-024-00980-9
  • Supplement: https://static-content.springer.com/esm/art%3A10.1007%2Fs40121-024-00980-9/MediaObjects/40121_2024_980_MOESM1_ESM.pdf

Bepirovirsen is a novel antisense oligonucleotide in development for chronic hepatitis B virus (HBV) infection. Youssef 2024 pooled three GSK studies into a single apparent (CL/F, V/F) three-compartment population PK model with first-order subcutaneous absorption and an absorption lag time.

Population

The analysis pooled 479 participants contributing 11,021 plasma concentrations across three studies (Youssef 2024 Table 1 and Supplementary Table 1): the phase 1 single- and multiple-ascending-dose study 213725 / NCT03020745 in healthy volunteers (N = 21 analysed, 75-450 mg SC), the phase 2a study 205695 / NCT02981602 in chronic HBV infection (N = 23, 150-300 mg SC), and the phase 2b B-Clear study 209668 / NCT04449029 in chronic HBV infection (N = 455, 300 mg SC weekly for up to 24 weeks with a 150 mg or placebo step-down in two arms).

96% of the analysis population had chronic HBV infection and 47.5% were on stable nucleos(t)ide-analogue therapy. The population was 63.7% male, 54.3% Asian, 39.3% Caucasian and 6.2% Black; median age 45 years (range 18-77) and median weight 70.2 kg (range 43.3-140.0). Participants with cirrhosis or liver failure, or with moderate-to-severe renal disease, were excluded.

Of the 12,140 available concentrations, 90.8% were used; 4.8% of those were below the 1 ng/mL LLOQ and were handled with the Beal M3 method.

The same information is available programmatically via the model’s population metadata:

pop <- rxode2::rxode(readModelDb("Youssef_2024_bepirovirsen"))$population
#> ℹ parameter labels from comments will be replaced by 'label()'
str(pop, max.level = 1)
#> List of 13
#>  $ species       : chr "human"
#>  $ n_subjects    : num 479
#>  $ n_studies     : num 3
#>  $ age_range     : chr "18-77 years"
#>  $ age_median    : chr "45 years"
#>  $ weight_range  : chr "43.3-140.0 kg"
#>  $ weight_median : chr "70.2 kg"
#>  $ sex_female_pct: num 36.3
#>  $ race_ethnicity: Named num [1:4] 54.3 39.3 6.2 0.2
#>   ..- attr(*, "names")= chr [1:4] "Asian" "Caucasian" "Black" "NativeIndianOrAlaskaNative"
#>  $ disease_state : chr "chronic hepatitis B virus infection (96%, both on and off stable nucleos(t)ide-analogue therapy) pooled with he"| __truncated__
#>  $ dose_range    : chr "75-450 mg subcutaneous: single ascending dose, six doses over 3 weeks, or 150-300 mg weekly for up to 24 weeks"
#>  $ regions       : chr "China Mainland, East Asia, Japan and Other (the four strata of Youssef 2024 Supplementary Table 7)"
#>  $ notes         : chr "Youssef 2024 Table 1 (demographics by study and overall). Three study cohorts: 213725 / NCT03020745 (phase 1, h"| __truncated__

Source trace

The per-parameter origin is recorded as an in-file comment next to each ini() entry in inst/modeldb/specificDrugs/Youssef_2024_bepirovirsen.R. The table below collects them in one place for review. Every number below is from Youssef 2024 Table 2 (“Parameter estimates of the final population PK model”) or its footnote.

Equation / parameter Value Source location
lcl (CL/F) 3.10 L/h Table 2, “CL/F (l/h)” (RSE 2.89%)
lvc (V2/F) 11.7 L Table 2, “V2/F (L)” (RSE 2.47%)
lka (KA) 0.237 1/h Table 2, “KA (h-1)” (RSE 3.48%)
ltlag (ALAG1) 0.232 h (healthy) Table 2, “ALAG1 (h), healthy” (RSE 4.68%)
lq (Q3/F) 0.107 L/h Table 2, “Q3/F (l/h)” (RSE 2.91%)
lvp (V3/F) 33.2 L (healthy) Table 2, “V3/F (L), healthy” (RSE 1.98%)
lq2 (Q4/F) 0.0428 L/h Table 2, “Q4/F (l/h)” (RSE 2.69%)
lvp2 (V4/F) 63.7 L Table 2, “V4/F (L)” (RSE 2.25%)
e_wt_cl 0.494 Table 2, “WT on CL/F”; footnote CL/F = 3.10 x (WT/73)^0.494
e_wt_vc 1.01 Table 2, “WT on V2/F”; footnote V2/F = 11.7 x (WT/73)^1.01
e_dis_chb_tlag -0.466 Table 2, “Chronic HBV infection on ALAG1”; footnote ALAG1 = 0.232 * (1-0.466)
e_dis_chb_vp -0.337 Table 2, “Chronic HBV infection on V3/F”; footnote V3/F = 33.2 * (1-0.337)
IIV variances (8 etas) 0.096 / 0.154 / 0.192 / 0.136 / 0.233 / 0.448 / 1.140 / 1.590 Table 2, “IIV (CV%)” column, first number of each cell
propSd sqrt(0.0449) = 0.21190 Table 2, “RV (CV%)” = 0.0449 (21.2) (RSE 1.59%)
Reference weight 73 kg Table 2 footnote (see Errata)
Structural model three compartments, first-order SC absorption, lag Fig. 1; Methods “Base Structural Model Development”
Residual error form proportional only Methods; Supplementary Table 2 run Poppk-base-v2
Retained covariates WT on CL/F and V2/F; chronic HBV on ALAG1 and V3/F Supplementary Tables 4 and 5

The IIV column is a variance, not an SD

Youssef 2024 Table 2 prints each IIV cell as <estimate> (<CV%>) without saying which scale the estimate is on. The paper’s own two columns settle it. Supplementary Methods Equation 2 defines the reported CV% for a log-normally distributed parameter, so CV = sqrt(exp(omega^2) - 1); back-transforming the eight estimates through that identity must reproduce the printed CV% column.

iiv <- tibble::tribble(
  ~parameter, ~estimate, ~published_cv,
  "KA",       0.192,     45.9,
  "ALAG1",    0.136,     38.2,
  "CL/F",     0.096,     31.7,
  "V2/F",     0.154,     40.8,
  "Q3/F",     0.233,     51.2,
  "V3/F",     0.448,     75.2,
  "Q4/F",     1.140,     146.0,
  "V4/F",     1.590,     198.0
) |>
  mutate(
    cv_if_variance = 100 * sqrt(exp(estimate) - 1),
    cv_if_sd       = 100 * sqrt(exp(estimate^2) - 1)
  )

iiv |>
  dplyr::rename(
    "Parameter"                = parameter,
    "Table 2 estimate"         = estimate,
    "Published CV%"            = published_cv,
    "CV% if estimate is omega^2" = cv_if_variance,
    "CV% if estimate is omega"   = cv_if_sd
  ) |>
  knitr::kable(digits = 1, caption = "Reading the Table 2 IIV estimates as variances reproduces the published CV%; reading them as SDs does not.")
Reading the Table 2 IIV estimates as variances reproduces the published CV%; reading them as SDs does not.
Parameter Table 2 estimate Published CV% CV% if estimate is omega^2 CV% if estimate is omega
KA 0.2 45.9 46.0 19.4
ALAG1 0.1 38.2 38.2 13.7
CL/F 0.1 31.7 31.7 9.6
V2/F 0.2 40.8 40.8 15.5
Q3/F 0.2 51.2 51.2 23.6
V3/F 0.4 75.2 75.2 47.1
Q4/F 1.1 146.0 145.8 163.3
V4/F 1.6 198.0 197.6 339.6 

# The variance reading reproduces every published CV% to within 0.4 percentage
# points; the SD reading is off by tens of points. This is arithmetic on
# transcribed constants -- no simulation, no randomness -- so the bound is tight.
stopifnot(
  max(abs(iiv$cv_if_variance - iiv$published_cv)) < 0.5,
  min(abs(iiv$cv_if_sd       - iiv$published_cv)) > 5
)

# The residual error is on the same footing: Table 2 prints "0.0449 (21.2)".
stopifnot(abs(100 * sqrt(0.0449) - 21.2) < 0.1)

Structural check 1: the model reproduces the printed covariate equations

The Table 2 footnote prints the covariate model in closed form. Solving the packaged model with the random effects zeroed must return individual parameters that match those equations exactly, for any weight and either disease status. This checks the transcription of all ten fixed effects and the covariate algebra in model(). The equations on the right-hand side below are typed from the paper, independently of the model file, so a mis-transcribed value makes this fail.

mod <- readModelDb("Youssef_2024_bepirovirsen")

grid <- tidyr::crossing(WT = c(43.3, 60, 70, 73, 100, 140), DIS_CHB = c(0, 1)) |>
  mutate(id = dplyr::row_number())

check_ev <- grid |>
  tidyr::crossing(time = c(0, 1)) |>
  mutate(amt = NA_real_, evid = 0L, cmt = "central") |>
  arrange(id, time)

chk <- rxode2::rxSolve(
  rxode2::zeroRe(mod), check_ev,
  keep = c("WT", "DIS_CHB"), maxsteps = 200000L
) |>
  as.data.frame() |>
  dplyr::distinct(WT, DIS_CHB, .keep_all = TRUE) |>
  mutate(
    # Youssef 2024 Table 2 footnote, transcribed here independently of the model.
    ref_cl   = 3.10 * (WT / 73)^0.494,
    ref_vc   = 11.7 * (WT / 73)^1.01,
    ref_ka   = 0.237,
    ref_tlag = 0.232 * (1 - 0.466 * DIS_CHB),
    ref_q    = 0.107,
    ref_vp   = 33.2 * (1 - 0.337 * DIS_CHB),
    ref_q2   = 0.0428,
    ref_vp2  = 63.7
  )
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalka', 'etaltlag', 'etalq', 'etalvp', 'etalq2', 'etalvp2'
#> Warning: multi-subject simulation without without 'omega'

rel_err <- c(
  cl = max(abs(chk$cl / chk$ref_cl - 1)),   vc  = max(abs(chk$vc  / chk$ref_vc  - 1)),
  ka = max(abs(chk$ka / chk$ref_ka - 1)),   tlag = max(abs(chk$tlag / chk$ref_tlag - 1)),
  q  = max(abs(chk$q  / chk$ref_q  - 1)),   vp  = max(abs(chk$vp  / chk$ref_vp  - 1)),
  q2 = max(abs(chk$q2 / chk$ref_q2 - 1)),   vp2 = max(abs(chk$vp2 / chk$ref_vp2 - 1))
)
print(signif(rel_err, 3))
#>       cl       vc       ka     tlag        q       vp       q2      vp2 
#> 7.77e-16 5.55e-16 4.11e-15 3.00e-15 8.88e-16 3.00e-15 4.22e-15 2.22e-15

# Deterministic: no IIV, no residual error, no random draw. Any disagreement
# beyond floating-point noise is a transcription or wiring error.
stopifnot(max(rel_err) < 1e-10)

# And the parameters collapse to the bare Table 2 values at the reference
# subject -- 73 kg, healthy.
ref_subject <- chk |> dplyr::filter(WT == 73, DIS_CHB == 0)
stopifnot(
  nrow(ref_subject) == 1L,
  abs(ref_subject$cl   - 3.10)   < 1e-10,
  abs(ref_subject$vc   - 11.7)   < 1e-10,
  abs(ref_subject$tlag - 0.232)  < 1e-10,
  abs(ref_subject$vp   - 33.2)   < 1e-10
)

Structural check 2: the reference weight is 73 kg, not the 70.2 kg median

Supplementary Equation 5 says continuous covariates were centred on their median values, and Table 1 gives an overall median weight of 70.2 kg – so the 73 in the Table 2 footnote looks like a typo. It is not, and the paper’s own simulations prove it.

For a linear model with an apparent clearance, steady-state exposure over a dosing interval is exactly AUCtau = Dose / (CL/F) – no other parameter enters. Table 4 reports simulated steady-state AUCtau for 300 mg weekly at three fixed weights, which turns that identity into a direct test of the reference weight.

ref_wt_test <- tibble::tibble(
  WT           = c(40, 70, 100),
  published    = c(131.0, 99.2, 83.2)   # Youssef 2024 Table 4, AUCtau median (ug*h/mL)
) |>
  mutate(
    auc_ref73 = 300 / (3.10 * (WT / 73)^0.494),
    auc_ref70 = 300 / (3.10 * (WT / 70)^0.494),
    pct_ref73 = 100 * (auc_ref73 / published - 1),
    pct_ref70 = 100 * (auc_ref70 / published - 1)
  )

ref_wt_test |>
  dplyr::rename(
    "Weight (kg)"                = WT,
    "Published AUCtau (ug*h/mL)" = published,
    "Dose/CL, 73 kg reference"   = auc_ref73,
    "% diff (73 kg)"             = pct_ref73,
    "Dose/CL, 70 kg reference"   = auc_ref70,
    "% diff (70 kg)"             = pct_ref70
  ) |>
  knitr::kable(digits = c(0, 1, 1, 2, 1, 2),
               caption = "A 73 kg reference reproduces Table 4; a 70 kg reference leaves a systematic negative bias at every weight.")
A 73 kg reference reproduces Table 4; a 70 kg reference leaves a systematic negative bias at every weight.
Weight (kg) Published AUCtau (ug*h/mL) Dose/CL, 73 kg reference Dose/CL, 70 kg reference % diff (73 kg) % diff (70 kg)
40 131.0 130.3 127.59 -0.6 -2.60
70 99.2 98.8 96.77 -0.4 -2.45
100 83.2 82.8 81.14 -0.4 -2.48 

# Deterministic arithmetic on published constants. The 73 kg reading sits inside
# the Monte Carlo error of a 1000-subject median (about 1.2% RSE); the 70 kg
# reading is biased the same way at all three weights, which sampling error
# cannot produce.
stopifnot(
  max(abs(ref_wt_test$pct_ref73)) < 1.5,
  min(abs(ref_wt_test$pct_ref70)) > 2.0,
  all(ref_wt_test$pct_ref70 < 0)
)

Virtual cohort

Original observed data are not publicly available. The simulations below use virtual cohorts whose covariate distributions approximate the published trial demographics (Youssef 2024 Table 1): a log-normal body weight with median 70.2 kg truncated to the observed 43.3-140.0 kg range, and chronic HBV infection for every subject except in the healthy-volunteer arm.

Each arm holds 200 participants, and rxode2::rxSetSeed() is called once per arm with the same seed so that all arms draw the same individual random effects. That pairing is deliberate: every comparison below is a within-subject contrast (the same virtual patient on a different regimen, at a different weight, or with a different disease status), so common random numbers are both the correct science and a large variance reduction.

set.seed(20240525)

# 1000, not 200. The Ctrough gate below asserts that the under-prediction is
# PRESENT in every row (see the Errata), so the cohort has to estimate it
# precisely enough that the deviation itself is what is being tested rather
# than Monte Carlo noise. rxSetSeed() fixes rxode2's stream per solver
# thread, so the cohort is redrawn when the thread count changes, and at 200
# per arm the worst Ctrough row swung -2.1 / -8.2 / -12.2 percent on
# 2 / 4 / 8 threads -- crossing the -5 bound at 2. At 1000 the same row is
# -8.8 to -11.5, and at 2000 it is -12.9, so the true deviation is around
# -13 percent and the 200-subject estimate was simply noisy.
N_ARM   <- 1000L
WEEK    <- 168             # hours in a weekly dosing interval
WK12    <- 11 * WEEK       # time of the 12th dose
WK24    <- 23 * WEEK       # time of the 24th (final) dose

# Observation grid: weekly troughs across treatment, a fine grid inside the two
# dosing intervals the paper reports NCA for, and weekly samples through the
# 12-week follow-up.
OBS_T <- sort(unique(c(
  seq(0, 24 * WEEK, by = WEEK),
  WK12 + c(seq(0, 24, by = 0.5), seq(30, WEEK, by = 6)),
  WK24 + c(seq(0, 24, by = 0.5), seq(30, WEEK, by = 6)),
  seq(24 * WEEK, 36 * WEEK, by = WEEK)
)))

# Body weights: log-normal, median matched to Table 1, truncated to the observed
# range. sigma = 0.20 on the log scale gives a spread close to the reported
# mean (SD) of 71.2 (14.5) kg.
trial_weights <- function(n) {
  pmin(pmax(stats::rlnorm(n, meanlog = log(70.2), sdlog = 0.20), 43.3), 140.0)
}

# One arm = a dosing schedule plus a covariate set, as a self-contained event
# table. `id_offset` keeps subject IDs disjoint across arms.
make_arm <- function(label, amts, wt, dis_chb = 1, n = N_ARM, id_offset = 0L) {
  subj <- tibble::tibble(
    id = id_offset + seq_len(n), arm = label, WT = wt, DIS_CHB = dis_chb
  )
  doses <- subj |>
    tidyr::crossing(dose_index = seq_along(amts)) |>
    mutate(time = (dose_index - 1) * WEEK, amt = amts[dose_index],
           evid = 1L, cmt = "depot") |>
    select(-dose_index)
  obs <- subj |>
    tidyr::crossing(time = OBS_T) |>
    mutate(amt = NA_real_, evid = 0L, cmt = "central")
  bind_rows(doses, obs) |> arrange(id, time, desc(evid))
}

arm_spec <- tibble::tribble(
  ~label,                                 ~amts,                        ~dis_chb,
  "450 mg weekly x 24 wk",                rep(450, 24),                 1,
  "300 mg weekly x 24 wk",                rep(300, 24),                 1,
  "300 mg x 12 wk, then 150 mg x 12 wk",  c(rep(300, 12), rep(150, 12)), 1,
  "300 mg weekly x 24 wk (healthy)",      rep(300, 24),                 0
)

# Same weights in every arm, so a cross-arm difference is never a weight artefact.
arm_wt <- trial_weights(N_ARM)

events <- do.call(bind_rows, lapply(seq_len(nrow(arm_spec)), function(i) {
  make_arm(arm_spec$label[i], arm_spec$amts[[i]], wt = arm_wt,
           dis_chb = arm_spec$dis_chb[i], id_offset = (i - 1L) * N_ARM)
}))

stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid")])))

Simulation

rxSolve() on an rxUi scales quadratically in the number of subjects passed in one call, so each arm is solved separately. The seed is set inside the loop so every arm draws the same individual random effects.

solve_arm <- function(label) {
  rxode2::rxSetSeed(20240525)
  rxode2::rxSolve(
    mod, events |> dplyr::filter(arm == label),
    keep = c("WT", "DIS_CHB", "arm"),
    # 24 lagged doses over 4032 h push the default step budget.
    maxsteps = 200000L
  ) |>
    as.data.frame()
}

sim <- do.call(bind_rows, lapply(arm_spec$label, solve_arm))
#> ℹ parameter labels from comments will be replaced by 'label()'

# Assert the pairing actually held -- it is load-bearing for the tightened
# cross-arm gates below, so do not assume it.
paired <- sim |>
  group_by(arm, id) |>
  summarise(cl = first(cl), .groups = "drop") |>
  mutate(subject = (id - 1L) %% N_ARM) |>
  group_by(subject) |>
  summarise(spread = diff(range(cl)), .groups = "drop")
stopifnot(max(paired$spread) == 0)

stopifnot(all(sim$Cc[!is.na(sim$Cc)] >= 0))

Structural check 3: mass balance at steady state

For this linear model AUCtau x (CL/F) = Dose holds exactly once steady state is reached, whatever the absorption and distribution parameters do. Checking that identity per subject, against each subject’s own clearance, tests the ODE wiring and the approach to steady state.

Note what this does and does not prove. Because it reads the model’s own cl, a mis-transcribed lcl would move both sides together and leave it green – it is a structural check, not a transcription check. Transcription is covered deterministically by structural check 1 above, and cross-checked at the cohort level at the end of this chunk against the paper’s printed clearance equation.

auc_trapz <- function(time, conc) {
  o <- order(time); time <- time[o]; conc <- conc[o]
  sum(diff(time) * (head(conc, -1) + tail(conc, -1)) / 2)
}

mb <- sim |>
  dplyr::filter(arm == "300 mg weekly x 24 wk", time >= WK24, time <= WK24 + WEEK,
                !is.na(Cc)) |>
  group_by(id, WT) |>
  summarise(auc = auc_trapz(time, Cc), cl_i = first(cl), .groups = "drop") |>
  mutate(
    pct_diff = 100 * (auc * cl_i / 300 - 1),
    # Dose / (CL/F) from Youssef 2024 Table 2's printed equation, transcribed
    # here independently of the model file.
    auc_typical = 300 / (3.10 * (WT / 73)^0.494)
  )

summary(mb$pct_diff)
#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#> -3.00379 -0.12871  0.10327 -0.01212  0.23844  0.46031

# By week 24 the deep compartment (Q4/F = 0.0428 L/h into V4/F = 63.7 L) is very
# close to, but not exactly at, steady state, and how close depends on the
# subject's own draw from a 146% / 198% CV pair -- so assert on the centre and a
# robust quantile, never on the extreme. Realised across three seeds: median
# +0.07% to +0.09%, 90th percentile of |diff| 0.41-0.46%, worst subject 3.4%.
stopifnot(
  abs(median(mb$pct_diff)) < 1,
  stats::quantile(abs(mb$pct_diff), 0.9) < 2
)

# Cohort-level transcription cross-check: the median individual AUCtau should
# track Dose / (CL/F) evaluated at the typical clearance, because the clearance
# eta is log-normal with median 1. This is a cohort statistic on 200 subjects
# with a 31.7% CV on clearance, so it is intrinsically noisy -- realised +2.7%,
# -1.3% and -6.0% across three seeds. 12 keeps headroom while still breaking on
# a mis-transcribed clearance or reference weight.
median_ratio_pct <- 100 * (median(mb$auc) / median(mb$auc_typical) - 1)
cat(sprintf("Cohort-median AUCtau vs the printed equation: %+.2f%%\n",
            median_ratio_pct))
#> Cohort-median AUCtau vs the printed equation: +0.44%
stopifnot(abs(median_ratio_pct) < 12)

Replicate published figures 

sim |>
  dplyr::filter(arm == "300 mg weekly x 24 wk", !is.na(Cc), Cc > 0) |>
  group_by(time) |>
  summarise(
    Q05 = quantile(Cc, 0.05), Q50 = quantile(Cc, 0.50), Q95 = quantile(Cc, 0.95),
    .groups = "drop"
  ) |>
  mutate(week = time / WEEK) |>
  ggplot(aes(week, Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.25, fill = "steelblue") +
  geom_line(colour = "steelblue4") +
  geom_vline(xintercept = 24, linetype = "dashed") +
  scale_y_log10() +
  labs(
    x = "Time (weeks)", y = "Bepirovirsen (ug/mL)",
    title = "300 mg SC weekly for 24 weeks, then follow-up",
    subtitle = "Median with 5th-95th percentile band; dashed line = end of treatment",
    caption = "Replicates Figure 2 of Youssef 2024 (phase 2b B-Clear pcVPC)."
  )
Replicates Figure 2 of Youssef 2024: simulated concentration-time profile for bepirovirsen 300 mg SC weekly for 24 weeks in participants with chronic HBV infection, followed by 12 weeks off treatment.

Replicates Figure 2 of Youssef 2024: simulated concentration-time profile for bepirovirsen 300 mg SC weekly for 24 weeks in participants with chronic HBV infection, followed by 12 weeks off treatment. 

sim |>
  dplyr::filter(time >= WK24, time <= WK24 + WEEK, !is.na(Cc), Cc > 0,
                arm != "300 mg weekly x 24 wk (healthy)") |>
  group_by(arm, time) |>
  summarise(Q50 = quantile(Cc, 0.50), .groups = "drop") |>
  ggplot(aes(time - WK24, Q50, colour = arm)) +
  geom_line() +
  scale_y_log10() +
  scale_colour_brewer(palette = "Dark2") +
  labs(x = "Time since final dose (h)", y = "Bepirovirsen (ug/mL)",
       colour = NULL, title = "Final dosing interval, median profile by regimen") +
  theme(legend.position = "bottom", legend.direction = "vertical")
Simulated profile over the final dosing interval by regimen, showing the dose-proportionality reported in Table 3.

Simulated profile over the final dosing interval by regimen, showing the dose-proportionality reported in Table 3.

Washout after the end of treatment

Youssef 2024 reports that in participants given 300 mg weekly for 24 weeks the median predicted concentration at week 36 was 93% lower than at week 24 (0.0224 vs 0.0017 ug/mL).

washout <- sim |>
  dplyr::filter(arm == "300 mg weekly x 24 wk", time %in% c(24 * WEEK, 36 * WEEK)) |>
  select(id, time, Cc) |>
  tidyr::pivot_wider(names_from = time, values_from = Cc,
                     names_prefix = "t")

reduction_pct <- 100 * (1 - median(washout[[paste0("t", 36 * WEEK)]]) /
                          median(washout[[paste0("t", 24 * WEEK)]]))
cat(sprintf("Median week-24 concentration: %.4f ug/mL\n",
            median(washout[[paste0("t", 24 * WEEK)]])))
#> Median week-24 concentration: 0.0198 ug/mL
cat(sprintf("Median week-36 concentration: %.4f ug/mL\n",
            median(washout[[paste0("t", 36 * WEEK)]])))
#> Median week-36 concentration: 0.0008 ug/mL
cat(sprintf("Reduction: %.1f%% (Youssef 2024 reports 93%%)\n", reduction_pct))
#> Reduction: 95.8% (Youssef 2024 reports 93%)

# The washout depends on the deep compartment, whose IIV is the largest in the
# model, so the realised reduction moves with the cohort draw (94.9% here).
# Bound the claim, not the draw: the paper's qualitative result is "nearly
# completely washed out by 12 weeks", and an upper bound keeps the gate able to
# go red if elimination were ever wired to run away.
stopifnot(reduction_pct > 85, reduction_pct < 99.5)

PKNCA validation

Steady-state NCA over a weekly dosing interval, computed with PKNCA. Two intervals are requested: the interval following the 12th dose (which Table 3 reports for the step-down arm) and the interval following the 24th and final dose.

sim_nca <- sim |>
  dplyr::filter(!is.na(Cc)) |>
  select(id, time, Cc, arm)

conc_obj <- PKNCA::PKNCAconc(as.data.frame(sim_nca), Cc ~ time | arm + id)

dose_df <- events |>
  dplyr::filter(evid == 1) |>
  select(id, time, amt, arm) |>
  as.data.frame()

dose_obj <- PKNCA::PKNCAdose(dose_df, amt ~ time | arm + id)

intervals <- data.frame(
  start   = c(WK12, WK24),
  end     = c(WK12 + WEEK, WK24 + WEEK),
  cmax    = TRUE,
  cmin    = TRUE,
  auclast = TRUE
)

nca_res <- PKNCA::pk.nca(PKNCA::PKNCAdata(conc_obj, dose_obj, intervals = intervals))
nca_tbl <- as.data.frame(nca_res$result) |>
  dplyr::filter(PPTESTCD %in% c("cmax", "cmin", "auclast")) |>
  mutate(visit = ifelse(start == WK12, "Week 12", "Week 24"),
         group = paste0(arm, ", ", visit))

Comparison against published NCA

Youssef 2024 Table 3 reports model-predicted steady-state exposures for the three active regimens simulated in the paper, and Table 4 adds a disease-status contrast. auclast over the dosing interval is the paper’s AUCtau and cmin at the end of the interval is its Ctrough.

published <- tibble::tribble(
  ~group,                                                   ~cmax, ~cmin,  ~auclast,
  "450 mg weekly x 24 wk, Week 24",                          12.6, 0.0336, 147.0,
  "300 mg weekly x 24 wk, Week 24",                           8.4, 0.0224,  98.2,
  "300 mg x 12 wk, then 150 mg x 12 wk, Week 12",             8.4, 0.0205,  97.9,
  "300 mg x 12 wk, then 150 mg x 12 wk, Week 24",             4.2, 0.0120,  49.2,
  "300 mg weekly x 24 wk (healthy), Week 24",                 8.4, 0.0243,  98.3
) |>
  as.data.frame()

cmp <- nlmixr2lib::ncaComparisonTable(
  simulated = nca_tbl |> dplyr::filter(group %in% published$group) |> as.data.frame(),
  reference = published,
  by        = "group",
  units     = c(cmax = "ug/mL", cmin = "ug/mL", auclast = "ug*h/mL"),
  tolerance_pct = 20
)

knitr::kable(
  cmp,
  caption = "Simulated vs. published steady-state exposures (Youssef 2024 Tables 3 and 4). * differs from reference by >20%.",
  align = c("l", "l", "r", "r", "r")
)
Simulated vs. published steady-state exposures (Youssef 2024 Tables 3 and 4). * differs from reference by >20%.
NCA parameter group Reference Simulated % diff
Cmax (ug/mL) 450 mg weekly x 24 wk, Week 24 12.6 12.6 +0.3%
Cmax (ug/mL) 300 mg weekly x 24 wk, Week 24 8.4 8.42 +0.3%
Cmax (ug/mL) 300 mg x 12 wk, then 150 mg x 12 wk, Week 12 8.4 8.42 +0.2%
Cmax (ug/mL) 300 mg x 12 wk, then 150 mg x 12 wk, Week 24 4.2 4.21 +0.3%
Cmax (ug/mL) 300 mg weekly x 24 wk (healthy), Week 24 8.4 8.41 +0.2%
Cmin (ug/mL) 450 mg weekly x 24 wk, Week 24 0.0336 0.0297 -11.6%
Cmin (ug/mL) 300 mg weekly x 24 wk, Week 24 0.0224 0.0198 -11.6%
Cmin (ug/mL) 300 mg x 12 wk, then 150 mg x 12 wk, Week 12 0.0205 0.0184 -10.2%
Cmin (ug/mL) 300 mg x 12 wk, then 150 mg x 12 wk, Week 24 0.012 0.0106 -11.8%
Cmin (ug/mL) 300 mg weekly x 24 wk (healthy), Week 24 0.0243 0.0222 -8.8%
AUClast (ug*h/mL) 450 mg weekly x 24 wk, Week 24 147 148 +0.7%
AUClast (ug*h/mL) 300 mg weekly x 24 wk, Week 24 98.2 98.7 +0.5%
AUClast (ug*h/mL) 300 mg x 12 wk, then 150 mg x 12 wk, Week 12 97.9 98.3 +0.4%
AUClast (ug*h/mL) 300 mg x 12 wk, then 150 mg x 12 wk, Week 24 49.2 49.6 +0.8%
AUClast (ug*h/mL) 300 mg weekly x 24 wk (healthy), Week 24 98.3 98.7 +0.4% 
if (!is.null(attr(cmp, "footnote"))) cat(attr(cmp, "footnote"), "\n")

AUCtau reproduces to within about 3% and Cmax to within about 6% across all five rows, including the step-down arm’s week-12 interval and the healthy-cohort contrast. Ctrough is consistently 14-19% low on this cohort. That last deviation is real and reproducible rather than a cohort artefact; it is characterised in the Errata below and gated separately, rather than being absorbed into a widened tolerance.

# ncaComparisonTable formats "% diff" as text, so parse it back for the gate.
gate <- cmp |>
  mutate(pct = suppressWarnings(as.numeric(gsub("[^0-9.+-]", "", `% diff`))))

auc_cmax <- gate |> dplyr::filter(grepl("AUC|Cmax", .data[["NCA parameter"]]))
ctrough  <- gate |> dplyr::filter(grepl("Cmin|Ctrough", .data[["NCA parameter"]]))

stopifnot(nrow(auc_cmax) == 10L, nrow(ctrough) == 5L)

# AUCtau and Cmax: the worst row realised 5.6%, 9.1% and 5.5% across three
# seeds, so 15 keeps headroom while still breaking on a mis-transcribed
# clearance, dose or volume (any of which shifts these by tens of percent).
stopifnot(max(abs(auc_cmax$pct)) < 15)

# Ctrough: a real, reproducible under-prediction (see Errata). With the 1000-per
# -arm cohort above, all five rows realised -8.8% to -16.9% across 2 and 8
# solver threads and 1000/2000 subjects -- never near zero, and never positive.
# Gated as a two-sided band so the vignette fails both if the deviation vanishes
# silently (which would mean something else changed) and if the terminal phase
# breaks outright, but NOT widened to hide the disagreement: the Cmin rows of
# the table above report it in full. The bound is unchanged from when this gate
# was written; what changed is the cohort size behind the estimate.
stopifnot(max(ctrough$pct) < -5, min(ctrough$pct) > -45)

Simulated exposure by body weight

Table 4 also reports steady-state exposures at three fixed body weights. Because AUCtau = Dose / (CL/F) is exact here, the weight rows are reproduced deterministically from the typical-value model, with no cohort draw involved.

wt_grid <- tibble::tibble(WT = c(40, 70, 100), DIS_CHB = 1, id = 1:3)

wt_ev <- wt_grid |>
  tidyr::crossing(dose_index = 1:24) |>
  mutate(time = (dose_index - 1) * WEEK, amt = 300, evid = 1L, cmt = "depot") |>
  select(-dose_index) |>
  bind_rows(
    wt_grid |>
      tidyr::crossing(time = WK24 + c(seq(0, 24, by = 0.5), seq(30, WEEK, by = 6))) |>
      mutate(amt = NA_real_, evid = 0L, cmt = "central")
  ) |>
  arrange(id, time, desc(evid))

wt_sim <- rxode2::rxSolve(rxode2::zeroRe(mod), wt_ev, keep = c("WT"),
                          maxsteps = 200000L) |>
  as.data.frame()
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalka', 'etaltlag', 'etalq', 'etalvp', 'etalq2', 'etalvp2'
#> Warning: multi-subject simulation without without 'omega'

wt_tbl <- wt_sim |>
  dplyr::filter(!is.na(Cc)) |>
  group_by(WT) |>
  summarise(
    auctau  = auc_trapz(time, Cc),
    cmax    = max(Cc),
    ctrough = Cc[which.max(time)],
    .groups = "drop"
  ) |>
  left_join(
    tibble::tibble(WT = c(40, 70, 100),
                   pub_auctau = c(131.0, 99.2, 83.2),
                   pub_cmax   = c(12.5, 8.5, 6.6),
                   pub_ctrough = c(0.0383, 0.0225, 0.0160)),
    by = "WT"
  ) |>
  mutate(
    pct_auctau = 100 * (auctau / pub_auctau - 1),
    pct_cmax   = 100 * (cmax / pub_cmax - 1)
  )

wt_tbl |>
  select(WT, pub_auctau, auctau, pct_auctau, pub_cmax, cmax, pct_cmax,
         pub_ctrough, ctrough) |>
  dplyr::rename(
    "Weight (kg)"                 = WT,
    "Published AUCtau (ug*h/mL)"  = pub_auctau,
    "Simulated AUCtau"            = auctau,
    "% diff AUCtau"               = pct_auctau,
    "Published Cmax (ug/mL)"      = pub_cmax,
    "Simulated Cmax"              = cmax,
    "% diff Cmax"                 = pct_cmax,
    "Published Ctrough (ug/mL)"   = pub_ctrough,
    "Simulated Ctrough"           = ctrough
  ) |>
  knitr::kable(digits = c(0, 1, 1, 1, 1, 2, 1, 4, 4),
               caption = "Typical-value steady-state exposures at 300 mg weekly versus Youssef 2024 Table 4. The published values are medians of a 1000-subject cohort, so a few percent of typical-value-versus-median difference is expected.")
Typical-value steady-state exposures at 300 mg weekly versus Youssef 2024 Table 4. The published values are medians of a 1000-subject cohort, so a few percent of typical-value-versus-median difference is expected.
Weight (kg) Published AUCtau (ug*h/mL) Simulated AUCtau % diff AUCtau Published Cmax (ug/mL) Simulated Cmax % diff Cmax Published Ctrough (ug/mL) Simulated Ctrough
40 131.0 130.3 -0.5 12.5 13.36 6.9 0.0383 0.0355
70 99.2 98.9 -0.3 8.5 8.99 5.7 0.0225 0.0208
100 83.2 83.0 -0.2 6.6 6.91 4.7 0.0160 0.0148 

# Deterministic (zeroRe): no cohort draw, so these bounds are tight. Realised
# -0.3% to +0.7% on AUCtau and +1.2% to +5.8% on Cmax.
stopifnot(
  max(abs(wt_tbl$pct_auctau)) < 5,
  max(abs(wt_tbl$pct_cmax))   < 10
)

# The paper's headline conclusion: over the 40-100 kg body-weight range the
# difference in median AUCtau is "<40%", which it judged not clinically
# relevant. Its own Table 4 values give 100 * (1 - 83.2/131.0) = 36.5%.
auc_drop_pct <- 100 * (1 - min(wt_tbl$auctau) / max(wt_tbl$auctau))
cat(sprintf("AUCtau decrease from 40 kg to 100 kg: %.1f%% (Youssef 2024: <40%%; its Table 4 gives 36.5%%)\n",
            auc_drop_pct))
#> AUCtau decrease from 40 kg to 100 kg: 36.3% (Youssef 2024: <40%; its Table 4 gives 36.5%)

# Deterministic: this is fully determined by the weight exponent, so bound it
# on both sides. The paper's own claim caps it at 40%; the lower bound is what
# makes the gate able to go red -- a theoretical 0.75 exponent would give 49.7%
# and no weight effect at all would give 0%.
stopifnot(auc_drop_pct < 40, auc_drop_pct > 30)

Assumptions and deviations

Errata in the source

  1. The abstract is wrong about race. It states that “Chronic HBV infection status, body weight, and Asian versus non-Asian race were key covariates included in the final model.” Race is not in the final model. Table 2 has no race parameter; the Results text says “race was not a statistically significant covariate”; Supplementary Table 4 step 3 shows race on CL/F and on V2/F both failing forward selection (dMVOF +8.15 and +5.32, both p = 1); and Supplementary Table 5 lists exactly four surviving covariate effects. The four real effects are extracted and race is recorded in covariatesDataExcluded. The paper did simulate Asian versus non-Asian exposures post hoc (Table 4), but that contrast is driven entirely by the two groups’ different mean weights acting through the retained weight effect.

  2. Reference weight 73 kg versus the 70.2 kg median. Supplementary Equation 5 states that continuous covariates were centred on their median values and Table 1 reports a median weight of 70.2 kg, yet the Table 2 footnote prints (WT/73). The printed equation is used, per the standing rule that a printed equation outranks prose, and the “Structural check 2” section above shows the paper’s own Table 4 simulations confirm it: a 70 kg reference leaves a systematic 2.4-2.6% negative bias at all three published weights, while 73 kg lands within 0.6%.

  3. Base-model versus final weight exponents. The main-text Methods quotes “0.48 and 0.92” and Supplementary Table 2 quotes 0.477 and 0.92. Those are the base model, fitted before the two disease-status effects entered. The final estimates in Table 2 are 0.494 and 1.01, and those are what the model file uses.

  4. Loading doses are excluded. Supplementary Table 1 footnote b records that phase 2b arms 1-3 received additional loading doses on Days 4 and 11 on top of weekly dosing. Table 3’s simulations are explicitly labelled “no loading dose”, so the regimens simulated here are weekly-only, matching the published table.

  5. Baseline HBsAg category threshold. Supplementary Table 3 describes the covariate as “>1000 IU/mL vs <=1000 IU/mL” while the Supplementary Table 4 footnote describes the same row as “>3 log10 IU/mL vs <=3 log10 IU/mL”. These are the same threshold written two ways. The covariate was rejected in any case.

Known deviation: Ctrough is under-predicted

AUCtau reproduces to within about 1% and Cmax to within about 3%, but Ctrough is consistently 10-27% below the published values, in every regimen and at every body weight, and in the same proportion at 450 mg and at 300 mg. The deviation is therefore structural rather than dose-related or a transcription error in a clearance or volume, both of which would have moved AUCtau too.

Two things make this the most fragile of the three metrics. First, Ctrough is the only one that sits entirely in the terminal phase, which is governed by the deep compartment – whose IIV is the largest in the model (146% CV on Q4/F, 198% on V4/F), so the cohort median moves by more than 10 percentage points between seeds. Second, the paper never states the sampling time behind “trough plasma concentration”; the value used here is the concentration exactly 168 h after the final dose, and the typical-value profile passes through the published 0.0225 ug/mL about 24 h earlier in the interval. A pre-dose visit sample taken slightly short of a full week would account for most of the gap.

This is recorded rather than tuned away. The comparison table above reports the Cmin rows in full, and the gate on Ctrough is a two-sided band: loose enough to admit the known deviation, but tight enough to fail both if the terminal phase breaks outright and if the deviation were ever to disappear silently.

Modelling assumptions

  • Weight distribution. Youssef 2024 reports the median and range of body weight but not its distribution. A log-normal with median 70.2 kg and sdlog = 0.20, truncated to the observed 43.3-140.0 kg range, is used; that gives a spread close to the reported mean (SD) of 71.2 (14.5) kg.
  • Disease status in the weight strata. Table 4’s fixed-weight rows do not state a disease status. DIS_CHB = 1 is assumed, matching the 96% of the analysis population with chronic HBV infection. The assumption is immaterial for AUCtau and Cmax, since the disease effect acts only on the absorption lag time and the shallow peripheral volume, neither of which enters Dose / (CL/F).
  • Cohort size. 200 participants per arm, against the paper’s 1000. Every cohort-derived assertion above is bounded with headroom measured across multiple seeds rather than pinned to a single draw.
  • Apparent parameters. All clearances and volumes are apparent (CL/F, V/F). Bepirovirsen was given only subcutaneously in these three studies, so bioavailability is not identifiable and no f(depot) term is applied; F is folded into the tabulated values exactly as the paper’s own notation states.
  • No time-varying covariates. Body weight is held at its baseline value, as in the source analysis.
  • Every parameter value in the model file comes from Youssef 2024 Table 2 or its printed footnote equations. No value was digitised from a figure, obtained by correspondence, or carried from an upstream model.