Skip to contents

Lecanemab plasma pathophysiology biomarkers (Bhagunde 2026)

Source: vignettes/articles/Bhagunde_2026_lecanemab_biomarkers.Rmd
Bhagunde_2026_lecanemab_biomarkers.Rmd

Model and source

Bhagunde 2026 develops three independent indirect-response PK/PD models, one per plasma biomarker, all driven by the same published lecanemab population PK model. They are packaged as three model files:

model_names <- c(
  "Bhagunde_2026_lecanemab_abeta4240",
  "Bhagunde_2026_lecanemab_ptau181",
  "Bhagunde_2026_lecanemab_gfap"
)
uis <- lapply(model_names, function(n) rxode2::rxode(readModelDb(n)))
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ parameter labels from comments will be replaced by 'label()'
names(uis) <- model_names
  • Citation: Bhagunde P, Penner N, Willis BA, Bell R, Sachdev P, Charil A, Irizarry MC, Hersch S, Reyderman L (2026). Pharmacokinetic/pharmacodynamic analyses of plasma pathophysiology biomarkers in subjects with early Alzheimer’s disease following lecanemab treatment. Alzheimer’s & Dementia: Translational Research & Clinical Interventions 12(2):e70246. doi:10.1002/trc2.70246. Lecanemab population PK driver (CL, V1, V2, Q and their covariate effects) taken from the cited upstream model: Majid O, Cao Y, Willis BA, et al. (2024) Population pharmacokinetics and exposure-response analyses of safety (ARIA-E and isolated ARIA-H) of lecanemab in subjects with early Alzheimer’s disease. CPT Pharmacometrics Syst Pharmacol 13(12):2111-2123. doi:10.1002/psp4.13224.
  • Article: https://doi.org/10.1002/trc2.70246
  • Supplement (Appendix A, Tables S1-S5, Figures S1-S7): https://www.ebi.ac.uk/europepmc/webservices/rest/PMC13095857/supplementaryFiles
  • Upstream lecanemab population PK (Majid 2024): https://doi.org/10.1002/psp4.13224
  • Upstream amyloid-PET plaque model (Bhagunde 2026, CPT:PSP): https://doi.org/10.1002/psp4.70173

Bhagunde_2026_lecanemab_abeta4240

Indirect-response PK/PD model for the absolute plasma amyloid-beta 42/40 ratio in subjects with early Alzheimer’s disease receiving the anti-protofibril monoclonal antibody lecanemab. Serum lecanemab concentration Cc drives a linear stimulation of the Abeta42/40 production rate: dR/dt = Kin * (1 + SLOPE * Cc) - Kout * R, with Kin = baseline * Kout so the untreated pool sits at its baseline. Covariate effects (all log-linear) are age, APOE4-carrier status, Japanese race and female sex on the baseline ratio, and age and body weight on SLOPE. Interindividual variability is exponential on baseline and SLOPE with an estimated correlation, and residual error is proportional. Fit by NONMEM FOCE-I to 12,468 plasma Abeta42/40 observations from 1994 subjects pooled across the lecanemab phase 2 Study 201 (Core and open-label extension) and phase 3 Clarity AD (Study 301, Core and open-label extension) trials. The lecanemab serum-concentration driver is the two-compartment population PK model of Majid 2024 (zero-order IV input, linear elimination from the central compartment); its parameters are held fixed here because the present paper did not re-estimate them.

Bhagunde_2026_lecanemab_ptau181

Indirect-response PK/PD model for absolute plasma tau phosphorylated at threonine 181 (p-tau181) in subjects with early Alzheimer’s disease receiving the anti-protofibril monoclonal antibody lecanemab. Serum lecanemab concentration Cc inhibits p-tau181 production through an Emax function: dR/dt = Kin * (1 - Emax * Cc / (EC50 + Cc)) - Kout * R, with Kin = baseline * Kout so the untreated pool sits at its baseline. Covariate effects act on the baseline only: power terms for body weight and baseline Mini-Mental State Examination score, and a ratio for APOE4-carrier status. Interindividual variability is exponential on all four PD parameters (baseline, Emax, EC50, Kout) with an estimated baseline ~ Emax correlation, and residual error is proportional. Fit by NONMEM FOCE-I to 7909 plasma p-tau181 observations from 2179 subjects pooled across the lecanemab phase 2 Study 201 (Core and open-label extension) and phase 3 Clarity AD (Study 301, Core and open-label extension) trials. The lecanemab serum-concentration driver is the two-compartment population PK model of Majid 2024 (zero-order IV input, linear elimination from the central compartment); its parameters are held fixed here because the present paper did not re-estimate them.

Bhagunde_2026_lecanemab_gfap

Indirect-response PK/PD model with a disease-progression term for plasma glial fibrillary acidic protein (GFAP) in subjects with early Alzheimer’s disease receiving the anti-protofibril monoclonal antibody lecanemab. Unlike the companion Abeta42/40 and p-tau181 models, the drug effect on GFAP is driven not by lecanemab concentration but by the relative reduction in brain amyloid plaque from baseline, which had the lowest objective function value of the structures tested: dGFAP/dt = KinG * (1 - SLP * (Plaque_baseline - Plaque(t)) / Plaque_baseline + DP * TIME) - KoutG * GFAP, with KoutG = KinG / BGFAP at baseline. The amyloid-plaque state is the semi-mechanistic amyloid-PET turnover model of Bhagunde 2026 (CPT:PSP), in which lecanemab stimulates plaque elimination linearly in serum concentration, and lecanemab serum concentration comes from the two-compartment population PK model of Majid 2024. Both upstream layers are held fixed because the present paper did not re-estimate them. Covariates on baseline GFAP are power terms for age, body weight and baseline Mini-Mental State Examination score plus a study-effect term for Study 201 (an assay batch effect). Variability is exponential on BGFAP and KinG, additive on the logit of SLP, and residual error is proportional. Fit by Monolix SAEM to 3098 plasma GFAP observations from 736 subjects pooled across the lecanemab phase 2 Study 201 and phase 3 Clarity AD (Study 301) trials.

Population

Bhagunde 2026 pooled the Core and open-label-extension parts of the lecanemab phase 2 Study 201 and the phase 3 Clarity AD trial (Study 301) in subjects with early Alzheimer’s disease (mild cognitive impairment due to Alzheimer’s disease or mild Alzheimer’s disease dementia). Study 201 randomised 856 subjects across placebo and five lecanemab regimens for an 18-month Core phase, followed by an off-treatment gap of 9-59 months (24 months on average) and then open-label 10 mg/kg biweekly; Clarity AD randomised 1:1 to placebo or lecanemab 10 mg/kg biweekly for 18 months with an open-label extension (Bhagunde 2026 Appendix A).

Each biomarker model was fit to a different subset (Bhagunde 2026 Table S1):

Analysis populations (Bhagunde 2026 Table S1).
Model Biomarker N subjects N observations
Bhagunde_2026_lecanemab_abeta4240 plasma Abeta42/40 ratio 1994 12468
Bhagunde_2026_lecanemab_ptau181 plasma p-tau181 2179 7909
Bhagunde_2026_lecanemab_gfap plasma GFAP 736 3098

The GFAP dataset is the smallest because it was restricted to subjects who also had amyloid PET data – the GFAP drug effect is driven by relative change in amyloid plaque rather than by lecanemab concentration.

The same information is available programmatically via each model’s population metadata (rxode2::rxode(readModelDb("<model>"))$population).

Source trace

The per-parameter origin is recorded as an in-file comment next to each ini() entry in inst/modeldb/specificDrugs/Bhagunde_2026_lecanemab_*.R. The table below collects the equations and parameters in one place.

Equation / parameter Value Source location
Lecanemab PK (all three models; held fixed())
d/dt(central), d/dt(peripheral1) two-compartment, zero-order IV input, linear elimination n/a Bhagunde 2026 Section 2.2.1; structure and values from Majid 2024
lcl (CL) 0.0154 L/h Majid 2024 Table 1
lvc (V1) 3.24 L Majid 2024 Table 1
lvp (V2) 2.00 L Majid 2024 Table 1
lq (Q) 0.00718 L/h Majid 2024 Table 1
e_wt_cl, e_alb_cl, e_sexf_cl, e_ada_pos_cl 0.353, -0.374, 0.791, 1.13 Majid 2024 Table 1 + CL = 0.0154*(BW/72)^0.353*(ALB/43)^-0.374*0.791^SEX*1.13^ADA
e_wt_vc, e_sexf_vc, e_race_japanese_vc 0.513, 0.868, 0.920 Majid 2024 Table 1 + V1 = 3.24*(BW/72)^0.513*0.868^SEX*0.920^JPN
e_race_japanese_vp 0.671 Majid 2024 Table 1 + V2 = 2.00*0.671^JPN
etalcl, etalvc (block), etalvp 34.9 %CV, 12.2 %CV, R = 0.144; 94.6 %CV Majid 2024 Table 1
Plasma Abeta42/40 model
dR/dt = Kin*(1 + SLOPE*Conc) - R*Kout n/a Bhagunde 2026 Section 2.2.2, first displayed equation
lrbase (baseline ratio) 0.0864 Bhagunde 2026 Table 1
lkout (Kout) 1.51 /year Bhagunde 2026 Table 1
lslope (SLOPE) 0.000704 per ug/mL Bhagunde 2026 Table 1
e_age_rbase, e_apoe4_carrier_rbase, e_race_japanese_rbase, e_sexf_rbase -0.00181, 0.0138, 0.0346, -0.00888 Bhagunde 2026 Table 1 (covariate effects on baseline)
e_age_slope, e_wt_slope 0.0193, -0.00602 Bhagunde 2026 Table 1 (covariate effects on slope)
etalrbase + etalslope block 7.21 %CV, 71.7 %CV, R = -0.573 Bhagunde 2026 Table 1
propSd 7.72 %CV Bhagunde 2026 Table 1
Plasma p-tau181 model
dR/dt = Kin*(1 - Emax*Conc/(EC50+Conc)) - R*Kout n/a Bhagunde 2026 Section 2.2.2, second displayed equation
lrbase, lemax, lec50, lkout 3.40 pg/mL, 0.480, 31.4 ug/mL, 0.428 /year Bhagunde 2026 Table 1
e_wt_rbase, e_mmse_rbase, e_apoe4_carrier_rbase -0.235, -0.468, 1.07 Bhagunde 2026 Table 1
etalrbase + etalemax block; etalec50; etalkout 39.7 / 28.0 %CV, R = 0.628; 252 %CV; 63.2 %CV Bhagunde 2026 Table 1
propSd 19.2 %CV Bhagunde 2026 Table 1
Plasma GFAP model
dGFAP/dt = KinG*(1 - SLP*(Plaque0-Plaque)/Plaque0 + DP*TIME) - KoutG*GFAP n/a Bhagunde 2026 Section 2.2.2, third displayed equation
KinG = BGFAP * KoutG at baseline n/a Bhagunde 2026 Section 2.2.2 text
dPlaque/dt = Kin - Plaque*Kout*(1 + DESLP*Conc) n/a Bhagunde 2026 Section 2.2.2, fourth displayed equation
lrbase (BGFAP), lkin (KinG), logitslope_gfap (SLP), lprog (DP) 315 pg/mL, 5.6 pg/mL/h, 0.237, 0.0307 /year Bhagunde 2026 Table 1
e_age_rbase, e_wt_rbase, e_mmse_rbase, e_study_lec201_rbase 1.03, -0.66, -0.63, 0.388 Bhagunde 2026 Table 1
etalrbase, etalkin, etalogitslope_gfap 38.2 %CV, 133 %CV, SD 0.931 on logit Bhagunde 2026 Table 1
propSd 0.193 (SD) Bhagunde 2026 Table 1
logitrbase_plaque, lkout_plaque, lslope_plaque (held fixed()) -1.12 (logit), 0.0572 /year, 0.154 per ug/mL Bhagunde 2026 CPT:PSP Table 1
e_apoe4_carrier_logitrbase_plaque, e_age_slope_plaque (held fixed()) 0.629, 3.10 Bhagunde 2026 CPT:PSP Table 1
Plaque logit scale BSL = 250 * plogis(phi) 250 CL max Bhagunde 2026 CPT:PSP Section 2.4.1

Simulation helpers 

HR_YR  <- 8766          # hours per year, as used inside model()
HR_MO  <- HR_YR / 12    # hours per month
DOSE_MGKG <- 10         # LEC10 regimens

# Reference subject: the centring values used by every covariate term.
cov_ref <- list(
  WT = 72, ALB = 43, SEXF = 0, ADA_POS = 0, RACE_JAPANESE = 0,
  AGE = 72, APOE4_CARRIER = 0, SCORE_MMSE = 26, STUDY_LEC201 = 0
)

# Build one subject's event table. `switch_mo = Inf` keeps biweekly dosing for
# the whole horizon; a finite value switches to monthly (LEC10-M) at that month;
# `stop_mo` discontinues dosing entirely.
make_events <- function(id, covs, horizon_mo = 48, switch_mo = Inf,
                        stop_mo = Inf, obs_step_mo = 1) {
  wt <- covs$WT
  dose_times <- numeric(0)
  t_h <- 0
  end_h <- horizon_mo * HR_MO
  repeat {
    if (t_h >= end_h || t_h >= stop_mo * HR_MO) break
    dose_times <- c(dose_times, t_h)
    ii <- if (t_h >= switch_mo * HR_MO) 24 * 28 else 24 * 14
    t_h <- t_h + ii
  }
  obs_times <- seq(0, end_h, by = obs_step_mo * HR_MO)
  ev <- dplyr::bind_rows(
    tibble::tibble(time = dose_times, amt = DOSE_MGKG * wt, evid = 1L,
                   dur = 1, cmt = "central"),
    tibble::tibble(time = obs_times, amt = NA_real_, evid = 0L,
                   dur = NA_real_, cmt = "central")
  ) |>
    dplyr::arrange(time, dplyr::desc(evid)) |>
    dplyr::mutate(id = id)
  for (nm in names(covs)) ev[[nm]] <- covs[[nm]]
  ev
}

# Typical-value (zeroRe) solve of one arm, returned as a data frame with `id`
# restored (rxSolve drops it for a single subject).
solve_typical <- function(model_name, ev) {
  out <- rxode2::rxSolve(
    rxode2::zeroRe(readModelDb(model_name)), ev, returnType = "data.frame"
  )
  if (is.null(out$id)) out$id <- ev$id[1]
  out
}

at_month <- function(df, value_col, months) {
  d <- df[!is.na(df[[value_col]]), c("time", value_col)]
  d <- d[!duplicated(d$time), ]
  d <- d[order(d$time), ]
  stats::approx(d$time, d[[value_col]], xout = months * HR_MO)$y
}

first_value <- function(df, value_col) {
  d <- df[!is.na(df[[value_col]]), c("time", value_col)]
  d[[value_col]][which.min(d$time)]
}

Validation 1 – the lecanemab PK layer

The PK layer is inherited from Majid 2024 and is the driver of all three PD models, so it is validated first with a per-subject NCA identity: for a linear two-compartment model with IV input, AUCinf must equal Dose / CL for every subject. cl is returned per subject by rxSolve, so the check is exact rather than a comparison against a summary statistic.

set.seed(20260901)
n_pk <- 50
pk_cohort <- tibble::tibble(
  id  = seq_len(n_pk),
  WT  = round(stats::rlnorm(n_pk, log(72), 0.20), 1),
  ALB = round(stats::rnorm(n_pk, 43, 3), 1),
  SEXF = stats::rbinom(n_pk, 1, 0.52),
  ADA_POS = 0,
  RACE_JAPANESE = stats::rbinom(n_pk, 1, 0.06),
  AGE = round(stats::rnorm(n_pk, 72, 8)),
  APOE4_CARRIER = stats::rbinom(n_pk, 1, 0.69),
  SCORE_MMSE = pmin(30, pmax(20, round(stats::rnorm(n_pk, 25.5, 2)))),
  STUDY_LEC201 = 0
)

# Single 10 mg/kg 1-h infusion; a fine early grid resolves the distribution
# phase, and the window is truncated at 1680 h (about 7 terminal half-lives) so
# lambda-z is fitted on real signal rather than solver noise.
pk_grid <- c(seq(0, 24, by = 0.5), seq(30, 168, by = 6), seq(180, 1680, by = 24))
pk_events <- dplyr::bind_rows(lapply(seq_len(n_pk), function(i) {
  covs <- as.list(pk_cohort[i, setdiff(names(pk_cohort), "id")])
  ev <- dplyr::bind_rows(
    tibble::tibble(time = 0, amt = DOSE_MGKG * covs$WT, evid = 1L, dur = 1,
                   cmt = "central"),
    tibble::tibble(time = pk_grid, amt = NA_real_, evid = 0L, dur = NA_real_,
                   cmt = "central")
  ) |>
    dplyr::arrange(time, dplyr::desc(evid)) |>
    dplyr::mutate(id = pk_cohort$id[i])
  for (nm in names(covs)) ev[[nm]] <- covs[[nm]]
  ev
}))
stopifnot(!anyDuplicated(unique(pk_events[, c("id", "time", "evid")])))

pk_sim <- rxode2::rxSolve(
  readModelDb("Bhagunde_2026_lecanemab_abeta4240"), pk_events,
  keep = c("WT"), returnType = "data.frame"
)
#> ℹ parameter labels from comments will be replaced by 'label()'
sim_nca <- pk_sim |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::mutate(treatment = "LEC 10 mg/kg single dose") |>
  dplyr::select(id, time, Cc, treatment)

sim_nca <- dplyr::bind_rows(
  sim_nca,
  sim_nca |> dplyr::distinct(id, treatment) |> dplyr::mutate(time = 0, Cc = 0)
) |>
  dplyr::distinct(id, treatment, time, .keep_all = TRUE) |>
  dplyr::arrange(id, treatment, time)
stopifnot(all(sim_nca$Cc >= 0))

conc_obj <- PKNCA::PKNCAconc(sim_nca, Cc ~ time | treatment + id)

dose_df <- pk_events |>
  dplyr::filter(evid == 1) |>
  dplyr::mutate(treatment = "LEC 10 mg/kg single dose") |>
  dplyr::select(id, time, amt, treatment)
dose_obj <- PKNCA::PKNCAdose(dose_df, amt ~ time | treatment + id)

intervals <- data.frame(
  start = 0, end = Inf,
  cmax = TRUE, tmax = TRUE, aucinf.obs = TRUE, half.life = TRUE
)
nca_res <- PKNCA::pk.nca(PKNCA::PKNCAdata(conc_obj, dose_obj, intervals = intervals))

nca_wide <- as.data.frame(nca_res) |>
  dplyr::filter(PPTESTCD %in% c("cmax", "tmax", "aucinf.obs", "half.life")) |>
  dplyr::select(id, PPTESTCD, PPORRES) |>
  tidyr::pivot_wider(names_from = PPTESTCD, values_from = PPORRES)

# Per-subject exact identity: AUCinf == Dose / CL.
per_subject_cl <- pk_sim |> dplyr::group_by(id) |> dplyr::summarise(cl = dplyr::first(cl), .groups = "drop")
ident <- nca_wide |>
  dplyr::left_join(per_subject_cl, by = "id") |>
  dplyr::left_join(dose_df |> dplyr::select(id, amt), by = "id") |>
  dplyr::mutate(auc_theory = amt / cl,
                pct_diff = 100 * (aucinf.obs - auc_theory) / auc_theory)
stopifnot(nrow(ident) == n_pk, !anyNA(ident$pct_diff))

tibble::tibble(
  Check = "AUCinf vs Dose/CL (per subject)",
  `Median % difference` = round(stats::median(ident$pct_diff), 3),
  `90th pct |% difference|` = round(stats::quantile(abs(ident$pct_diff), 0.9), 3),
  `Max |% difference|` = round(max(abs(ident$pct_diff)), 3)
) |>
  knitr::kable(caption = "Structural identity check on the inherited lecanemab PK layer.")
Structural identity check on the inherited lecanemab PK layer.
Check Median % difference 90th pct |% difference| Max |% difference|
AUCinf vs Dose/CL (per subject) -0.008 0.257 2.055 

stopifnot(
  abs(stats::median(ident$pct_diff)) < 1,
  stats::quantile(abs(ident$pct_diff), 0.9) < 2
)

The identity holds to well under a percent, confirming that the two-compartment disposition, the covariate equations and the dose bookkeeping were transcribed correctly.

Validation 2 – typical-value biomarker trajectories

Bhagunde 2026 Section 4 states the model-predicted typical-value trajectories explicitly, which gives a direct numeric answer key for all three biomarkers:

“Plasma Abeta42/40 levels increase from ~0.087 at baseline to ~0.095 after 18 months of treatment … Plasma p-tau181 levels decrease from 3.37 pg/mL at baseline to 2.80 pg/mL after 18 months … plasma GFAP levels decrease from 316 pg/mL at baseline to 275 pg/mL after 18 months.”

arms <- tibble::tribble(
  ~arm,                    ~switch_mo, ~stop_mo,
  "LEC10-BW (continuous)",       Inf,      Inf,
  "LEC10-M from month 18",        18,      Inf,
  "LEC10-M from month 24",        24,      Inf,
  "LEC10-M from month 30",        30,      Inf,
  "Discontinue at month 18",     Inf,       18
)

state_of <- c(Bhagunde_2026_lecanemab_abeta4240 = "abeta4240",
              Bhagunde_2026_lecanemab_ptau181   = "ptau181",
              Bhagunde_2026_lecanemab_gfap      = "gfap")

typ <- dplyr::bind_rows(lapply(seq_len(nrow(arms)), function(a) {
  ev <- make_events(id = a, covs = cov_ref, horizon_mo = 48,
                    switch_mo = arms$switch_mo[a], stop_mo = arms$stop_mo[a])
  dplyr::bind_rows(lapply(model_names, function(mn) {
    s <- solve_typical(mn, ev)
    tibble::tibble(
      arm = arms$arm[a], model = mn, time = s$time, Cc = s$Cc,
      value = s[[state_of[[mn]]]],
      plaque = if (mn == "Bhagunde_2026_lecanemab_gfap") s$plaque else NA_real_
    )
  }))
}))
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalrbase', 'etalslope'
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalrbase', 'etalemax', 'etalec50', 'etalkout'
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalrbase', 'etalkin', 'etalogitslope_gfap'
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalrbase', 'etalslope'
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalrbase', 'etalemax', 'etalec50', 'etalkout'
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalrbase', 'etalkin', 'etalogitslope_gfap'
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalrbase', 'etalslope'
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalrbase', 'etalemax', 'etalec50', 'etalkout'
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalrbase', 'etalkin', 'etalogitslope_gfap'
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalrbase', 'etalslope'
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalrbase', 'etalemax', 'etalec50', 'etalkout'
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalrbase', 'etalkin', 'etalogitslope_gfap'
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalrbase', 'etalslope'
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalrbase', 'etalemax', 'etalec50', 'etalkout'
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalrbase', 'etalkin', 'etalogitslope_gfap'

bw <- dplyr::filter(typ, arm == "LEC10-BW (continuous)")

answer_key <- tibble::tribble(
  ~model,                              ~Biomarker,                ~Units,  ~`Baseline (paper)`, ~`Month 18 (paper)`,
  "Bhagunde_2026_lecanemab_abeta4240", "Plasma Abeta42/40 ratio", "ratio", 0.087,               0.095,
  "Bhagunde_2026_lecanemab_ptau181",   "Plasma p-tau181",         "pg/mL", 3.37,                2.80,
  "Bhagunde_2026_lecanemab_gfap",      "Plasma GFAP",             "pg/mL", 316,                 275
)

bw_by_model <- split(bw, bw$model)
traj_chk <- answer_key |>
  dplyr::mutate(
    `Baseline (model)` = vapply(model, function(m) first_value(bw_by_model[[m]], "value"), numeric(1)),
    `Month 18 (model)` = vapply(model, function(m) at_month(bw_by_model[[m]], "value", 18), numeric(1))
  ) |>
  dplyr::mutate(
    `Baseline % diff` = 100 * (`Baseline (model)` - `Baseline (paper)`) / `Baseline (paper)`,
    `Month 18 % diff` = 100 * (`Month 18 (model)` - `Month 18 (paper)`) / `Month 18 (paper)`
  )

traj_chk |>
  dplyr::select(Biomarker, Units, `Baseline (paper)`, `Baseline (model)`,
                `Baseline % diff`, `Month 18 (paper)`, `Month 18 (model)`,
                `Month 18 % diff`) |>
  knitr::kable(
    digits = c(NA, NA, 4, 4, 1, 4, 4, 1),
    caption = paste(
      "Typical-value (zeroRe) reference subject vs the trajectories stated in",
      "Bhagunde 2026 Section 4."
    )
  )
Typical-value (zeroRe) reference subject vs the trajectories stated in Bhagunde 2026 Section 4.
Biomarker Units Baseline (paper) Baseline (model) Baseline % diff Month 18 (paper) Month 18 (model) Month 18 % diff
Plasma Abeta42/40 ratio ratio 0.087 0.0864 -0.7 0.095 0.0940 -1.1
Plasma p-tau181 pg/mL 3.370 3.4000 0.9 2.800 2.7883 -0.4
Plasma GFAP pg/mL 316.000 315.0000 -0.3 275.000 268.9202 -2.2 

stopifnot(
  # Baselines are the tabulated typical values; the paper's Section 4 numbers
  # are cohort-level and round differently, so allow 1.5%.
  max(abs(traj_chk$`Baseline % diff`)) < 1.5,
  # The month-18 value is the load-bearing prediction: it combines the PK layer,
  # the drug-effect function and Kout. A mis-transcribed CL, EC50, SLOPE or
  # Kout moves it by tens of percent.
  max(abs(traj_chk$`Month 18 % diff`)) < 5
)

Degradation half-lives

The paper derives a half-life from each biomarker’s Kout. These are exact identities on the ini() values, so they are asserted tightly.

kout_ab <- exp(uis[["Bhagunde_2026_lecanemab_abeta4240"]]$theta[["lkout"]])
kout_pt <- exp(uis[["Bhagunde_2026_lecanemab_ptau181"]]$theta[["lkout"]])
hl <- tibble::tibble(
  Biomarker = c("Plasma Abeta42/40 ratio", "Plasma p-tau181"),
  `Kout (1/year)` = c(kout_ab, kout_pt),
  `t1/2 (years, model)` = log(2) / c(kout_ab, kout_pt),
  `t1/2 (years, paper)` = c(0.5, 1.6),
  Source = c("Section 3.2.1", "Section 3.2.2")
)
knitr::kable(hl, digits = 3, caption = "Biomarker degradation half-lives.")
Biomarker degradation half-lives.
Biomarker Kout (1/year) t1/2 (years, model) t1/2 (years, paper) Source
Plasma Abeta42/40 ratio 1.510 0.459 0.5 Section 3.2.1
Plasma p-tau181 0.428 1.620 1.6 Section 3.2.2 
stopifnot(abs(log(2) / kout_ab - 0.5) < 0.05, abs(log(2) / kout_pt - 1.6) < 0.05)

Re-accumulation after discontinuation

Section 4 also quantifies the first-year rate of return toward baseline after stopping lecanemab at 18 months. This exercises a different part of each model (the unstimulated turnover, and for GFAP the plaque re-accumulation plus the disease-progression drift).

disc <- dplyr::filter(typ, arm == "Discontinue at month 18")

reacc <- tibble::tribble(
  ~model,                              ~Biomarker,                ~Units,           ~`Rate (paper)`, ~Direction,
  "Bhagunde_2026_lecanemab_abeta4240", "Plasma Abeta42/40 ratio", "ratio/year",     0.006,           "decrease",
  "Bhagunde_2026_lecanemab_ptau181",   "Plasma p-tau181",         "pg/mL per year", 0.19,            "increase",
  "Bhagunde_2026_lecanemab_gfap",      "Plasma GFAP",             "pg/mL per year", 12,              "increase"
)
disc_by_model <- split(disc, disc$model)
reacc <- reacc |>
  dplyr::mutate(
    `Rate (model)` = vapply(model, function(m) {
      d <- disc_by_model[[m]]
      abs(at_month(d, "value", 30) - at_month(d, "value", 18))
    }, numeric(1)),
    `% diff` = 100 * (`Rate (model)` - `Rate (paper)`) / `Rate (paper)`
  )

reacc |>
  dplyr::select(Biomarker, Units, Direction, `Rate (paper)`, `Rate (model)`, `% diff`) |>
  knitr::kable(digits = c(NA, NA, NA, 4, 4, 1),
               caption = "First-year change after stopping lecanemab at month 18 (Bhagunde 2026 Section 4).")
First-year change after stopping lecanemab at month 18 (Bhagunde 2026 Section 4).
Biomarker Units Direction Rate (paper) Rate (model) % diff
Plasma Abeta42/40 ratio ratio/year decrease 0.006 0.0057 -4.7
Plasma p-tau181 pg/mL per year increase 0.190 0.1726 -9.2
Plasma GFAP pg/mL per year increase 12.000 12.0537 0.4 

stopifnot(max(abs(reacc$`% diff`)) < 25)

All three first-year re-accumulation rates land within a quarter of the published values. The looser tolerance here reflects that Section 4 reports these rates to one or two significant figures.

Validation 3 – the amyloid-plaque layer inside the GFAP model

The GFAP drug effect is driven by relative plaque reduction, so the plaque ODE is validated separately against its own source paper, which reports the mean decrease from baseline in amyloid PET at 18, 30 and 42 months of LEC10-BW. Because baseline plaque depends strongly on APOE4-carrier status (61.5 CL for non-carriers, 82.7 CL for carriers), both are shown; the published cohort mean is a mixture of the two.

plaque_arms <- dplyr::bind_rows(lapply(c(0, 1), function(ap) {
  ev <- make_events(id = 100 + ap, covs = modifyList(cov_ref, list(APOE4_CARRIER = ap)),
                    horizon_mo = 48)
  s <- solve_typical("Bhagunde_2026_lecanemab_gfap", ev)
  tibble::tibble(APOE4 = ap, time = s$time, plaque = s$plaque)
}))
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalrbase', 'etalkin', 'etalogitslope_gfap'
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalrbase', 'etalkin', 'etalogitslope_gfap'

plaque_by_apoe <- split(plaque_arms, plaque_arms$APOE4)
plaque_tab <- tidyr::expand_grid(APOE4 = c(0, 1), Month = c(18, 30, 42)) |>
  dplyr::mutate(
    `Baseline (CL)` = vapply(as.character(APOE4),
                             function(a) first_value(plaque_by_apoe[[a]], "plaque"),
                             numeric(1)),
    `Decrease (CL, model)` = mapply(function(a, mo) {
      d <- plaque_by_apoe[[as.character(a)]]
      first_value(d, "plaque") - at_month(d, "plaque", mo)
    }, APOE4, Month)
  ) |>
  dplyr::left_join(
    tibble::tibble(Month = c(18, 30, 42), `Decrease (CL, paper)` = c(59.9, 67, 71.9)),
    by = "Month"
  ) |>
  dplyr::mutate(`APOE4 carrier` = ifelse(APOE4 == 1, "yes", "no")) |>
  dplyr::select(`APOE4 carrier`, Month, `Baseline (CL)`, `Decrease (CL, model)`,
                `Decrease (CL, paper)`)

knitr::kable(plaque_tab, digits = 1,
             caption = paste(
               "Amyloid plaque decrease from baseline on LEC10-BW vs the mean",
               "values reported in Bhagunde 2026 (CPT:PSP) Section 3.1."
             ))
Amyloid plaque decrease from baseline on LEC10-BW vs the mean values reported in Bhagunde 2026 (CPT:PSP) Section 3.1.
APOE4 carrier Month Baseline (CL) Decrease (CL, model) Decrease (CL, paper)
no 18 61.5 50.0 59.9
no 30 61.5 56.3 67.0
no 42 61.5 58.1 71.9
yes 18 82.7 67.2 59.9
yes 30 82.7 75.7 67.0
yes 42 82.7 78.1 71.9 

# The published cohort mean must fall between the two genotype strata at every
# reported time point. This is a genuine two-sided bracket, not a one-sided
# bound: it fails if either the plaque turnover or the DESLP drug effect is
# mis-transcribed in either direction.
bracket <- plaque_tab |>
  dplyr::group_by(Month) |>
  dplyr::summarise(lo = min(`Decrease (CL, model)`), hi = max(`Decrease (CL, model)`),
                   ref = dplyr::first(`Decrease (CL, paper)`), .groups = "drop")
stopifnot(all(bracket$ref >= bracket$lo, bracket$ref <= bracket$hi))

# The typical baseline plaque values printed in the upstream paper (62 and 83
# Centiloids) are exact consequences of the logit-scale baseline and the APOE4
# ratio, so they are asserted tightly.
b0 <- plaque_tab |> dplyr::distinct(`APOE4 carrier`, `Baseline (CL)`)
stopifnot(
  abs(b0$`Baseline (CL)`[b0$`APOE4 carrier` == "no"] - 62) < 1,
  abs(b0$`Baseline (CL)`[b0$`APOE4 carrier` == "yes"] - 83) < 1
)

Validation 4 – covariate forest plots (Figure 4)

Bhagunde 2026 Figure 4 reports fold changes relative to the reference subject for each retained covariate. Panel A is the fold change in Abeta42/40 change from baseline at 18 months; panel C is the fold change in the absolute month-18 GFAP level.

ab_scen <- tibble::tibble(
  Covariate = c("APOE4", "Race", "Sex", "Age", "Age",
                "Baseline weight", "Baseline weight"),
  Level     = c("carrier", "Japanese", "female", "83 years", "57 years",
                "97.5 kg", "48.2 kg"),
  override  = list(list(APOE4_CARRIER = 1), list(RACE_JAPANESE = 1), list(SEXF = 1),
                   list(AGE = 83), list(AGE = 57), list(WT = 97.5), list(WT = 48.2)),
  `Fold change (Fig 4A)` = c(0.99, 1.03, 1.01, 1.21, 0.77, 1.18, 0.78)
)

cfb18 <- function(model_name, override) {
  ev <- make_events(id = 1, covs = modifyList(cov_ref, override), horizon_mo = 20)
  s <- solve_typical(model_name, ev)
  st <- state_of[[model_name]]
  at_month(s, st, 18) - s[[st]][which.min(s$time)]
}
ref_cfb <- cfb18("Bhagunde_2026_lecanemab_abeta4240", list())
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalrbase', 'etalslope'

ab_forest <- ab_scen |>
  dplyr::mutate(`Fold change (model)` = vapply(
    override,
    function(o) cfb18("Bhagunde_2026_lecanemab_abeta4240", o) / ref_cfb,
    numeric(1)
  )) |>
  dplyr::mutate(Difference = `Fold change (model)` - `Fold change (Fig 4A)`) |>
  dplyr::select(Covariate, Level, `Fold change (Fig 4A)`, `Fold change (model)`, Difference)
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalrbase', 'etalslope'
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalrbase', 'etalslope'
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalrbase', 'etalslope'
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalrbase', 'etalslope'
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalrbase', 'etalslope'
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalrbase', 'etalslope'
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalrbase', 'etalslope'

knitr::kable(ab_forest, digits = 3,
             caption = "Abeta42/40 change-from-baseline fold changes vs Bhagunde 2026 Figure 4A.")
Abeta42/40 change-from-baseline fold changes vs Bhagunde 2026 Figure 4A.
Covariate Level Fold change (Fig 4A) Fold change (model) Difference
APOE4 carrier 0.99 1.014 0.024
Race Japanese 1.03 1.037 0.007
Sex female 1.01 1.252 0.242
Age 83 years 1.21 1.212 0.002
Age 57 years 0.77 0.769 -0.001
Baseline weight 97.5 kg 1.18 1.044 -0.136
Baseline weight 48.2 kg 0.78 0.890 0.110

The age, race and APOE4 rows reproduce the published panel to within 0.03, which confirms both the tabulated baseline / slope coefficients and the 72-year age centring the paper does not print. Sex and body weight – the only two covariates in this panel that also act on lecanemab PK – do not reproduce, and no single reading of the panel reconciles them; the discrepancy is analysed under Assumptions and deviations. Both groups are asserted below: the agreeing rows against the published values, the discrepant rows against the model’s own values, so neither can drift silently.

pk_free <- dplyr::filter(ab_forest, !Covariate %in% c("Sex", "Baseline weight"))
stopifnot(nrow(pk_free) == 4, max(abs(pk_free$Difference)) < 0.05)

# Sex and body weight also act through the PK layer; the published panel does
# not reproduce (see Assumptions and deviations). Pin the model's own values so
# the known deviation cannot silently change.
pk_bound <- dplyr::filter(ab_forest, Covariate %in% c("Sex", "Baseline weight"))
expected <- c(female = 1.252, `97.5 kg` = 1.044, `48.2 kg` = 0.890)
stopifnot(nrow(pk_bound) == 3,
          all(abs(pk_bound$`Fold change (model)` - expected[pk_bound$Level]) < 0.02))
gfap_scen <- tibble::tibble(
  Covariate = c("Age", "Age", "Baseline MMSE", "Baseline MMSE", "Weight", "Weight"),
  Level     = c("83 years", "57 years", "29 units", "22 units", "50 kg", "100 kg"),
  override  = list(list(AGE = 83), list(AGE = 57), list(SCORE_MMSE = 29),
                   list(SCORE_MMSE = 22), list(WT = 50), list(WT = 100)),
  `Fold change (Fig 4C)` = c(1.12, 0.836, 0.929, 1.11, 1.27, 0.805)
)

gfap18 <- function(override) {
  ev <- make_events(id = 1, covs = modifyList(cov_ref, override), horizon_mo = 20)
  at_month(solve_typical("Bhagunde_2026_lecanemab_gfap", ev), "gfap", 18)
}
ref_gfap18 <- gfap18(list())
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalrbase', 'etalkin', 'etalogitslope_gfap'

gfap_forest <- gfap_scen |>
  dplyr::mutate(`Fold change (model)` =
                  vapply(override, function(o) gfap18(o) / ref_gfap18, numeric(1))) |>
  dplyr::mutate(Difference = `Fold change (model)` - `Fold change (Fig 4C)`) |>
  dplyr::select(Covariate, Level, `Fold change (Fig 4C)`, `Fold change (model)`, Difference)
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalrbase', 'etalkin', 'etalogitslope_gfap'
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalrbase', 'etalkin', 'etalogitslope_gfap'
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalrbase', 'etalkin', 'etalogitslope_gfap'
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalrbase', 'etalkin', 'etalogitslope_gfap'
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalrbase', 'etalkin', 'etalogitslope_gfap'
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalrbase', 'etalkin', 'etalogitslope_gfap'

knitr::kable(gfap_forest, digits = 3,
             caption = "Month-18 GFAP fold changes vs Bhagunde 2026 Figure 4C.")
Month-18 GFAP fold changes vs Bhagunde 2026 Figure 4C.
Covariate Level Fold change (Fig 4C) Fold change (model) Difference
Age 83 years 1.120 1.124 0.004
Age 57 years 0.836 0.840 0.004
Baseline MMSE 29 units 0.929 0.934 0.005
Baseline MMSE 22 units 1.110 1.111 0.001
Weight 50 kg 1.270 1.299 0.029
Weight 100 kg 0.805 0.792 -0.013 

# MMSE acts only on baseline GFAP through the tabulated power term, so it
# recovers the panel almost exactly and pins the MMSE = 26 centring value the
# paper does not print. Weight and age additionally move lecanemab exposure and
# therefore plaque removal, so they are allowed a slightly wider band.
mmse_rows <- dplyr::filter(gfap_forest, Covariate == "Baseline MMSE")
stopifnot(nrow(mmse_rows) == 2, max(abs(mmse_rows$Difference)) < 0.01)
other_rows <- dplyr::filter(gfap_forest, Covariate != "Baseline MMSE")
stopifnot(nrow(other_rows) == 4, max(abs(other_rows$Difference)) < 0.035)

All six rows land within 0.03 of the published panel. The MMSE rows agree to within 0.005 and the body-weight rows to within 0.03, which is what pins the two normalisation constants the paper never prints: 26 MMSE units and 72 kg.

Replicate Figure 5 – maintenance dosing

Figure 5 shows the model-predicted plasma biomarkers for continuous LEC10-BW over 48 months compared with a switch to monthly LEC10-M at 18, 24 or 30 months.

lbl <- c(Bhagunde_2026_lecanemab_abeta4240 = "(A) Plasma Abeta42/40 ratio",
         Bhagunde_2026_lecanemab_ptau181   = "(B) Plasma p-tau181 (pg/mL)",
         Bhagunde_2026_lecanemab_gfap      = "(C) Plasma GFAP (pg/mL)")

typ |>
  dplyr::filter(arm != "Discontinue at month 18") |>
  dplyr::mutate(panel = lbl[model], month = time / HR_MO) |>
  ggplot(aes(month, value, colour = arm)) +
  geom_line(linewidth = 0.7) +
  facet_wrap(~panel, ncol = 1, scales = "free_y") +
  scale_x_continuous(breaks = seq(0, 48, by = 6)) +
  labs(x = "Time (months)", y = NULL, colour = NULL,
       title = "Figure 5 -- transition to monthly maintenance dosing",
       caption = "Replicates Figure 5 of Bhagunde 2026 (typical-value predictions).") +
  theme(legend.position = "bottom")

typ_keep <- dplyr::filter(typ, arm != "Discontinue at month 18")
final48 <- typ_keep |>
  dplyr::distinct(model, arm) |>
  dplyr::mutate(month48 = mapply(function(m, a) {
    at_month(typ_keep[typ_keep$model == m & typ_keep$arm == a, ], "value", 48)
  }, model, arm)) |>
  dplyr::group_by(model) |>
  dplyr::mutate(ratio_to_bw = month48 / month48[arm == "LEC10-BW (continuous)"]) |>
  dplyr::ungroup()

final48 |>
  dplyr::mutate(Biomarker = lbl[model]) |>
  dplyr::select(Biomarker, arm, month48, ratio_to_bw) |>
  dplyr::rename("Arm" = arm, "Month 48 value" = month48,
                "Ratio to LEC10-BW" = ratio_to_bw) |>
  knitr::kable(digits = 4, caption = "Month-48 biomarker levels by maintenance-dose arm.")
Month-48 biomarker levels by maintenance-dose arm.
Biomarker Arm Month 48 value Ratio to LEC10-BW
(A) Plasma Abeta42/40 ratio LEC10-BW (continuous) 0.0949 1.0000
(B) Plasma p-tau181 (pg/mL) LEC10-BW (continuous) 2.3390 1.0000
(C) Plasma GFAP (pg/mL) LEC10-BW (continuous) 282.7263 1.0000
(A) Plasma Abeta42/40 ratio LEC10-M from month 18 0.0907 0.9565
(B) Plasma p-tau181 (pg/mL) LEC10-M from month 18 2.5267 1.0803
(C) Plasma GFAP (pg/mL) LEC10-M from month 18 286.6662 1.0139
(A) Plasma Abeta42/40 ratio LEC10-M from month 24 0.0908 0.9575
(B) Plasma p-tau181 (pg/mL) LEC10-M from month 24 2.5013 1.0694
(C) Plasma GFAP (pg/mL) LEC10-M from month 24 285.9018 1.0112
(A) Plasma Abeta42/40 ratio LEC10-M from month 30 0.0911 0.9605
(B) Plasma p-tau181 (pg/mL) LEC10-M from month 30 2.4714 1.0566
(C) Plasma GFAP (pg/mL) LEC10-M from month 30 285.2149 1.0088 

# The paper's two quantitative Figure 5 claims:
#  (1) Abeta42/40 "re-equilibrates to a LOWER plateau" under monthly dosing;
#  (2) p-tau181 and GFAP are "sustained at levels similar to IV LEC10-BW",
#      and results "did not vary based on when the transition occurred".
ab48 <- dplyr::filter(final48, model == "Bhagunde_2026_lecanemab_abeta4240")
stopifnot(all(ab48$ratio_to_bw[ab48$arm != "LEC10-BW (continuous)"] < 0.995))

pd48 <- dplyr::filter(final48, model != "Bhagunde_2026_lecanemab_abeta4240")
stopifnot(all(abs(pd48$ratio_to_bw - 1) < 0.10))

switch_spread <- final48 |>
  dplyr::filter(grepl("^LEC10-M", arm)) |>
  dplyr::group_by(model) |>
  dplyr::summarise(spread = (max(month48) - min(month48)) / mean(month48), .groups = "drop")
stopifnot(all(switch_spread$spread < 0.03))

The three arms behave as the paper describes, and the model puts numbers on its qualitative statements. Halving the dosing frequency at month 18 leaves the Abeta42/40 ratio about 4% below the continuous-biweekly plateau at month 48 (the paper’s “re-equilibrating to a lower plateau level”), p-tau181 about 8% above it and GFAP about 1% above it (the paper’s “sustained … at levels similar to IV LEC10-BW”). The ordering follows directly from how saturated each drug-effect function is at clinical exposure: the GFAP effect runs through the plaque model, whose 1 + DESLP * Conc term is far into its linear-in-log regime, so halving exposure barely moves it; the p-tau181 Emax function is only about three-quarters saturated, so it loses the most; the Abeta42/40 effect is exactly linear in concentration and loses in proportion to the exposure it gives up. Across the three switch times the month-48 value varies by under 3% – the paper’s “results did not vary based on when the transition occurred”.

Between-subject variability

A stochastic cohort on LEC10-BW illustrates the interindividual variability each model carries. Note the extremely wide EC50 variability the p-tau181 model reports (252 %CV), which dominates the spread of that panel.

set.seed(20260902)
n_vpc <- 150
vpc_cov <- tibble::tibble(
  id = seq_len(n_vpc),
  WT = round(stats::rlnorm(n_vpc, log(72), 0.20), 1),
  ALB = round(stats::rnorm(n_vpc, 43, 3), 1),
  SEXF = stats::rbinom(n_vpc, 1, 0.52),
  ADA_POS = 0,
  RACE_JAPANESE = stats::rbinom(n_vpc, 1, 0.06),
  AGE = round(stats::rnorm(n_vpc, 72, 8)),
  APOE4_CARRIER = stats::rbinom(n_vpc, 1, 0.69),
  SCORE_MMSE = pmin(30, pmax(20, round(stats::rnorm(n_vpc, 25.5, 2)))),
  STUDY_LEC201 = 0
)

vpc_events <- dplyr::bind_rows(lapply(seq_len(n_vpc), function(i) {
  make_events(id = vpc_cov$id[i],
              covs = as.list(vpc_cov[i, setdiff(names(vpc_cov), "id")]),
              horizon_mo = 48, obs_step_mo = 3)
}))
stopifnot(!anyDuplicated(unique(vpc_events[, c("id", "time", "evid")])))
vpc <- dplyr::bind_rows(lapply(model_names, function(mn) {
  s <- rxode2::rxSolve(readModelDb(mn), vpc_events, returnType = "data.frame")
  tibble::tibble(model = mn, id = s$id, time = s$time, value = s[[state_of[[mn]]]])
}))
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ parameter labels from comments will be replaced by 'label()'

vpc |>
  dplyr::filter(!is.na(value)) |>
  dplyr::group_by(model, time) |>
  dplyr::summarise(Q05 = stats::quantile(value, 0.05),
                   Q50 = stats::median(value),
                   Q95 = stats::quantile(value, 0.95), .groups = "drop") |>
  dplyr::mutate(panel = lbl[model], month = time / HR_MO) |>
  ggplot(aes(month, Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.25) +
  geom_line(linewidth = 0.7) +
  facet_wrap(~panel, ncol = 1, scales = "free_y") +
  scale_x_continuous(breaks = seq(0, 48, by = 6)) +
  labs(x = "Time (months)", y = NULL,
       title = "Simulated LEC10-BW cohort: median and 5th-95th percentiles",
       caption = paste("n =", n_vpc, "virtual subjects; covariate distributions are illustrative."))

# Direction-of-effect check on the simulated cohort medians: over 18 months the
# Abeta42/40 ratio must rise and the two damage markers must fall, which is the
# paper's central qualitative finding (Section 3.1).
dir_chk <- vpc |>
  dplyr::filter(!is.na(value)) |>
  dplyr::group_by(model, id) |>
  dplyr::summarise(rel18 = value[which.min(abs(time - 18 * HR_MO))] / value[which.min(time)],
                   .groups = "drop") |>
  dplyr::group_by(model) |>
  dplyr::summarise(`Median month-18 / baseline` = stats::median(rel18),
                   `10th pct` = stats::quantile(rel18, 0.1),
                   `90th pct` = stats::quantile(rel18, 0.9), .groups = "drop") |>
  dplyr::mutate(Biomarker = lbl[model]) |>
  dplyr::select(Biomarker, `Median month-18 / baseline`, `10th pct`, `90th pct`)

knitr::kable(dir_chk, digits = 3,
             caption = "Relative month-18 change on LEC10-BW across the simulated cohort.")
Relative month-18 change on LEC10-BW across the simulated cohort.
Biomarker Median month-18 / baseline 10th pct 90th pct
(A) Plasma Abeta42/40 ratio 1.105 1.036 1.271
(C) Plasma GFAP (pg/mL) 0.851 0.634 0.986
(B) Plasma p-tau181 (pg/mL) 0.822 0.691 0.923 

med <- dplyr::pull(dir_chk, `Median month-18 / baseline`)
names(med) <- dir_chk$Biomarker
stopifnot(
  med[["(A) Plasma Abeta42/40 ratio"]] > 1.05,
  med[["(B) Plasma p-tau181 (pg/mL)"]] < 0.95,
  med[["(C) Plasma GFAP (pg/mL)"]] < 0.95
)

Assumptions and deviations

Sources that had to be acquired

  • The lecanemab population PK model is not in this paper. Bhagunde 2026 Section 2.2.1 says “Final parameter estimates can be found in Table S1”, but Table S1 of the published supplement is Summary of Plasma Biomarker Observations and Subjects by Dosing Treatment – a data-count table with no PK parameters. The cross-reference is broken. The PK layer here is therefore taken from Majid 2024 (doi:10.1002/psp4.13224), which Bhagunde 2026 itself cites (reference 18) as the published source of the model, and every PK value is wrapped in fixed() to mark that this paper did not estimate it.
  • The paper’s PK model is a modification of Majid 2024. Section 2.2.1 states the previously developed model was “modified to include a more physiologically plausible relationship between lecanemab PK and manufacturing process”, moving the process effect from relative bioavailability onto V1 and V2. The modified estimates were not published. This extraction uses Majid 2024 as printed and omits the manufacturing-process comparability factor entirely (equivalent to F = 1, the Process A reference), because the F-based parameterisation is explicitly not the one this paper used. With F = 0.904 (Process B) instead, the month-18 amyloid plaque decrease for an APOE4 carrier would be 64.8 CL rather than 67.1 CL; both bracket the published 59.9 CL cohort mean, so this choice does not drive any conclusion above.
  • The amyloid-plaque model is not in this paper either. The plaque ODE is printed in Section 2.2.2 but its parameters are not tabulated. They come from the companion Bhagunde 2026 CPT:PSP paper (doi:10.1002/psp4.70173), Table 1, and are likewise fixed().

Centring / normalisation constants the paper does not print

Every covariate in Bhagunde 2026 Table 1 is reported without its normalisation value. Three constants were recovered, and each recovery is itself asserted above:

  • Body weight, 72 kg. Printed in the Majid 2024 PK equations, and independently recovered from Figure 4C: the plotted GFAP fold changes (50 kg -> 1.27, 100 kg -> 0.805) imply reference weights of 71.8 and 72.0 kg under the tabulated power of -0.66.
  • Age, 72 years. Section 3.2.1 states “Compared to a reference subject aged 72 years”. Centring the Abeta42/40 age-on-slope term at 72 years reproduces the Figure 4A age panel (83 years: model 1.24 vs plotted 1.21; 57 years: model 0.75 vs plotted 0.77). The same 72 years is used for the plaque model’s age-on-DESLP power, whose normalisation that paper also does not print.
  • Baseline MMSE, 26 units. Recovered from Figure 4C: 22 units -> 1.11 and 29 units -> 0.929 under the tabulated GFAP MMSE power of -0.63 imply 25.96 and 25.80 respectively. The same value is used for the p-tau181 MMSE power, whose reference the paper also omits. This matches the centring already used for SCORE_MMSE in Delor_2013_alzheimer.R.

Errata and internal inconsistencies in the published tables

  • EC50 units. Table 1 labels p-tau181 EC50 = 31.4 as “mg/mL”. Taken literally the drug effect would be about 0.2% at clinical exposures and p-tau181 would not move at all; read as ug/mL the model reproduces the paper’s own 3.37 -> 2.80 pg/mL trajectory to within 1% (see Validation 2). The model encodes ug/mL.
  • GFAP KinG versus the stated GFAP half-life. Table 1 reports KinG = 5.6 pg/mL/h, and the text defines KoutG = KinG / BGFAP, giving KoutG = 0.0178 /h (half-life 39 h). Section 3.2.3 separately states “The typical t1/2 of plasma GFAP to return to baseline following lecanemab treatment discontinuation was estimated to be 1.7 years”, which would require KoutG = 0.408 /year and hence KinG = 128 pg/mL/year, a value 23-fold different in any consistent time unit. The two statements cannot both hold. The tabulated KinG = 5.6 pg/mL/h is encoded, because it is the parameter table and because it reproduces the paper’s own quantitative GFAP trajectory (316 -> 275 pg/mL at 18 months; the model gives 315 -> 269, within 2.2%), whereas KoutG = 0.408 /year would require a relative plaque reduction above 100% to reach 275 pg/mL and is therefore arithmetically excluded. The consequence is that the packaged model returns toward baseline faster than 1.7 years after discontinuation.
  • GFAP proportional residual error. Table 1 reports 0.193 under a “%CV” heading; taken literally that is 0.193%, which is implausible for a Simoa biomarker assay. The table footnote defines proportional error %CV as “square root of variance x 100”, so 0.193 is the standard deviation, i.e. 19.3%. propSd = 0.193 is encoded.
  • Figure 4A is not self-consistent about the PK covariate model. The two covariates in the Abeta42/40 forest panel that also act on lecanemab PK – sex and body weight – are the only two the packaged model fails to reproduce, and they fail in opposite directions.
    • Age, Japanese race and APOE4 do reproduce: 83 years gives 1.212 against a plotted 1.21, 57 years gives 0.769 against 0.77, Japanese 1.037 against 1.03, APOE4 carrier 1.014 against 0.99. This confirms the tabulated baseline and slope coefficients and the 72-year centring.
    • Sex: the packaged model gives 1.25 against a plotted 1.01. A female’s clearance is 0.791 times a male’s (Majid 2024), which raises steady-state exposure – and hence the Abeta42/40 change from baseline – by about 26%. The panel shows no such effect, i.e. it behaves as though the PK covariate model were switched off.
    • Body weight: the packaged model gives 1.044 (97.5 kg) and 0.890 (48.2 kg) against plotted 1.18 and 0.78. The plotted values are instead reproduced to within 0.01 by the PK-exposure effect of weight alone: a 10 mg/kg dose with CL proportional to WT^0.353 makes steady-state exposure proportional to WT^0.647, giving 1.22 and 0.77. That is, the weight row behaves as though only the PK covariate model were on and the tabulated Weight on slope (exponential) = -0.00602 had not been applied.
    No single reading of the panel reproduces both rows. Table 1 is the parameter table and outranks a figure, so the tabulated coefficients are encoded as printed – exp(theta * (COV - reference)), the form the neighbouring Age on slope (exponential) = 0.0193 row unambiguously takes and which that row confirms numerically. Both discrepancies are reported in the Figure 4A comparison above and pinned by assertions.
  • Figure 4A 48.2 kg annotation. The panel prints “0.7 [0.73, 0.84]”, a point estimate outside its own confidence interval. The interval midpoint (0.78) is used as the reference value above.

Modelling choices

  • Time unit. All three models run in hours so the inherited PK parameters (L/h) are used exactly as published. Kout, DP and the plaque Kout, reported per year, are divided by 8766 h/year inside model().
  • Infusion duration. Neither Bhagunde 2026 nor Majid 2024 states the zero-order input duration. A 1-hour infusion is used in this vignette. It has no material effect on biomarkers whose turnover half-lives are months to years.
  • Cc carries no residual error. Bhagunde 2026 drove the PD models with post hoc individual PK predictions, so Cc is encoded as an algebraic individual prediction and no PK residual-error model is attached. The lecanemab PK residual error reported by Majid 2024 (proportional 21.0 %CV, additive SD 1.12 ug/mL) is therefore deliberately not carried.
  • Plaque layer variability. The upstream amyloid-PET model estimates IIV on baseline plaque, on an additive “amyloid floor” term and on DESLP. Because Bhagunde 2026 used the plaque model only as a predictor of GFAP production, the plaque layer is encoded as a typical-value trajectory with no IIV, and the amyloid floor (typical value fixed to 0 in the upstream paper) is omitted.
  • Paper-specific output states. The three observation variables (abeta4240, ptau181, gfap) and the plaque state are declared through paper_specific_compartments, following vanMaanen_2025_amyloid.R. checkModelConventions() raises one warning per model that the single-output observation variable is not a registered canonical; that is the same accepted warning vanMaanen_2025_amyloid.R carries for plaque.
  • Virtual cohort covariates are illustrative. Bhagunde 2026 does not tabulate baseline demographics (Appendix A refers to the study-design publications). The distributions used for the stochastic cohort figure are chosen to be plausible for an early-Alzheimer’s trial population and are not taken from the paper; none of the numeric assertions in this vignette depends on them, which all use the reference subject.
  • ADA status. ADA_POS is held at 0 throughout. Lecanemab ADA incidence is not reported in this paper, and the effect on CL is a 13% increase, well inside the covariate bands examined above.