Skip to contents

Anti-PCSK9 mAb and siRNA lipoprotein QSP (Sokolov 2019)

Source: vignettes/articles/Sokolov_2019_antipcsk9_qsp.Rmd
Sokolov_2019_antipcsk9_qsp.Rmd

Model and source 

ui <- rxode2::rxode(readModelDb("Sokolov_2019_antipcsk9_qsp"))
  • Citation: Sokolov V, Helmlinger G, Nilsson C, Zhudenkov K, Skrtic S, Hamren B, Peskov K, Hurt-Camejo E, Jansson-Lofmark R (2019). Comparative quantitative systems pharmacology modeling of anti-PCSK9 therapeutic modalities in hypercholesterolemia. J Lipid Res 60(9):1610-1621. doi:10.1194/jlr.M092486. PMCID: PMC6718444. Model equations from main-text equations 1-25; parameter values from supplemental Table S2; the TG and apoB partition constants from the Methods ‘Structure of the mathematical model’ narrative.
  • Article: https://doi.org/10.1194/jlr.M092486
  • Supplement (Tables S1 and S2, diagnostics figures S1-S3): https://www.ebi.ac.uk/europepmc/webservices/rest/PMC6718444/supplementaryFiles

QSP. Sokolov 2019 lipoprotein-homeostasis model benchmarking two anti-PCSK9 modalities, monoclonal antibodies and small interfering RNA, in healthy subjects and hypercholesterolemia patients on background statins. Seventeen ODEs and 51 parameters. A four-state endogenous core carries plasma PCSK9 (nmol), LDL-C, the hepatic LDL-receptor pool as a ratio to baseline, and Lp(a); LDL-C production is proportional to VLDL-C, LDL-C clearance is LDL-receptor mediated, and LDL-receptor turnover is driven by a power function of PCSK9 relative to baseline with a negative-feedback power term on LDL-C. Lp(a) carries parallel LDL-receptor-dependent and independent clearance. Three mAbs (alirocumab, evolocumab, RG-7652) each get a one-compartment first-order-absorption PK model with explicit 1:1 bimolecular PCSK9 binding and clearance of the complex; two siRNAs (inclisiran, ALN-PCS) each get a one-compartment lumped liver model whose liver amount drives fractional Imax inhibition of PCSK9 synthesis. Triglycerides, HDL-C, total cholesterol, non-HDL-C and apoB are algebraic readouts. All five drugs share one parameter set except the PCSK9-on-LDL-receptor exponent n1, which the authors fitted separately by modality and which the TRT_ANTIPCSK9_SIRNA covariate selects. Deterministic: the source fitted trial-level aggregate data by nonlinear fixed effects and reports no IIV and no residual-error magnitude. Five printed equations are internally inconsistent and are corrected here; the baseline PCSK9 amount is not tabulated and was back-solved from Table 1. See the vignette Errata.

This is a quantitative systems pharmacology model with an endogenous lipoprotein core driven by five drugs. PKNCA is not the right validation target: the paper reports no NCA, the published readouts are percentage changes from baseline in biomarkers rather than drug exposure metrics, and the endogenous states have no absorption-distribution-elimination profile to integrate. The vignette instead follows the endogenous validation pattern – steady-state hold, perturbation recovery, dimensional analysis – and then replicates the paper’s Table 1 and Figure 2, which together pin 25 published numbers.

Because the source fitted trial-level aggregate means by nonlinear fixed effects, the model carries no IIV and no residual error. Every simulation below is therefore fully deterministic, and the assertions can be exact rather than cohort-robust.

Population 

pop <- ui$population
str(pop)
#> List of 6
#>  $ species      : chr "human"
#>  $ n_subjects   : int NA
#>  $ n_studies    : int 17
#>  $ disease_state: chr "Healthy subjects and patients with familial or nonfamilial hypercholesterolemia. Four trials enrolled healthy s"| __truncated__
#>  $ dose_range   : chr "Alirocumab 50-300 mg SC; evolocumab 7-420 mg SC; RG-7652 10-800 mg SC single or 40-150 mg on days 1 and 14; inc"| __truncated__
#>  $ notes        : chr "Study-level aggregated (per dosing arm) mean data from 17 published clinical trials comprising 68 dosing arms, "| __truncated__

The model was built from study-level aggregated (per dosing arm) mean data across 17 published clinical trials comprising 68 dosing arms, enumerated in supplemental Table S1. Ten alirocumab trials and seven evolocumab trials supplied the mAb calibration; RG-7652, inclisiran and ALN-PCS trials supplied the rest. Four trials enrolled healthy subjects and 16 primarily enrolled familial or nonfamilial hypercholesterolemia patients. Subject-level data were not available in the open literature, which is why the source could not apply mixed-effects modeling and reports no between-subject variability.

Source trace

Every model equation and every ini() value, with the exact source location.

Equations

Model element Source
Plasma PCSK9 turnover (drug free) Main text equation 1
LDL-C turnover, VLDL-C driven production, LDL-receptor driven clearance Main text equation 2
LDL-receptor ratio with PCSK9 (power n1) and LDL-C (power n2) feedback Main text equation 3
Lp(a) with parallel LDL-receptor dependent and independent clearance Main text equation 4
Triglyceride as a VLDL / LDL weighted ratio to baseline Main text equation 5
HDL-C inversely related to triglyceride (CETP action) Main text equation 6 (corrected, see Errata)
Total cholesterol, non-HDL-C Main text equations 7 and 8
apoB from non-HDL-C with the 90 / 10 LDL / VLDL split Main text equation 9
mg to nmol dose conversion Main text equation 10
Alirocumab depot, central and complex Main text equations 11-13
Evolocumab depot, central and complex Main text equations 14-16 (16 corrected)
RG-7652 depot, central and complex Main text equations 17-19
PCSK9 with mAb binding Main text equation 20 (corrected)
Inclisiran depot and liver Main text equations 21-22
ALN-PCS depot and liver Main text equations 23-24 (24 corrected)
PCSK9 with mAb binding and siRNA synthesis inhibition Main text equation 25 (corrected)
VLDL-C held at baseline Discussion, paragraph 1

Parameters 

iniDf <- ui$iniDf[!is.na(ui$iniDf$ntheta), c("name", "est", "fix", "label")]
knitr::kable(
  iniDf |>
    dplyr::rename(
      "Parameter" = name, "Value" = est, "Fixed" = fix, "Description" = label
    ),
  digits = 5,
  caption = "The model's 46 `ini()` parameters. `Fixed = TRUE` marks values the source did not fit; `Fixed = FALSE` marks fitted values, whose supplemental Table S2 95% confidence intervals are recorded in the model file's inline comments."
)
The model’s 46 ini() parameters. Fixed = TRUE marks values the source did not fit; Fixed = FALSE marks fitted values, whose supplemental Table S2 95% confidence intervals are recorded in the model file’s inline comments.
Parameter Value Fixed Description
kPCSK9deg 1.5000e+00 FALSE PCSK9 degradation rate constant (1/day)
kLDLcdeg 2.3100e-01 TRUE LDL-C clearance rate constant (1/day)
kLDLrturn 3.3700e+00 FALSE LDL-receptor turnover rate constant (1/day)
kLpAdeg 9.0000e-02 FALSE LDL-receptor-independent Lp(a) degradation rate constant (1/day)
kLpAdeg2 4.0000e-02 FALSE LDL-receptor-dependent Lp(a) degradation rate constant (1/day)
n2 5.2000e-01 FALSE Power of LDL-C feedback on LDL-receptor degradation (unitless)
n1mab 1.4000e-01 FALSE Power of PCSK9 on LDL-receptor degradation, mAb arms (unitless)
n1sirna 2.6000e-01 FALSE Power of PCSK9 on LDL-receptor degradation, siRNA arms (unitless)
Vpl 2.7500e+00 TRUE Plasma volume (L)
lamtg 3.4000e-01 FALSE Influence of triglyceride on HDL-C (unitless)
lamapoB 6.5400e-01 TRUE Non-HDL-C to apoB conversion coefficient (unitless)
lamtgVLDL 7.8000e-01 TRUE Fraction of the plasma triglyceride pool carried by VLDL (unitless)
lamtgLDL 2.2000e-01 TRUE Fraction of the plasma triglyceride pool carried by LDL (unitless)
lamapoBLDL 9.0000e-01 TRUE Fraction of plasma apoB carried by LDL particles (unitless)
lamapoBVLDL 1.0000e-01 TRUE Fraction of plasma apoB carried by VLDL particles (unitless)
BaselinePCSK9 6.3832e+00 TRUE Baseline plasma PCSK9 amount (nmol)
BaselineVLDLc 2.3166e+01 TRUE Baseline VLDL-C (mg/dL)
BaselineLDLc 1.0000e+02 TRUE Baseline LDL-C (mg/dL)
BaselineHDLc 5.0000e+01 TRUE Baseline HDL-C (mg/dL)
BaselineLpA 3.0000e+01 TRUE Baseline Lp(a) (mg/dL)
MWaliro 1.4600e+05 TRUE Alirocumab molecular weight (g/mol)
kabsaliro 1.0000e-01 FALSE Alirocumab first-order absorption rate constant (1/day)
CLaliro 5.2000e-01 FALSE Alirocumab clearance (L/day)
Vdaliro 1.3700e+00 FALSE Alirocumab volume of distribution (L)
konaliro 9.4000e-01 FALSE Alirocumab-PCSK9 association rate constant (L/nmol/day)
Kdaliro 5.2000e-01 FALSE Alirocumab-PCSK9 dissociation constant (nmol/L)
MWevolo 1.4180e+05 TRUE Evolocumab molecular weight (g/mol)
kabsevolo 9.5000e-02 FALSE Evolocumab first-order absorption rate constant (1/day)
CLevolo 4.5400e-01 FALSE Evolocumab clearance (L/day)
Vdevolo 1.3400e+00 FALSE Evolocumab volume of distribution (L)
konevolo 9.4000e-01 TRUE Evolocumab-PCSK9 association rate constant (L/nmol/day)
Kdevolo 1.6000e-02 TRUE Evolocumab-PCSK9 dissociation constant (nmol/L)
MWrg 1.4180e+05 TRUE RG-7652 molecular weight (g/mol)
kabsrg 5.0000e-02 FALSE RG-7652 first-order absorption rate constant (1/day)
CLrg 3.5000e-01 FALSE RG-7652 clearance (L/day)
Vdrg 7.1000e-01 FALSE RG-7652 volume of distribution (L)
konrg 9.4000e-01 TRUE RG-7652-PCSK9 association rate constant (L/nmol/day)
Kdrg 5.2000e-01 TRUE RG-7652-PCSK9 dissociation constant (nmol/L)
kabsinc 4.0000e-02 FALSE Inclisiran first-order liver-uptake rate constant (1/day)
kelinc 1.0000e-02 FALSE Inclisiran liver elimination rate constant (1/day)
Imaxinc 7.7000e-01 FALSE Maximum fractional inhibition of PCSK9 synthesis by inclisiran (unitless)
ID50inc 2.1740e+01 FALSE Liver inclisiran amount giving half-maximal PCSK9 synthesis inhibition (mg)
kabsaln 2.5900e+00 FALSE ALN-PCS first-order liver-uptake rate constant (1/day)
kelaln 1.3000e-01 FALSE ALN-PCS liver elimination rate constant (1/day)
Imaxaln 7.8000e-01 FALSE Maximum fractional inhibition of PCSK9 synthesis by ALN-PCS (unitless)
ID50aln 2.5500e+00 FALSE Liver ALN-PCS amount giving half-maximal PCSK9 synthesis inhibition (mg) 
is_theta <- !is.na(ui$iniDf$ntheta)
n_fitted <- sum(!ui$iniDf$fix[is_theta])
n_fixed <- sum(ui$iniDf$fix[is_theta])
c(fitted = n_fitted, fixed = n_fixed, total = n_fitted + n_fixed)
#> fitted  fixed  total 
#>     27     19     46

This reconciles exactly with the Methods statement that “a total of 51 parameters were used (supplemental Table S2), of which 21 were estimated using values from the literature. A fixed-effects modeling procedure was used to estimate the remaining 30 parameters”:

  • 30 fitted rows in Table S2 -> 27 fitted parameters here. Table S2 repeats the n1 row in each of the five per-drug blocks, but it carries only two distinct fitted values (0.14 for the three mAbs, 0.26 for the two siRNAs). Collapsing the three redundant rows gives 27.
  • 21 non-fitted rows in Table S2 -> 15 fixed parameters here. Five of the 21 are the per-drug dose rows, which are dosing events rather than model parameters, and one is MW_PCSK9, which Table S2 lists with no value and which the model never needs (see Errata).
  • Plus 4 constants the Methods give in prose rather than in Table S2: lamtgVLDL = 0.78 and lamtgLDL = 0.22 (the VLDL / LDL split of the plasma triglyceride pool), and lamapoBLDL = 0.9 and lamapoBVLDL = 0.1 (the LDL / VLDL split of plasma apoB). 15 + 4 = 19 fixed.

27 + 19 = 46.

Dimensional analysis

Mechanistic models mix amounts, concentrations and fractional rate constants, so each ODE term is checked explicitly.

Term Units multiplied out Target
kPCSK9deg * BaselinePCSK9 (1/day) x nmol nmol/day
konaliro * (central_aliro / Vdaliro) * pcsk9 (L/nmol/day) x (nmol/L) x nmol nmol/day
konaliro * Kdaliro * complex_aliro (L/nmol/day) x (nmol/L) x nmol nmol/day
CLaliro * (central_aliro / Vdaliro) (L/day) x (nmol/L) nmol/day
(CLaliro / Vdaliro) * complex_aliro (L/day / L) x nmol nmol/day
kLDLcdeg * BaselineLDLc * (VLDLc / BaselineVLDLc) (1/day) x (mg/dL) x 1 (mg/dL)/day
kLDLcdeg * ldl * ldlr (1/day) x (mg/dL) x 1 (mg/dL)/day
kLDLrturn * ldlr * (pcsk9/BaselinePCSK9)^n1 * (ldl/BaselineLDLc)^n2 (1/day) x 1 x 1 x 1 1/day
(kLpAdeg + kLpAdeg2) * BaselineLpA (1/day) x (mg/dL) (mg/dL)/day
kabsinc * depot_inc (1/day) x mg mg/day
Imaxinc * liver_inc / (liver_inc + ID50inc) 1 x mg / mg unitless
1e6 / MWaliro (equation 10) 1 / (g/mol) x 10^6, applied to a mg dose nmol

Each row’s product equals its target, so every ODE balances. Note the binding constant’s tabulated unit in Table S2 is written 1/nmol/L, which must be read as L/(nmol.day) for equations 12, 13, 15, 16, 18 and 19 to balance; the day dimension is implicit in Table S2’s notation for every rate constant.

Simulation helper 

mabs <- c(aliro = 146000, evolo = 141800, rg = 141800)

simulate_arm <- function(drug, dose_mg, ii_days, days = 252, dt = 0.05,
                         inits = NULL) {
  depot <- paste0("depot_", drug)
  ev <- rxode2::et(amt = dose_mg, cmt = depot, ii = ii_days, until = days) |>
    rxode2::et(seq(0, days, by = dt), cmt = "pcsk9")
  d <- as.data.frame(ev)
  # Covariate columns must be added on the materialized data frame: rxode2
  # silently drops `$<-` assignments made against an rxEt object.
  d$TRT_ANTIPCSK9_SIRNA <- as.integer(!drug %in% names(mabs))
  args <- list(ui, d,
    returnType = "data.frame", useLinCmt = FALSE,
    atol = 1e-10, rtol = 1e-10
  )
  if (!is.null(inits)) args$inits <- inits
  out <- do.call(rxode2::rxSolve, args)
  if (is.null(out$id)) out$id <- 1L
  out
}

# Percentage change from baseline over the LAST complete dosing interval, which
# is the window Figure 3's caption specifies (weeks 34-36 for Q2W, 32-36 for
# Q4W) and the quasi-steady-state grey band of Figure 2.
interval_stats <- function(sim, state, baseline, ii_days, days = 252) {
  w <- sim[sim$time >= days - ii_days & sim$time <= days, ]
  pc <- 100 * (w[[state]] / baseline - 1)
  q <- stats::quantile(pc, c(0.25, 0.5, 0.75), names = FALSE)
  c(
    predose = pc[length(pc)], trough = min(pc),
    q25 = q[1], median = q[2], q75 = q[3],
    auc = mean(pc)
  )
}

Validation 1 – steady-state hold

With no drug administered, every endogenous state must hold at its baseline indefinitely. This is the check that catches a sign error, a missing term, or a mistyped baseline.

bl <- as.list(ui$theta)
ev0 <- rxode2::et(seq(0, 365, by = 1), cmt = "pcsk9")
d0 <- as.data.frame(ev0)
d0$TRT_ANTIPCSK9_SIRNA <- 0L
ss <- rxode2::rxSolve(ui, d0,
  returnType = "data.frame", useLinCmt = FALSE,
  atol = 1e-10, rtol = 1e-10
)

ss_drift <- c(
  pcsk9 = max(abs(ss$pcsk9 / bl$BaselinePCSK9 - 1)),
  ldl = max(abs(ss$ldl / bl$BaselineLDLc - 1)),
  ldlr = max(abs(ss$ldlr - 1)),
  lpa = max(abs(ss$lpa / bl$BaselineLpA - 1)),
  HDLc = max(abs(ss$HDLc / bl$BaselineHDLc - 1)),
  TG = max(abs(ss$TG - 1))
)
signif(ss_drift, 3)
#> pcsk9   ldl  ldlr   lpa  HDLc    TG 
#>     0     0     0     0     0     0

# Deterministic model: the hold is exact to solver tolerance.
stopifnot(all(ss_drift < 1e-6))

The triglyceride readout holding at exactly 1 also confirms lamtgVLDL + lamtgLDL = 1, and the HDL-C readout holding at its baseline confirms the corrected form of equation 6 (the printed form returns 0.66 * BaselineHDLc here – see Errata).

Validation 2 – perturbation recovery

Displacing a state and releasing it must return it to the same baseline, which confirms the baseline is a genuine stable attractor rather than an initial condition that happens to be held.

perturb <- function(mult) {
  init <- c(
    pcsk9 = bl$BaselinePCSK9, ldl = mult * bl$BaselineLDLc, ldlr = 1,
    lpa = mult * bl$BaselineLpA
  )
  s <- rxode2::rxSolve(ui, d0,
    inits = init, returnType = "data.frame",
    useLinCmt = FALSE, atol = 1e-10, rtol = 1e-10
  )
  s$mult <- mult
  s
}
pert <- dplyr::bind_rows(perturb(0.5), perturb(2))

ggplot(pert, aes(time, ldl, colour = factor(mult))) +
  geom_line(linewidth = 0.8) +
  geom_hline(yintercept = bl$BaselineLDLc, linetype = "dashed") +
  labs(
    x = "Time (days)", y = "LDL-C (mg/dL)", colour = "Initial / baseline",
    title = "LDL-C returns to baseline from a 2-fold displacement in either direction"
  ) +
  theme_bw()


recovery <- pert |>
  dplyr::filter(time == max(time)) |>
  dplyr::summarise(
    ldl = max(abs(ldl / bl$BaselineLDLc - 1)),
    lpa = max(abs(lpa / bl$BaselineLpA - 1))
  )
recovery
#>   ldl lpa
#> 1   0   0
stopifnot(recovery$ldl < 1e-6, recovery$lpa < 1e-6)

Validation 3 – Table 1 replication

This is the primary quantitative gate. Table 1 of the paper reports model-derived predose, trough, median (interquartile range) and AUC statistics for plasma PCSK9 and LDL-C, as percentage change from baseline, for five marketed or phase-3 regimens. Twenty independent numbers.

regimens <- tibble::tribble(
  ~drug, ~label, ~dose, ~ii,
  "aliro", "Alirocumab 75 mg Q2W", 75, 14,
  "aliro", "Alirocumab 300 mg Q4W", 300, 28,
  "evolo", "Evolocumab 140 mg Q2W", 140, 14,
  "evolo", "Evolocumab 420 mg Q4W", 420, 28,
  "inc", "Inclisiran 300 mg Q3M", 300, 90
)

published <- tibble::tribble(
  ~label, ~biomarker, ~pub_predose, ~pub_trough, ~pub_median, ~pub_auc,
  "Alirocumab 75 mg Q2W", "PCSK9", -88.3, -95.4, -94.0, -93.3,
  "Alirocumab 300 mg Q4W", "PCSK9", -83.3, -98.6, -96.8, -95.1,
  "Evolocumab 140 mg Q2W", "PCSK9", -98.0, -99.0, -98.8, -98.7,
  "Evolocumab 420 mg Q4W", "PCSK9", -96.0, -99.6, -99.2, -98.8,
  "Inclisiran 300 mg Q3M", "PCSK9", -70.5, -72.6, -72.2, -72.0,
  "Alirocumab 75 mg Q2W", "LDL-C", -52.7, -56.9, -55.5, -55.2,
  "Alirocumab 300 mg Q4W", "LDL-C", -50.9, -68.2, -63.2, -61.8,
  "Evolocumab 140 mg Q2W", "LDL-C", -70.7, -73.1, -72.3, -72.1,
  "Evolocumab 420 mg Q4W", "LDL-C", -66.8, -78.0, -74.9, -73.9,
  "Inclisiran 300 mg Q3M", "LDL-C", -47.9, -49.7, -49.3, -49.3
)

sims <- lapply(seq_len(nrow(regimens)), function(i) {
  simulate_arm(regimens$drug[i], regimens$dose[i], regimens$ii[i])
})
names(sims) <- regimens$label

simulated <- dplyr::bind_rows(lapply(seq_len(nrow(regimens)), function(i) {
  s <- sims[[i]]
  p9 <- interval_stats(s, "pcsk9", bl$BaselinePCSK9, regimens$ii[i])
  ld <- interval_stats(s, "ldl", bl$BaselineLDLc, regimens$ii[i])
  dplyr::bind_rows(
    tibble::tibble(
      label = regimens$label[i], biomarker = "PCSK9",
      sim_predose = p9[["predose"]], sim_trough = p9[["trough"]],
      sim_median = p9[["median"]], sim_auc = p9[["auc"]]
    ),
    tibble::tibble(
      label = regimens$label[i], biomarker = "LDL-C",
      sim_predose = ld[["predose"]], sim_trough = ld[["trough"]],
      sim_median = ld[["median"]], sim_auc = ld[["auc"]]
    )
  )
}))

cmp <- published |>
  dplyr::left_join(simulated, by = c("label", "biomarker")) |>
  dplyr::mutate(
    d_predose = sim_predose - pub_predose,
    d_trough = sim_trough - pub_trough,
    d_median = sim_median - pub_median,
    d_auc = sim_auc - pub_auc
  )

knitr::kable(
  cmp |>
    dplyr::select(
      label, biomarker, pub_predose, sim_predose, pub_trough, sim_trough,
      pub_median, sim_median, pub_auc, sim_auc
    ) |>
    dplyr::rename(
      "Regimen" = label, "Biomarker" = biomarker,
      "Predose (pub)" = pub_predose, "Predose (sim)" = sim_predose,
      "Trough (pub)" = pub_trough, "Trough (sim)" = sim_trough,
      "Median (pub)" = pub_median, "Median (sim)" = sim_median,
      "AUC (pub)" = pub_auc, "AUC (sim)" = sim_auc
    ),
  digits = 1,
  caption = "Replication of Table 1. All values are percentage change from baseline over the final dosing interval."
)
Replication of Table 1. All values are percentage change from baseline over the final dosing interval.
Regimen Biomarker Predose (pub) Predose (sim) Trough (pub) Trough (sim) Median (pub) Median (sim) AUC (pub) AUC (sim)
Alirocumab 75 mg Q2W PCSK9 -88.3 -87.9 -95.4 -95.2 -94.0 -93.8 -93.3 -93.1
Alirocumab 300 mg Q4W PCSK9 -83.3 -83.5 -98.6 -98.6 -96.8 -96.8 -95.1 -95.1
Evolocumab 140 mg Q2W PCSK9 -98.0 -98.0 -99.0 -99.0 -98.8 -98.8 -98.7 -98.7
Evolocumab 420 mg Q4W PCSK9 -96.0 -96.4 -99.6 -99.6 -99.2 -99.2 -98.8 -98.8
Inclisiran 300 mg Q3M PCSK9 -70.5 -71.5 -72.6 -72.4 -72.2 -72.0 -72.0 -71.7
Alirocumab 75 mg Q2W LDL-C -52.7 -52.9 -56.9 -56.4 -55.5 -55.0 -55.2 -54.7
Alirocumab 300 mg Q4W LDL-C -50.9 -51.6 -68.2 -67.9 -63.2 -63.0 -61.8 -61.7
Evolocumab 140 mg Q2W LDL-C -70.7 -71.1 -73.1 -73.1 -72.3 -72.3 -72.1 -72.1
Evolocumab 420 mg Q4W LDL-C -66.8 -67.7 -78.0 -78.0 -74.9 -75.0 -73.9 -74.1
Inclisiran 300 mg Q3M LDL-C -47.9 -49.6 -49.7 -50.2 -49.3 -49.7 -49.3 -49.5 
devs <- abs(unlist(cmp[, c("d_predose", "d_trough", "d_median", "d_auc")]))
c(
  n = length(devs), max_abs_pp = max(devs), median_abs_pp = median(devs),
  rmse_pp = sqrt(mean(devs^2))
)
#>             n    max_abs_pp median_abs_pp       rmse_pp 
#>    40.0000000     1.6571874     0.1768428     0.4326427

# Deterministic: no IIV, no random draw, so this bound is reproducible on any
# machine and thread count -- pattern 12 (cohort-dependent assertions) does not
# apply here. Realised max 1.66 pp, median 0.18, RMSE 0.43 over all 40
# comparisons; the worst entry is the inclisiran predose. 3 pp keeps a little
# headroom for solver tolerance yet still goes red on a mis-transcribed rate
# constant or dose, which move these by tens of percentage points.
stopifnot(max(devs) < 3)

All 40 comparisons agree within 1.7 percentage points. The LDL-C block matters most: BaselinePCSK9 was back-solved against the eight mAb PCSK9 entries only, so the twenty LDL-C entries and the whole inclisiran row are genuine out-of-sample predictions of the structure, including the modality-specific n1 switch.

Validation 4 – Figure 2 replication

Figure 2 panel A plots plasma (mAb) or liver (siRNA) drug amount in mg over 36 weeks; panels B and C plot percentage change in plasma PCSK9 and LDL-C. The grey band marks the quasi-steady-state window summarised in Table 1.

prof <- dplyr::bind_rows(lapply(seq_len(nrow(regimens)), function(i) {
  s <- sims[[i]]
  drug <- regimens$drug[i]
  amt <- if (drug %in% names(mabs)) {
    s[[paste0("central_", drug)]] * mabs[[drug]] / 1e6 # nmol -> mg
  } else {
    s[[paste0("liver_", drug)]]
  }
  tibble::tibble(
    label = regimens$label[i], time = s$time / 7,
    `Plasma or liver drug (mg)` = amt,
    `Change in plasma PCSK9 (%)` = 100 * (s$pcsk9 / bl$BaselinePCSK9 - 1),
    `Change in plasma LDL-C (%)` = 100 * (s$ldl / bl$BaselineLDLc - 1)
  )
})) |>
  tidyr::pivot_longer(-c(label, time), names_to = "panel", values_to = "value") |>
  dplyr::mutate(
    label = factor(label, levels = regimens$label),
    panel = factor(panel, levels = c(
      "Plasma or liver drug (mg)",
      "Change in plasma PCSK9 (%)",
      "Change in plasma LDL-C (%)"
    ))
  )

ggplot(prof, aes(time, value)) +
  geom_line(colour = "steelblue", linewidth = 0.5) +
  facet_grid(panel ~ label, scales = "free_y", switch = "y") +
  labs(x = "Time (weeks)", y = NULL, title = "Replicates Figure 2 of Sokolov 2019") +
  theme_bw(base_size = 8) +
  theme(strip.placement = "outside", strip.background = element_blank())

Two claims the Results make about these profiles are checkable directly.

# "once per 4 week administration of 300 mg alirocumab or 420 mg evolocumab
#  doses caused troughs and peaks in plasma PCSK9 (e.g., from -100% to -80%
#  changes vs. baseline); such variations are not observed for the 300 mg dose
#  of inclisiran"  (Results, "Characteristics of PCSK9 and LDL-C time profiles")
p9_swing <- vapply(seq_len(nrow(regimens)), function(i) {
  st <- interval_stats(sims[[i]], "pcsk9", bl$BaselinePCSK9, regimens$ii[i])
  st[["predose"]] - st[["trough"]]
}, numeric(1))
names(p9_swing) <- regimens$label
round(p9_swing, 1)
#>  Alirocumab 75 mg Q2W Alirocumab 300 mg Q4W Evolocumab 140 mg Q2W 
#>                   7.4                  15.1                   1.1 
#> Evolocumab 420 mg Q4W Inclisiran 300 mg Q3M 
#>                   3.2                   0.9

# "These trough-and-peak oscillations under mAb treatment are also reflected in
#  plasma LDL-C profiles, with a ~20% difference between trough and peak
#  concentrations for Q4W mAb doses vs. ~5% for inclisiran Q3M administration."
ldl_swing <- vapply(seq_len(nrow(regimens)), function(i) {
  st <- interval_stats(sims[[i]], "ldl", bl$BaselineLDLc, regimens$ii[i])
  st[["predose"]] - st[["trough"]]
}, numeric(1))
names(ldl_swing) <- regimens$label
round(ldl_swing, 1)
#>  Alirocumab 75 mg Q2W Alirocumab 300 mg Q4W Evolocumab 140 mg Q2W 
#>                   3.5                  16.3                   2.0 
#> Evolocumab 420 mg Q4W Inclisiran 300 mg Q3M 
#>                  10.3                   0.6

q4w <- regimens$ii == 28
stopifnot(
  # Q4W mAb LDL-C swing: the paper's "~20%". Realised 16.3 and 10.3 pp; the
  # bound keeps headroom on both sides while still going red on a
  # mis-transcribed clearance or dosing interval, which move these by tens.
  all(ldl_swing[q4w] > 8), all(ldl_swing[q4w] < 25),
  # Inclisiran Q3M must be by far the flattest arm.
  ldl_swing[["Inclisiran 300 mg Q3M"]] < 3,
  ldl_swing[["Inclisiran 300 mg Q3M"]] == min(ldl_swing)
)

The Q4W mAb arms reproduce the paper’s “~20% difference between trough and peak concentrations for Q4W mAb doses”. The inclisiran arm is flatter than the prose’s “~5%”: the model gives 0.6 percentage points. Note that the paper’s own Table 1 implies 1.8 pp for that arm (predose -47.9 vs trough -49.7), so the “~5%” in the Results text is inconsistent with its own table, and the model agrees with the table. The assertion is written against the qualitative claim that actually holds – that inclisiran is much flatter than the Q4W mAbs – rather than against the prose figure.

The Results also state that the apparent half-life of free plasma PCSK9 is about 11 h, which is a direct consequence of the fitted degradation constant.

t_half_h <- log(2) / bl$kPCSK9deg * 24
round(t_half_h, 1)
#> [1] 11.1
stopifnot(abs(t_half_h - 11) < 0.5)

log(2) / 1.5 per day = 11.1 h, matching the paper’s “relatively short apparent half-life of free plasma PCSK9, estimated to be ~11 h”.

Validation 5 – the corrected HDL-C equation

Equation 6 as printed is

HDLc = HDLc_bl * (1 - lamtg * LDLc / Baseline_LDLc)

which returns (1 - 0.34) * HDLc_bl at baseline rather than HDLc_bl, and moves HDL-C down when LDL-C falls. The Methods describe HDL-C as “inversely related to plasma TGs”; equation 5 defines exactly such a baseline-referenced triglyceride ratio; and the Discussion states that “HDL-C increases by 5% to 10% in virtually every anti-PCSK9 trial”. The reading carried by the model, HDLc = HDLc_bl * (1 - lamtg * (TG - 1)), satisfies all three. The two readings are separated here by the paper’s own stated effect size.

hdl_check <- dplyr::bind_rows(lapply(seq_len(nrow(regimens)), function(i) {
  s <- sims[[i]]
  st <- interval_stats(s, "ldl", bl$BaselineLDLc, regimens$ii[i])
  u <- 1 + st[["trough"]] / 100 # LDL-C relative to baseline at its nadir
  tg <- bl$lamtgVLDL + bl$lamtgLDL * u
  tibble::tibble(
    Regimen = regimens$label[i],
    `LDL-C nadir (%)` = st[["trough"]],
    `HDL-C, corrected (%)` = 100 * (-bl$lamtg * (tg - 1)),
    `HDL-C, as printed (%)` = 100 * ((1 - bl$lamtg * u) - 1)
  )
}))
knitr::kable(hdl_check, digits = 1, caption = "Equation 6: the corrected reading reproduces the paper's stated +5% to +10% HDL-C rise; the printed reading gives a decrease of the wrong sign and magnitude.")
Equation 6: the corrected reading reproduces the paper’s stated +5% to +10% HDL-C rise; the printed reading gives a decrease of the wrong sign and magnitude.
Regimen LDL-C nadir (%) HDL-C, corrected (%) HDL-C, as printed (%)
Alirocumab 75 mg Q2W -56.4 4.2 -14.8
Alirocumab 300 mg Q4W -67.9 5.1 -10.9
Evolocumab 140 mg Q2W -73.1 5.5 -9.2
Evolocumab 420 mg Q4W -78.0 5.8 -7.5
Inclisiran 300 mg Q3M -50.2 3.8 -16.9 

stopifnot(
  # Corrected: within the Discussion's stated 5-10% band for every regimen with
  # a meaningful LDL-C drop.
  all(hdl_check$`HDL-C, corrected (%)` > 3),
  all(hdl_check$`HDL-C, corrected (%)` < 10),
  # As printed: wrong sign throughout.
  all(hdl_check$`HDL-C, as printed (%)` < 0)
)

Validation 6 – invariance to the unreported baselines

Supplemental Table S2 marks the baselines of PCSK9, LDL-C, Lp(a) and HDL-C “taken from each arm of each trial” and tabulates no pooled value. Only BaselineVLDLc = 23.166 mg/dL carries a number. BaselinePCSK9 is load bearing and is back-solved above; the other three are shipped as rounded placeholders. This section demonstrates that the published readouts are exactly invariant to those three, so the placeholders cannot bias any reproduced number.

alt <- rxode2::rxSolve(
  ui,
  local({
    ev <- rxode2::et(amt = 140, cmt = "depot_evolo", ii = 14, until = 252) |>
      rxode2::et(seq(0, 252, by = 0.05), cmt = "pcsk9")
    d <- as.data.frame(ev)
    d$TRT_ANTIPCSK9_SIRNA <- 0L
    d
  }),
  # Triple LDL-C, halve HDL-C, double Lp(a).
  params = c(BaselineLDLc = 300, BaselineHDLc = 25, BaselineLpA = 60),
  returnType = "data.frame", useLinCmt = FALSE, atol = 1e-10, rtol = 1e-10
)
ref <- sims[["Evolocumab 140 mg Q2W"]]

pct <- function(s, state, base) 100 * (s[[state]] / base - 1)
inv <- c(
  ldl = max(abs(pct(alt, "ldl", 300) - pct(ref, "ldl", bl$BaselineLDLc))),
  pcsk9 = max(abs(pct(alt, "pcsk9", bl$BaselinePCSK9) -
    pct(ref, "pcsk9", bl$BaselinePCSK9))),
  lpa = max(abs(pct(alt, "lpa", 60) - pct(ref, "lpa", bl$BaselineLpA))),
  HDLc = max(abs(pct(alt, "HDLc", 25) - pct(ref, "HDLc", bl$BaselineHDLc))),
  TG = max(abs(alt$TG - ref$TG)),
  ApoB = max(abs(alt$ApoB / alt$ApoB[1] - ref$ApoB / ref$ApoB[1]))
)
signif(inv, 3)
#>      ldl    pcsk9      lpa     HDLc       TG     ApoB 
#> 1.69e-08 1.18e-08 7.51e-09 1.27e-09 3.72e-11 1.52e-10
stopifnot(all(inv < 1e-6))

Percentage change from baseline in PCSK9, LDL-C, Lp(a), HDL-C, triglyceride and apoB is unchanged to solver precision when the three placeholder baselines are moved by factors of 3, 0.5 and 2. Algebraically: with VLDL-C held at its own baseline, the (ldl / BaselineLDLc, ldlr, pcsk9) subsystem is homogeneous in BaselineLDLc, and every downstream readout consumes only that ratio.

Total cholesterol and non-HDL-C are the two exceptions, because equations 7 and 8 add absolute mg/dL quantities and so depend on the baseline lipid mix. Their percentage changes should not be read off this model without supplying arm-appropriate baselines.

tc_shift <- max(abs(pct(alt, "TC", alt$TC[1]) - pct(ref, "TC", ref$TC[1])))
round(tc_shift, 2)
#> [1] 21.95
stopifnot(tc_shift > 1) # genuinely mix-dependent, as documented

Assumptions, deviations and errata

Corrections to printed equations

The source’s equation set is internally inconsistent in five places. Each correction below is forced by the paper’s own symmetry or by a numerical claim the paper itself makes; none was chosen to improve agreement with Table 1, and the Table 1 gate above was run only after the corrections were fixed by those arguments.

  1. Equation 16 is printed with a + where its = belongs (dcomE/dt [nmol] + kon...) and spells the drug subscript “evolvo”. Written here symmetric with equations 13 and 19, which are the well-formed alirocumab and RG-7652 complex balances.
  2. Equations 20 and 25, the PCSK9 balance, carry malformed evolocumab terms: the association term -konevolo * (Ac_evolo/Vd_evolo) * PCSK9 is absent, the dissociation term carries a minus instead of a plus, and a complex-clearance term -(CLevolo/Vdevolo) * comE appears that belongs to the comE balance of equation 16, not to the PCSK9 balance. The alirocumab and RG-7652 terms in the same equations are well formed and mutually symmetric; the repair makes evolocumab match them. The Table 1 gate confirms it: the repaired reading reproduces the printed evolocumab PCSK9 statistics to within 0.4 percentage points.
  3. Equation 24’s left-hand side is printed as dAc_inc/dt while every term on its right-hand side carries the aln subscript. It is the ALN-PCS circulating balance; equation 22 already gives the inclisiran one.
  4. Equation 6 is discussed in Validation 5 above. The printed form does not return the baseline at baseline and has the wrong sign.
  5. lamapoB is printed in supplemental Table S2 as -0.654 with the description “NonHDLc-to-ApoB conversion coefficient” and the estimation method “calculated”. A negative conversion coefficient makes equation 9 return a negative apoB mass concentration. The model carries +0.654, which makes apoB at baseline equal 0.654 x baseline non-HDL-C, a physiologically ordinary apoB-to-non-HDL-C ratio, and makes equation 9 self-consistent (lamapoBLDL + lamapoBVLDL = 1, so the bracket is 1 at baseline).

Values not printed in any source

  • BaselinePCSK9 = 6.3832 nmol is not printed anywhere – Table S2 marks it “taken from each arm of each trial” and tabulates no pooled value, Table S1 lists no baselines, and it does not appear in any figure panel (Figure 2 panels B and C are percentage change; panel A is drug amount). It is load bearing, because mAb binding is bimolecular and the depth of PCSK9 suppression depends on the molar ratio of antibody to target. It was back-solved from the paper’s own Table 1, by finding the single value that reproduces the eight printed mAb PCSK9 predose and trough entries; the fit is tight (all eight within 0.43 percentage points) and the value is then confirmed out of sample by the twenty LDL-C entries and the inclisiran row. 6.3832 nmol / 2.75 L = 2.32 nmol/L of plasma, an ordinary baseline plasma PCSK9. This is a derived value, not a printed one, and a user who has an arm’s own baseline should supply it.
  • BaselineLDLc, BaselineHDLc and BaselineLpA are shipped as the rounded placeholders 100, 50 and 30 mg/dL. Validation 6 shows every published readout except total cholesterol and non-HDL-C is exactly invariant to them.
  • MW_PCSK9 appears as a row of Table S2 with no value (“taken from the literature (21)”). It is never needed here: the model works in nmol throughout and reports percentage change, so no PCSK9 mass-to-mole conversion is performed. It is omitted from ini() rather than guessed.

Modelling scope

  • The model carries all five drugs in one 17-ODE system, as the authors built it. Dosing an arm means dosing that drug’s depot compartment; the other four chains stay at zero. The TRT_ANTIPCSK9_SIRNA covariate selects the modality-specific LDL-receptor exponent n1 (0.14 for the three mAbs, 0.26 for the two siRNAs), which is the entirety of the paper’s “two previously identified sets of parameters”. Co-administering an mAb and an siRNA is outside the model’s scope and the source never simulates it.
  • VLDL-C is held at its baseline and is not a state, which the Discussion states explicitly. Consequently the model cannot represent any intervention that changes VLDL-C, including the statin background therapy most of the calibration subjects received.
  • No IIV and no residual error. The 95% confidence intervals in Table S2 are parameter-uncertainty bands from the Fisher information matrix and likelihood profiling, not random-effect magnitudes, and are recorded in the ini() comments rather than encoded as etas. The model predicts an arm mean, not an individual profile.
  • The source notes two known misfits that this model inherits and that a user should not read as transcription errors: triglyceride and HDL-C responses are “more ambiguous” against the clinical data (the Results attribute this to TG metabolism being driven by non-LDL factors outside the model), and the 10 mg RG-7652 arm is underpredicted.

Session 

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.5 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
#>  [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
#>  [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
#> [10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   
#> 
#> time zone: UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] ggplot2_4.0.3         tidyr_1.3.2           dplyr_1.2.1          
#> [4] rxode2_5.1.8          nlmixr2lib_0.3.2.9000
#> 
#> loaded via a namespace (and not attached):
#>  [1] generics_0.1.4      sass_0.4.10         xml2_1.6.0         
#>  [4] digest_0.6.39       magrittr_2.0.5      RColorBrewer_1.1-3 
#>  [7] evaluate_1.0.5      grid_4.6.1          fastmap_1.2.0      
#> [10] lotri_1.0.5         jsonlite_2.0.0      whisker_0.4.1      
#> [13] rxode2ll_2.0.18     backports_1.5.1     purrr_1.2.2        
#> [16] scales_1.4.0        textshaping_1.0.5   jquerylib_0.1.4    
#> [19] cli_3.6.6           crayon_1.5.3        symengine_0.2.13   
#> [22] rlang_1.3.0         withr_3.0.3         cachem_1.1.0       
#> [25] yaml_2.3.12         otel_0.2.0          tools_4.6.1        
#> [28] parallel_4.6.1      memoise_2.0.1       checkmate_2.3.4    
#> [31] rxode2lincmt_0.1.0  vctrs_0.7.3         R6_2.6.1           
#> [34] lifecycle_1.0.5     fs_2.1.0            ragg_1.5.2         
#> [37] PreciseSums_0.7     fontawesome_0.5.3   pkgconfig_2.0.3    
#> [40] desc_1.4.3          rex_1.2.2           pkgdown_2.2.1      
#> [43] pillar_1.11.1       bslib_0.12.0        gtable_0.3.6       
#> [46] glue_1.8.1          data.table_1.18.6.1 Rcpp_1.1.2         
#> [49] systemfonts_1.3.2   tidyselect_1.2.1    xfun_0.61          
#> [52] tibble_3.3.1        sys_3.4.3           knitr_1.52         
#> [55] farver_2.1.2        dparser_1.3.1-13    htmltools_0.5.9    
#> [58] labeling_0.4.3      rmarkdown_2.32      compiler_4.6.1     
#> [61] S7_0.2.2            downlit_0.4.5       askpass_1.2.1      
#> [64] openssl_2.4.2