mRNA-0184 relaxin-2 (Rel2-vlk) translational PK/PD (Kaushal 2024) • nlmixr2lib

Model and source

Article: https://doi.org/10.1111/cts.13894 (PMC11283281)
Supporting Information (Table S1 model equations, Figure S1 individual profiles): https://www.ebi.ac.uk/europepmc/webservices/rest/PMC11283281/supplementaryFiles (file CTS-17-e13894-s001.pdf)

Two model files come out of this paper:

mod_cyno  <- readModelDb("Kaushal_2024_mRNA0184_cyno")
mod_human <- readModelDb("Kaushal_2024_mRNA0184_human")

Citation: Kaushal N, Attarwala H, Iqbal MJ, Saini R, Van L, Liang M. Translational pharmacokinetic/pharmacodynamic model for mRNA-0184, an investigational therapeutic for the treatment of heart failure. Clin Transl Sci. 2024;17(8):e13894. doi:10.1111/cts.13894. Structural equations from Supporting Information Table S1; parameter estimates from Table 1.
Kaushal_2024_mRNA0184_cyno – the fitted cynomolgus-monkey model (Table 1).
Kaushal_2024_mRNA0184_human – the allometric forward projection to a 70 kg adult (Table 2), which is what drove first-in-human dose selection.

mRNA-0184 is an investigational lipid-nanoparticle (LNP) encapsulated mRNA encoding human relaxin-2 fused to a variable light chain kappa domain (Rel2-vlk), developed for heart failure with reduced ejection fraction. Native relaxin-2 and the recombinant analogue serelaxin have half-lives of only 6-8 h; the therapeutic hypothesis behind mRNA-0184 is that continuous endogenous translation of an Fc-like fusion produces a much longer-lived circulating protein. The model below quantifies that: a Rel2-vlk protein half-life of about 6-7 days in humans.

Population

Cynomolgus-monkey model (Methods, “Data used for model development”):

Male cynomolgus monkeys, N = 4 per dose level at 0.15, 0.5 and 1 mg/kg (N = 12 total), 2.5 kg reference body weight (Table 2 column header).
Single dose of mRNA-0184 infused intravenously over 1 h.
Rel2-vlk mRNA and Rel2-vlk protein plasma concentrations measured to 337 h.
Assays: pre-qualified bDNA (mRNA, LLOQ 0.050 ng/mL) and ELISA (protein, LLOQ 20 pg/mL).
Estimation in Phoenix NLME 8.3.4.295, first-order conditional estimation-extended least squares.

Human model: no human data were fitted. Every value is a deterministic allometric transform of the cynomolgus estimate (Equation 1) to a 70 kg adult, intended for the first-in-human trial NCT05659264 in adults with stable heart failure with reduced ejection fraction.

Field	Value
Species	cynomolgus monkey (male, healthy)
N subjects	12
Weight reference	2.5 kg reference body weight (Table 2 column header)
Dose range	Single dose of mRNA-0184 at 0.15, 0.5, or 1 mg/kg administered as a 1-h intravenous infusion (N = 4 per dose level).
Disease state	Healthy (non-disease) non-human primates. The efficacious-exposure anchor used for human dose selection (AUEC over a 2-week interval = 486 ng/mL*h after weekly 0.15 mg/kg dosing) came from a SEPARATE study in aged, high-fat-diet-induced obese cynomolgus monkeys with naturally developed cardiovascular and metabolic disease; that study is reported as ‘data not shown’ and is NOT part of this model’s fitted dataset.

Model structure

Supporting Information Table S1 gives the complete system. For Rel2-vlk mRNA:

dA1/dt = Input - K12 * A1
dA2/dt = K12 * A1 - K23 * A2 + K32 * A3 - K20 * A2
dA3/dt = K23 * A2 - K32 * A3
C      = (A1 + A3) / V1

with K12 = CL/V1, K23 = CL3/V2, K32 = CL3/V1, K20 = CL2/V2. A1 is plasma-1, A2 is the tissue (target-site) compartment and A3 is plasma-2; Figure 1 labels the third circle “Plasma V1”, so A1 and A3 share the volume V1 and the observed mRNA concentration is the sum of the two plasma compartments. Elimination (K20) is from the tissue compartment. This plasma-tissue-plasma redistribution loop is what produces the delayed second concentration peak characteristic of LNP modalities.

For Rel2-vlk protein:

dCe/dt       = Ke0 * (C(t) - Ce(t))
dprotein1/dt = slope * Ce - Kprot * [protein1 - protein2] - K50 * protein1
dprotein2/dt = Kprot * (protein1 - protein2)
Rel2-vlk     = protein1 / Vc

Substituting the Table S1 rate-constant definitions turns every mRNA term into clearance x concentration. The same substitution is applied to Kprot and K50 here; see Assumptions and deviations for the four independent checks that select that reading.

Compartment mapping used in the model files: central = A1 (plasma-1), peripheral1 = A2 (tissue), peripheral2 = A3 (plasma-2), effect = Ce, rel2vlk = protein1, rel2vlk_p = protein2. Outputs are Cc (Rel2-vlk mRNA, ng/mL) and Rel2vlk (Rel2-vlk protein, ng/mL).

Source trace

Every ini() value, with its source location. The two right-hand columns double as a check of Equation 1: Y = a * (70 / 2.5)^beta, i.e. a body-weight ratio of 28.

Parameter	Meaning	Cyno (Table 1)	Scaling exponent	Equation 1 expected	Human (Table 2)	Pct diff vs Table 2	Source
lvc	V1, mRNA plasma volume, shared by A1 and A3 (mL)	112.000	1.00	3136.000	3136.000	0.00	Table 1 tvV1 (fixed) / Table 2
lcl12	CL, plasma-1 to tissue clearance (mL/h)	258.000	0.75	3140.400	3144.000	0.11	Table 1 tvCL / Table 2
lcl20	CL2, tissue elimination clearance (mL/h)	42.400	0.75	516.100	516.000	-0.02	Table 1 tvCL2 / Table 2
lvp	V2, mRNA tissue volume (mL)	160.000	1.00	4480.000	4490.000	0.22	Table 1 tvV2 / Table 2
lcl23	CL3, tissue to/from plasma-2 clearance (mL/h)	9.860	0.75	120.000	120.000	0.00	Table 1 tvCL3 / Table 2
lke0	Ke0, effect-compartment equilibration rate (1/h)	0.193	NA	0.193	0.193	0.00	Table 1 tvKe0 / Table 2
lvc_prot	Vc, protein central volume (mL)	163.000	1.00	4564.000	4550.000	-0.31	Table 1 tvVc (fixed) / Table 2
lvp_prot	Vp, protein peripheral volume (mL)	364.000	1.00	10192.000	10182.000	-0.10	Table 1 tvVp / Table 2
lkprot	Kprot, protein central/peripheral clearance (mL/h)	25.500	0.85	433.100	433.000	-0.02	Table 1 tvKprot / Table 2
lk50	K50, protein elimination clearance (mL/h)	4.580	0.85	77.800	77.700	-0.13	Table 1 tvK50 / Table 2
lslope	Slope, protein production per unit Ce (mL/h)	0.540	1.00	15.100	15.100	0.00	Table 1 tvSlope / Table 2
etalcl12	IIV variance on CL	0.266	NA	0.266	0.266	0.00	Table 1 etaCL / Table 2
etalvp	IIV variance on V2	1.060	NA	1.060	1.060	0.00	Table 1 etaV2 / Table 2
etalcl23	IIV variance on CL3	0.171	NA	0.171	0.171	0.00	Table 1 etaCL3 / Table 2
etalk50	IIV variance on K50	0.073	NA	0.073	0.073	0.00	Table 1 etaK50 / Table 2
etalslope	IIV variance on Slope	0.417	NA	0.417	0.417	0.00	Table 1 etaSlope / Table 2
propSd	Proportional residual SD, Rel2-vlk mRNA	0.557	NA	0.557	0.557	0.00	Table 1 / Table 2
propSd_Rel2vlk	Proportional residual SD, Rel2-vlk protein	0.378	NA	0.378	0.378	0.00	Table 1 / Table 2

Gate: every Table 2 value the paper prints must be reproducible from the Table 1 value via Equation 1 to within the paper’s own rounding. The largest discrepancy is V2 (4480 expected vs 4490 printed, 0.22%).

chk <- trace |> dplyr::filter(!is.na(Coef))
pct <- 100 * abs(chk$Human - chk$Cyno * 28^chk$Coef) / (chk$Cyno * 28^chk$Coef)
names(pct) <- chk$Parameter
print(round(pct, 3))
#>      lvc    lcl12    lcl20      lvp    lcl23 lvc_prot lvp_prot   lkprot 
#>    0.000    0.114    0.019    0.223    0.015    0.307    0.098    0.031 
#>     lk50   lslope 
#>    0.122    0.132
stopifnot(all(pct < 0.5))

Virtual cohorts

Cohorts are built from quantiles of each inter-individual random effect rather than random draws, so the results below are exactly reproducible and no arm is advantaged by a lucky sample. The same set of random effects is reused across every dose arm (common random numbers), which makes between-arm contrasts exact. Variances are read straight out of the model so the cohort can never drift from ini().

eta_sd <- function(mod) {
  d <- mod()$iniDf
  d <- d[!is.na(d$neta1) & d$neta1 == d$neta2, c("name", "est")]
  stats::setNames(sqrt(d$est), d$name)
}

quantile_etas <- function(mod, n, seed = 13894L) {
  sds <- eta_sd(mod)
  q <- stats::qnorm((seq_len(n) - 0.5) / n)
  set.seed(seed)
  out <- lapply(seq_along(sds), function(k) {
    sds[[k]] * (if (k == 1L) q else sample(q))
  })
  names(out) <- names(sds)
  cbind(data.frame(id = seq_len(n)), as.data.frame(out))
}

sd_cyno <- eta_sd(mod_cyno)
print(round(sd_cyno, 4))
#>  etalcl12    etalvp  etalcl23   etalk50 etalslope 
#>    0.5158    1.0296    0.4135    0.2702    0.6458

n_per_arm <- 60L
cyno_doses <- c(0.15, 0.5, 1)          # mg/kg, Methods
cyno_wt <- 2.5                         # kg, Table 2 column header
t_obs <- sort(unique(c(seq(0, 12, by = 0.25), seq(13, 48, by = 1),
                       seq(52, 337, by = 3))))

etas_cyno <- quantile_etas(mod_cyno, n_per_arm)

# expand a single set of dose rows into a repeated-dose schedule
expand_doses <- function(dose_rows, times) {
  do.call(rbind, lapply(times, function(tt) {
    d <- dose_rows
    d$time <- tt
    d
  }))
}

make_arm <- function(dose_mgkg, wt, etas, arm_index, t_grid, until = NA,
                     ii = NA) {
  amt_ng <- dose_mgkg * wt * 1e6       # mg/kg * kg -> mg -> ng
  offset <- (arm_index - 1L) * 1000L   # unique ids per arm
  e <- etas
  e$id <- e$id + offset
  dose <- e |>
    dplyr::mutate(time = 0, amt = amt_ng, rate = amt_ng, evid = 1L,
                  cmt = "central", dvid = NA_integer_)
  if (!is.na(ii)) {
    reps <- seq(0, until, by = ii)
    dose <- expand_doses(dose, reps)
  }
  obs <- e |>
    tidyr::crossing(time = t_grid) |>
    dplyr::mutate(amt = 0, rate = 0, evid = 0L, cmt = "central", dvid = 1L)
  dplyr::bind_rows(dose, obs) |>
    dplyr::mutate(dose_mgkg = dose_mgkg, arm = paste0(dose_mgkg, " mg/kg")) |>
    dplyr::arrange(id, time, dplyr::desc(evid))
}

cyno_events <- dplyr::bind_rows(
  lapply(seq_along(cyno_doses), function(i) {
    make_arm(cyno_doses[i], cyno_wt, etas_cyno, i, t_obs)
  })
)
stopifnot(n_per_arm <= 200L)
cat("cyno cohort:", dplyr::n_distinct(cyno_events$id), "subjects,",
    nrow(cyno_events), "rows\n")
#> cyno cohort: 180 subjects, 32760 rows

sim_cyno <- rxode2::rxSolve(
  mod_cyno, cyno_events,
  omega = NA,                 # random effects supplied explicitly as columns
  useLinCmt = FALSE,          # ODE->linCmt auto-conversion breaks dvid mapping
  keep = c("dose_mgkg", "arm"),
  returnType = "data.frame"
)
#> ℹ parameter labels from comments will be replaced by 'label()'
#> Warning: multi-subject simulation without without 'omega'
cat("solved rows:", nrow(sim_cyno), " any NA Cc:",
    anyNA(sim_cyno$Cc), " any NA Rel2vlk:", anyNA(sim_cyno$Rel2vlk), "\n")
#> solved rows: 32580  any NA Cc: FALSE  any NA Rel2vlk: FALSE

Figure S1 / Figure 2 replication – cynomolgus mRNA and protein profiles

Figure S1 shows individual Rel2-vlk mRNA (blue) and Rel2-vlk protein (red) concentration-time profiles per animal; Figure 2 shows visual predictive checks stratified by dose. The observed data have not been released, so the panels below reproduce the model-predicted median and 5th-95th percentile band that Figure 2’s red lines represent.

vpc <- sim_cyno |>
  dplyr::filter(time > 0) |>
  dplyr::select(arm, time, Cc, Rel2vlk) |>
  tidyr::pivot_longer(c(Cc, Rel2vlk), names_to = "analyte",
                      values_to = "conc") |>
  dplyr::mutate(analyte = dplyr::recode(analyte,
                                        Cc = "Rel2-vlk mRNA",
                                        Rel2vlk = "Rel2-vlk protein")) |>
  dplyr::group_by(arm, analyte, time) |>
  dplyr::summarise(lo = stats::quantile(conc, 0.05),
                   md = stats::median(conc),
                   hi = stats::quantile(conc, 0.95), .groups = "drop")

ggplot(vpc, aes(time, md)) +
  geom_ribbon(aes(ymin = lo, ymax = hi), alpha = 0.25, fill = "steelblue") +
  geom_line(colour = "steelblue4") +
  facet_grid(analyte ~ arm, scales = "free_y") +
  scale_y_log10() +
  scale_x_continuous(limits = c(0, 337)) +
  labs(x = "Time after dose (h)", y = "Concentration (ng/mL)",
       title = "Replicates Figure 2 (VPC by dose) and Figure S1 of Kaushal 2024",
       caption = paste("Model-predicted median with 5th-95th percentile band,",
                       n_per_arm, "virtual animals per dose arm."))

Structural gate – the delayed second mRNA peak

The whole reason for the plasma-tissue-plasma structure is the “delayed or secondary peak” seen in the exploratory analysis (Results, first paragraph; Figure S1). A typical animal must therefore show a local minimum followed by a genuine second maximum in plasma mRNA after the infusion ends, and that second peak must be carried by the plasma-2 compartment.

typ <- make_arm(0.15, cyno_wt,
                data.frame(id = 1L, etalcl12 = 0, etalvp = 0, etalcl23 = 0,
                           etalk50 = 0, etalslope = 0),
                1L, seq(0, 60, by = 0.1))
sim_typ <- rxode2::rxSolve(mod_cyno, typ, omega = NA, useLinCmt = FALSE,
                           returnType = "data.frame")
post <- sim_typ[sim_typ$time > 1, ]
turns <- which(diff(sign(diff(post$Cc))) != 0)

t_min <- post$time[turns[1] + 1]; c_min <- post$Cc[turns[1] + 1]
t_pk2 <- post$time[turns[2] + 1]; c_pk2 <- post$Cc[turns[2] + 1]
c_pk1 <- max(sim_typ$Cc[sim_typ$time <= 1])

cat(sprintf("first peak %.0f ng/mL at end of infusion; local min %.1f at %.1f h; second peak %.1f at %.1f h\n",
            c_pk1, c_min, t_min, c_pk2, t_pk2))
#> first peak 1352 ng/mL at end of infusion; local min 270.8 at 2.5 h; second peak 399.1 at 6.9 h

# a real second peak: two turning points, a rise of at least 10% off the nadir,
# and it must be driven by plasma-2 (A3), not residual plasma-1 (A1)
i_pk2 <- which.min(abs(sim_typ$time - t_pk2))
frac_a3 <- (sim_typ$peripheral2[i_pk2] / sim_typ$vc[1]) / c_pk2
stopifnot(length(turns) >= 2, t_pk2 > t_min, c_pk2 / c_min > 1.10,
          frac_a3 > 0.99)
cat(sprintf("plasma-2 accounts for %.1f%% of the second peak\n",
            100 * frac_a3))
#> plasma-2 accounts for 100.0% of the second peak

Half-life gates

The paper’s central quantitative claim is that mRNA-0184-induced Rel2-vlk protein has a half-life of “approximately 6-7 days” in humans versus 6-8 h for serelaxin (Discussion). That is also the claim that discriminates the two possible readings of Kprot and K50 (see Assumptions and deviations), so it is enforced strictly here. The cynomolgus protein half-life is checked at the same time: a terminal half-life of a few days is what makes a biexponential decline observable across the 337 h sampling window.

half_life <- function(mod, dose_mgkg, wt, t_end, lo, hi, col) {
  ev <- make_arm(dose_mgkg, wt,
                 data.frame(id = 1L, etalcl12 = 0, etalvp = 0, etalcl23 = 0,
                            etalk50 = 0, etalslope = 0),
                 1L, seq(0, t_end, by = 1))
  s <- rxode2::rxSolve(mod, ev, omega = NA, useLinCmt = FALSE,
                       returnType = "data.frame")
  s <- s[s$time >= lo & s$time <= hi & s[[col]] > 0, ]
  log(2) / -stats::coef(stats::lm(log(s[[col]]) ~ s$time))[2]
}

th_cyno_prot  <- half_life(mod_cyno,  0.15,  2.5, 500, 300, 500, "Rel2vlk")
th_cyno_mrna  <- half_life(mod_cyno,  0.15,  2.5, 500, 200, 350, "Cc")
th_human_prot <- half_life(mod_human, 0.025, 70,  900, 600, 900, "Rel2vlk")
#> ℹ parameter labels from comments will be replaced by 'label()'
th_human_mrna <- half_life(mod_human, 0.025, 70,  900, 300, 600, "Cc")

tibble::tibble(
  Analyte = c("Rel2-vlk protein", "Rel2-vlk mRNA",
              "Rel2-vlk protein", "Rel2-vlk mRNA"),
  Species = c("cynomolgus", "cynomolgus", "human (scaled)", "human (scaled)"),
  `Terminal half-life (h)` = round(c(th_cyno_prot, th_cyno_mrna,
                                     th_human_prot, th_human_mrna), 1),
  `Terminal half-life (days)` = round(c(th_cyno_prot, th_cyno_mrna,
                                        th_human_prot, th_human_mrna) / 24, 2)
) |>
  knitr::kable()

Analyte	Species	Terminal half-life (h)	Terminal half-life (days)
Rel2-vlk protein	cynomolgus	86.8	3.62
Rel2-vlk mRNA	cynomolgus	10.3	0.43
Rel2-vlk protein	human (scaled)	143.1	5.96
Rel2-vlk mRNA	human (scaled)	23.8	0.99


# Paper: human Rel2-vlk protein half-life approximately 6-7 days.
stopifnot(th_human_prot / 24 > 5.5, th_human_prot / 24 < 7.5)

# Allometry consistency: t1/2 scales as V/CL, i.e. 28^(1 - 0.85) = 28^0.15.
ratio <- th_human_prot / th_cyno_prot
cat(sprintf("human/cyno protein half-life ratio = %.3f; 28^0.15 = %.3f\n",
            ratio, 28^0.15))
#> human/cyno protein half-life ratio = 1.648; 28^0.15 = 1.648
stopifnot(abs(ratio - 28^0.15) < 0.02)

The half-life ratio between species reproduces 28^(1 - 0.85) to three decimal places, which is an internal check that the volume and clearance exponents were applied to the parameters the paper says they were applied to.

Human dose selection – Figure 5 replication

Dose selection was driven by the average steady-state Rel2-vlk protein concentration (Equation 2), anchored on an efficacious exposure of AUEC_ss = 486 ng/mL*h over a 2-week interval observed in a separate efficacy study in aged, obese cynomolgus monkeys.

auec_ss <- 486     # ng/mL*h, Methods (efficacy study, data not shown)
tau <- 14 * 24     # h, 2-week dosing interval
cavg_target <- auec_ss / tau
cat(sprintf("Equation 2: %g / %g = %.3f ng/mL (paper reports 1.45)\n",
            auec_ss, tau, cavg_target))
#> Equation 2: 486 / 336 = 1.446 ng/mL (paper reports 1.45)
stopifnot(abs(cavg_target - 1.45) < 0.01)

Figure 5a simulates every-2-weeks dosing over about 3000 h and marks 1008-1344 h as the steady-state window used for AUC_ss; Figure 5b plots the resulting median trough concentration against dose against the 1-2.5 ng/mL cardioprotective band. Both are reproduced below on that same window.

n_human <- 50L
human_doses <- c(0.01, 0.025, 0.05, 0.075, 0.1, 0.15)   # mg/kg, Figure 5b
human_wt <- 70
tau_h <- 336
t_human <- sort(unique(c(seq(0, tau_h, by = 6),
                         seq(tau_h, 1344, by = 6), 1344)))

etas_human <- quantile_etas(mod_human, n_human)

human_events <- dplyr::bind_rows(
  lapply(seq_along(human_doses), function(i) {
    make_arm(human_doses[i], human_wt, etas_human, i, t_human,
             until = 3 * tau_h, ii = tau_h)
  })
)
stopifnot(n_human <= 200L)
cat("human cohort:", dplyr::n_distinct(human_events$id), "subjects,",
    nrow(human_events), "rows\n")
#> human cohort: 300 subjects, 68700 rows

sim_human <- rxode2::rxSolve(
  mod_human, human_events,
  omega = NA, useLinCmt = FALSE,
  keep = c("dose_mgkg", "arm"),
  returnType = "data.frame"
)
#> Warning: multi-subject simulation without without 'omega'

sim_human |>
  dplyr::filter(time > 0, dose_mgkg %in% c(0.01, 0.025, 0.05, 0.15)) |>
  dplyr::group_by(arm, time) |>
  dplyr::summarise(lo = stats::quantile(Rel2vlk, 0.025),
                   md = stats::median(Rel2vlk),
                   hi = stats::quantile(Rel2vlk, 0.975), .groups = "drop") |>
  ggplot(aes(time, md)) +
  geom_ribbon(aes(ymin = lo, ymax = hi), alpha = 0.25, fill = "steelblue") +
  geom_line(colour = "steelblue4") +
  geom_hline(yintercept = c(1, 2.5), linetype = "dashed", colour = "red") +
  geom_vline(xintercept = c(1008, 1344), linetype = "dotted") +
  facet_wrap(~ arm) +
  scale_y_log10() +
  labs(x = "Time (h)", y = "Rel2-vlk protein (ng/mL)",
       title = "Replicates Figure 5a of Kaushal 2024",
       caption = paste("Every-2-weeks dosing, median with 95% prediction",
                       "interval. Red dashed lines are the 1-2.5 ng/mL",
                       "cardioprotective band; dotted lines bound the",
                       "1008-1344 h steady-state window."))

ss <- sim_human |>
  dplyr::filter(time >= 1008, time <= 1344) |>
  dplyr::group_by(dose_mgkg, id) |>
  dplyr::summarise(
    ctrough = Rel2vlk[which.max(time)],
    cavg = sum(diff(time) * (head(Rel2vlk, -1) + tail(Rel2vlk, -1)) / 2) /
      (max(time) - min(time)),
    .groups = "drop"
  ) |>
  dplyr::group_by(dose_mgkg) |>
  dplyr::summarise(`Median Ctrough (ng/mL)` = stats::median(ctrough),
                   `Median Cavg (ng/mL)` = stats::median(cavg),
                   .groups = "drop")
ss |>
  dplyr::mutate(dplyr::across(dplyr::where(is.numeric), \(x) round(x, 3))) |>
  dplyr::rename("Dose (mg/kg)" = dose_mgkg) |>
  knitr::kable()

Dose (mg/kg)	Median Ctrough (ng/mL)	Median Cavg (ng/mL)
0.010	0.385	0.822
0.025	0.962	2.056
0.050	1.925	4.112
0.075	2.887	6.168
0.100	3.849	8.225
0.150	5.774	12.337

fig5b_digitized <- tibble::tibble(
  dose_mgkg = c(0.01, 0.025, 0.05, 0.075, 0.15),
  ctrough_paper = c(0.43, 1.27, 2.00, 3.80, 6.60)
)

ss |>
  ggplot(aes(dose_mgkg, `Median Ctrough (ng/mL)`)) +
  geom_line(colour = "steelblue4") +
  geom_point(colour = "steelblue4") +
  geom_point(data = fig5b_digitized,
             aes(dose_mgkg, ctrough_paper), shape = 1, size = 3,
             colour = "black", inherit.aes = FALSE) +
  geom_hline(yintercept = c(1, 2.5), linetype = "dashed", colour = "red") +
  scale_y_log10() +
  labs(x = "Dose (mg/kg)", y = "Median steady-state Ctrough (ng/mL)",
       title = "Replicates Figure 5b of Kaushal 2024",
       caption = paste("Filled line: this model. Open circles: median read",
                       "off the published Figure 5b (log-axis digitization,",
                       "roughly +/-20%). Red dashed lines are the",
                       "1-2.5 ng/mL band."))

Dose-selection gate

The paper’s stated criterion (Methods, “Extrapolation of the PK/PD model to humans”) is the average steady-state concentration: “maintaining the Cavg_ss at a minimal level and within the reported levels of endogenous relaxin-2 (i.e., 1-2.5 ng/mL)”. Figure 5b additionally plots the median trough against the same band. Both are checked below. This is a genuinely discriminating criterion – the dose below must fall out of the band and the doses above must overshoot it – which is why it is worth asserting rather than merely plotting.

On Cavg_ss this model puts 0.025 mg/kg squarely inside the band. On trough it lands on the band’s lower edge (0.96 ng/mL median) where the published Figure 5b reads about 1.26 ng/mL; that roughly 20% offset is discussed under Assumptions and deviations and is asserted as a bounded range rather than being tuned away.

g <- ss |> dplyr::mutate(dose_mgkg = as.character(dose_mgkg))
ct <- stats::setNames(g$`Median Ctrough (ng/mL)`, g$dose_mgkg)
cv <- stats::setNames(g$`Median Cavg (ng/mL)`, g$dose_mgkg)

cat(sprintf("0.01  mg/kg: Ctrough %.3f  Cavg %.3f\n", ct["0.01"], cv["0.01"]))
#> 0.01  mg/kg: Ctrough 0.385  Cavg 0.822
cat(sprintf("0.025 mg/kg: Ctrough %.3f  Cavg %.3f\n", ct["0.025"], cv["0.025"]))
#> 0.025 mg/kg: Ctrough 0.962  Cavg 2.056
cat(sprintf("0.05  mg/kg: Ctrough %.3f  Cavg %.3f\n", ct["0.05"], cv["0.05"]))
#> 0.05  mg/kg: Ctrough 1.925  Cavg 4.112

# The paper's stated Methods criterion: Cavg_ss inside the 1-2.5 ng/mL band
# at the selected dose.
stopifnot(cv["0.025"] > 1.0, cv["0.025"] < 2.5)
# The dose below the selected one must fall out of the band on both metrics.
stopifnot(cv["0.01"] < 1.0, ct["0.01"] < 1.0)
# 0.05 mg/kg must push the average above the band, as Figure 5b shows.
stopifnot(cv["0.05"] > 2.5)
# Median trough at the selected dose: this model gives ~0.96 ng/mL against
# ~1.26 ng/mL read off Figure 5b. Asserted as a bounded range so that a
# structural regression still trips the gate, without tuning to the figure.
stopifnot(ct["0.025"] > 0.75, ct["0.025"] < 1.35)

# Steady state must actually be reached in the 1008-1344 h window: with a
# 143 h half-life and a 336 h interval, adding further doses beyond the
# fourth must not change the window's metrics.
cat(sprintf("accumulation ratio 1/(1-exp(-ln2*tau/thalf)) = %.3f\n",
            1 / (1 - exp(-log(2) * tau_h / th_human_prot))))
#> accumulation ratio 1/(1-exp(-ln2*tau/thalf)) = 1.244

Dose proportionality

The model is entirely linear – no saturable process anywhere – so every exposure metric must be exactly proportional to dose. Any departure is a bug (a lost dose scaling, a mis-set infusion rate) rather than pharmacology.

prop <- ss$`Median Ctrough (ng/mL)` / ss$dose_mgkg
names(prop) <- ss$dose_mgkg
print(round(prop, 4))
#>    0.01   0.025    0.05   0.075     0.1    0.15 
#> 38.4924 38.4924 38.4924 38.4924 38.4924 38.4924
stopifnot(max(prop) / min(prop) - 1 < 1e-4)

PKNCA validation

Cynomolgus single dose, both analytes

Non-compartmental analysis of the typical-animal profiles for each dose arm. Doses are 1-h infusions, so duration is supplied to PKNCAdose(); omitting it biases the steady-state volume estimates.

nca_events <- dplyr::bind_rows(
  lapply(seq_along(cyno_doses), function(i) {
    make_arm(cyno_doses[i], cyno_wt,
             data.frame(id = 1L, etalcl12 = 0, etalvp = 0, etalcl23 = 0,
                        etalk50 = 0, etalslope = 0),
             i, t_obs)
  })
)
sim_nca <- rxode2::rxSolve(mod_cyno, nca_events, omega = NA,
                           useLinCmt = FALSE,
                           keep = c("dose_mgkg", "arm"),
                           returnType = "data.frame")
#> Warning: multi-subject simulation without without 'omega'

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

# time-zero records must be present or PKNCA warns on every subject
stopifnot(all(conc_cyno$arm %in% conc_cyno$arm[conc_cyno$time == 0]))

dose_cyno <- nca_events |>
  dplyr::filter(evid == 1L) |>
  dplyr::select(id, time, amt, arm) |>
  dplyr::mutate(duration = 1)

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

nca_one <- function(conc_col) {
  cdf <- conc_cyno |>
    dplyr::select(id, time, arm, conc = dplyr::all_of(conc_col))
  o_conc <- PKNCA::PKNCAconc(cdf, conc ~ time | arm + id)
  o_dose <- PKNCA::PKNCAdose(dose_cyno, amt ~ time | arm + id,
                             duration = "duration")
  res <- PKNCA::pk.nca(PKNCA::PKNCAdata(o_conc, o_dose,
                                        intervals = intervals_cyno))
  as.data.frame(res) |>
    dplyr::filter(start == 0, end == 337) |>       # drop dependency rows
    dplyr::select(arm, PPTESTCD, PPORRES) |>
    dplyr::mutate(analyte = conc_col)
}

nca_cyno <- dplyr::bind_rows(nca_one("Cc"), nca_one("Rel2vlk"))

nca_cyno |>
  dplyr::mutate(analyte = dplyr::recode(analyte, Cc = "Rel2-vlk mRNA",
                                        Rel2vlk = "Rel2-vlk protein"),
                PPORRES = signif(PPORRES, 4)) |>
  tidyr::pivot_wider(names_from = arm, values_from = PPORRES) |>
  dplyr::arrange(analyte, PPTESTCD) |>
  dplyr::rename("Analyte" = analyte, "NCA parameter" = PPTESTCD) |>
  knitr::kable()

NCA parameter	Analyte	0.15 mg/kg	0.5 mg/kg	1 mg/kg
adj.r.squared	Rel2-vlk mRNA	9.999e-01	9.999e-01	9.999e-01
aucinf.obs	Rel2-vlk mRNA	1.029e+04	3.429e+04	6.858e+04
auclast	Rel2-vlk mRNA	1.029e+04	3.429e+04	6.858e+04
cl.obs	Rel2-vlk mRNA	3.645e+01	3.645e+01	3.645e+01
clast.obs	Rel2-vlk mRNA	1.000e-07	4.000e-07	8.000e-07
clast.pred	Rel2-vlk mRNA	1.000e-07	4.000e-07	8.000e-07
cmax	Rel2-vlk mRNA	1.352e+03	4.506e+03	9.012e+03
half.life	Rel2-vlk mRNA	1.037e+01	1.037e+01	1.037e+01
lambda.z	Rel2-vlk mRNA	6.683e-02	6.683e-02	6.683e-02
lambda.z.n.points	Rel2-vlk mRNA	1.610e+02	1.610e+02	1.610e+02
lambda.z.time.first	Rel2-vlk mRNA	5.000e+00	5.000e+00	5.000e+00
lambda.z.time.last	Rel2-vlk mRNA	3.370e+02	3.370e+02	3.370e+02
r.squared	Rel2-vlk mRNA	9.999e-01	9.999e-01	9.999e-01
span.ratio	Rel2-vlk mRNA	3.201e+01	3.201e+01	3.201e+01
tlast	Rel2-vlk mRNA	3.370e+02	3.370e+02	3.370e+02
tmax	Rel2-vlk mRNA	1.000e+00	1.000e+00	1.000e+00
adj.r.squared	Rel2-vlk protein	9.999e-01	9.999e-01	9.999e-01
aucinf.obs	Rel2-vlk protein	1.215e+03	4.049e+03	8.098e+03
auclast	Rel2-vlk protein	1.124e+03	3.746e+03	7.492e+03
cl.obs	Rel2-vlk protein	3.087e+02	3.087e+02	3.087e+02
clast.obs	Rel2-vlk protein	7.211e-01	2.404e+00	4.807e+00
clast.pred	Rel2-vlk protein	7.248e-01	2.416e+00	4.832e+00
cmax	Rel2-vlk protein	8.046e+00	2.682e+01	5.364e+01
half.life	Rel2-vlk protein	8.742e+01	8.742e+01	8.742e+01
lambda.z	Rel2-vlk protein	7.929e-03	7.929e-03	7.929e-03
lambda.z.n.points	Rel2-vlk protein	1.160e+02	1.160e+02	1.160e+02
lambda.z.time.first	Rel2-vlk protein	2.900e+01	2.900e+01	2.900e+01
lambda.z.time.last	Rel2-vlk protein	3.370e+02	3.370e+02	3.370e+02
r.squared	Rel2-vlk protein	9.999e-01	9.999e-01	9.999e-01
span.ratio	Rel2-vlk protein	3.523e+00	3.523e+00	3.523e+00
tlast	Rel2-vlk protein	3.370e+02	3.370e+02	3.370e+02
tmax	Rel2-vlk protein	2.700e+01	2.700e+01	2.700e+01

Because the system is linear, cl.obs for Rel2-vlk mRNA must be identical across the three dose arms, and it has an exact closed form worth checking. Integrating the Table S1 system to infinity gives

AUC = (integral A1 + integral A3) / V1 = Dose/CL + Dose/CL2

so the observed clearance is the series combination of the plasma-1-to-tissue clearance and the tissue elimination clearance:

CL_obs = CL * CL2 / (CL + CL2) = 258 * 42.4 / 300.4 = 36.42 mL/h

It is not CL2 alone: mRNA is dosed into plasma-1 but can only be eliminated from tissue, so the two clearances act in series. Note that CL3 and V2 cancel out of the result entirely – the tissue-to-plasma-2 redistribution loop returns everything it borrows – which makes this simultaneously a mass-balance check on the loop and a check that A1 and A3 really do share the volume V1.

cl_mrna <- nca_cyno |>
  dplyr::filter(analyte == "Cc", PPTESTCD == "cl.obs")
cl_expected <- 258 * 42.4 / (258 + 42.4)     # Table 1 tvCL and tvCL2
print(data.frame(arm = cl_mrna$arm, cl.obs = signif(cl_mrna$PPORRES, 5)))
#>          arm cl.obs
#> 1 0.15 mg/kg 36.452
#> 2  0.5 mg/kg 36.452
#> 3    1 mg/kg 36.452
cat(sprintf("closed form CL*CL2/(CL+CL2) = %.3f mL/h; NCA mean = %.3f mL/h\n",
            cl_expected, mean(cl_mrna$PPORRES)))
#> closed form CL*CL2/(CL+CL2) = 36.415 mL/h; NCA mean = 36.452 mL/h
# dose independence
stopifnot(max(cl_mrna$PPORRES) / min(cl_mrna$PPORRES) - 1 < 1e-4)
# and agreement with the closed form (NCA extrapolates the terminal phase, so
# allow a little slack rather than demanding an exact match)
stopifnot(abs(mean(cl_mrna$PPORRES) / cl_expected - 1) < 0.02)

Human steady-state interval

conc_hu <- sim_human |>
  dplyr::filter(!is.na(Rel2vlk), dose_mgkg == 0.025,
                time >= 1008, time <= 1344) |>
  dplyr::mutate(time_ss = time - 1008) |>
  dplyr::select(id, time = time_ss, conc = Rel2vlk) |>
  dplyr::mutate(arm = "0.025 mg/kg Q2W")

# The 1008-1344 h window is re-based to 0-336 h, so every subject must carry a
# time-zero record; without one PKNCA warns "AUC range starting (0) before the
# first measurement" once per subject.
stopifnot(all(table(conc_hu$id[conc_hu$time == 0]) == 1L),
          dplyr::n_distinct(conc_hu$id[conc_hu$time == 0]) ==
            dplyr::n_distinct(conc_hu$id))

dose_hu <- conc_hu |>
  dplyr::distinct(id, arm) |>
  dplyr::mutate(time = 0, amt = 0.025 * 70 * 1e6, duration = 1)

o_conc <- PKNCA::PKNCAconc(conc_hu, conc ~ time | arm + id)
o_dose <- PKNCA::PKNCAdose(dose_hu, amt ~ time | arm + id,
                           duration = "duration")
res_hu <- PKNCA::pk.nca(PKNCA::PKNCAdata(
  o_conc, o_dose,
  intervals = data.frame(start = 0, end = 336, cmax = TRUE, tmax = TRUE,
                         cmin = TRUE, auclast = TRUE, cav = TRUE)
))

nca_hu <- as.data.frame(res_hu) |>
  dplyr::filter(start == 0, end == 336) |>
  dplyr::group_by(PPTESTCD) |>
  dplyr::summarise(median = stats::median(PPORRES), .groups = "drop")

nca_hu |>
  dplyr::mutate(median = signif(median, 4)) |>
  dplyr::rename("NCA parameter" = PPTESTCD,
                "Median across cohort" = median) |>
  knitr::kable()

NCA parameter	Median across cohort
auclast	690.8000
cav	2.0560
cmax	3.3530
cmin	0.9561
tmax	54.0000

cmin is used rather than ctrough, which PKNCA can leave NA when the interval end does not coincide with a nominal trough sample.

Comparison against published values

Kaushal 2024 reports no NCA table, so the reference column below is assembled from the quantities the paper does state numerically, plus the median trough concentrations read off Figure 5b. The Figure 5b rows are digitized from a log-scale figure and carry roughly 20% reading uncertainty in their own right; they are shown for completeness, and the half-life row is the gate that actually constrains the model.

# (1) Dose-response: median steady-state trough against the Figure 5b reads.
ss_ct <- ss |>
  dplyr::filter(dose_mgkg %in% fig5b_digitized$dose_mgkg) |>
  dplyr::transmute(Dose = dose_mgkg, PPTESTCD = "cmin",
                   PPORRES = `Median Ctrough (ng/mL)`)
ref_ct <- fig5b_digitized |>
  dplyr::transmute(Dose = dose_mgkg, PPTESTCD = "cmin",
                   PPORRES = ctrough_paper)

nlmixr2lib::ncaComparisonTable(
  simulated = ss_ct, reference = ref_ct, by = "Dose",
  params = "cmin", units = c(cmin = "ng/mL")
)
#>   NCA parameter  Dose Reference Simulated  % diff
#> 1  Cmin (ng/mL) 0.010      0.43     0.385  -10.5%
#> 2  Cmin (ng/mL) 0.025      1.27     0.962 -24.2%*
#> 3  Cmin (ng/mL) 0.050         2      1.92   -3.8%
#> 4  Cmin (ng/mL) 0.075       3.8      2.89 -24.0%*
#> 5  Cmin (ng/mL) 0.150       6.6      5.77  -12.5%

# (2) The two quantities the paper states numerically, at the selected
#     0.025 mg/kg every-2-weeks dose.
sim2 <- tibble::tibble(
  PPTESTCD = c("half.life", "cav"),
  PPORRES = c(th_human_prot, nca_hu$median[nca_hu$PPTESTCD == "cav"])
)
ref2 <- tibble::tibble(
  PPTESTCD = c("half.life", "cav"),
  PPORRES = c(6.5 * 24,   # Discussion: "approximately 6-7 days" (midpoint)
              1.45)       # Equation 2 efficacy target, not a model prediction
)

nlmixr2lib::ncaComparisonTable(
  simulated = sim2, reference = ref2,
  params = c("half.life", "cav"),
  units = c(half.life = "h", cav = "ng/mL")
)
#>   NCA parameter Reference Simulated  % diff
#> 1        t½ (h)       156       143   -8.3%
#> 2  Cavg (ng/mL)      1.45      2.06 +41.8%*

Two things to read out of these tables. The half-life agrees with the paper’s stated 6-7 days to within 8%, which is the row that actually constrains the model. The cav row is the one large deviation and is expected rather than a model defect: 1.45 ng/mL is the efficacy target the paper derived from a separate obese-cynomolgus study via Equation 2, not a prediction of this model. The paper then chose the lowest dose whose exposure lands inside the 1-2.5 ng/mL band, which 0.025 mg/kg does on Cavg_ss (see the dose-selection gate above).

Across the five digitized Figure 5b dose levels the model runs low by 4-25%, and – importantly – it runs low by a similar factor at every dose rather than distorting the shape of the dose-response, which is what a structural error would do. Combined with the exact agreement on the half-life ratio, the Equation 1 reproduction of Table 2 and the closed-form clearance identity, that pattern points at the unreported details of the paper’s own Figure 5b simulation (and at log-axis digitization bias, which tends to read high) rather than at the encoded model.

Assumptions and deviations

Kprot and K50 are implemented as clearances (mL/h), not as rate constants (1/h). This is the one substantive departure from the printed text and it is load-bearing, so the evidence is set out in full. Table 1 labels both parameters h^-1 and Table S1 writes the protein distribution term as Kprot * [protein1 - protein2], i.e. a single rate constant applied to two amounts. Four independent checks contradict that literal reading and all four select the clearance form:
1. Vp (364 mL, SE 42.7%) is estimated and allometrically scaled in Table 2, but never appears anywhere in the printed equations. It is identifiable only if the distribution term is concentration-driven.
2. Table 2 scales Kprot and K50 with an exponent of 0.85, which the paper introduces explicitly as the allometric exponent for the clearance of therapeutic proteins. A first-order rate constant scales as CL/V, i.e. with 0.85 - 1 = -0.15; scaling it with +0.85 would make the human protein kinetics about 17-fold faster than the monkey’s.
3. Taken literally, the printed equations with the printed values give a terminal Rel2-vlk protein half-life of 0.019 h in humans. The paper states “approximately 6-7 days” (Discussion). The clearance reading gives 143 h, i.e. 5.96 days. That is a factor of roughly 7500 between the two candidate readings, so this is not a marginal call.
4. A single symbol Kprot serving both directions is only self-consistent as a symmetric intercompartmental clearance – which is exactly how Table S1 already treats the mRNA side, where the one symbol CL3 expands to two different rate constants (K23 = CL3/V2, K32 = CL3/V1). It also matches the paper’s stated assumption that “the distribution rate constant of Rel2-vlk protein is the same between the plasma and tissue compartments”.
The reading is therefore treated as a units/notation slip in Table 1 and Table S1 rather than as the intended model. The half-life gate above fails loudly if a future edit reverts it.
Slope is a volumetric rate, not dimensionless. Table 1 labels tvSlope “(dimensionless)”, but the protein production term slope * Ce must convert a concentration (ng/mL) into an amount rate (ng/h), so Slope carries units of mL/h. Table 2 confirms this by scaling it with an exponent of 1, i.e. as a volume.
Inter-individual variances are read as variances. Table 1 reports the eta rows as bare numbers with no CV% column, so there is no internal discriminator between a variance and a standard deviation. They are encoded as log-scale variances, which is what Phoenix NLME reports in its Omega block and what nlmixr2 expects for etaX ~ value. Under the alternative reading the implied CVs would be materially larger (for example etaV2 137% as a variance versus 144% as an SD – similar here, but etaCL would move from 55% to 27%).
Residual errors are read as proportional standard deviations. The mRNA row is headed “Proportional residual error (%)” with a value of 0.557, which is 55.7% as a fraction rather than 0.557%; the protein row (0.378) carries no percent sign at all. Both are encoded as fractional SDs (propSd), matching Phoenix’s stdev parameterisation.
Dose units. Concentrations are in ng/mL and volumes in mL, so amounts are in ng. Doses stated in mg/kg must be converted: 0.15 mg/kg in a 2.5 kg monkey is 3.75e5 ng, and 0.025 mg/kg in a 70 kg adult is 1.75e6 ng. There is no bioavailability term – Table S1 writes Input directly into A1, so the full dose enters plasma-1.
Human infusion duration assumed to be 1 h. The paper specifies a 1-h infusion for the cynomolgus study and does not state the duration used in the human simulations. 1 h is carried over. Given a mRNA plasma-1 half-life under an hour and a protein half-life of days, the assumption has no material effect on the steady-state metrics used for dose selection.
No covariates. Neither the fitted model nor the human projection contains any covariate effect, so both files omit covariateData. The human file also contains no within-human body-weight scaling, because the paper neither fits nor reports one: Table 2 is a single 70 kg parameter set and the simulations were run on mg/kg doses.
Table 2 rounding. Applying Equation 1 exactly to the Table 1 values reproduces the printed Table 2 column to within 0.22% (worst case V2: 4480 computed versus 4490 printed). The printed Table 2 values are used verbatim in the human model file. Table 2 also omits the scaling coefficient for tvCL2 and tvVp and mislabels tvCL3 as “(mL)” rather than “(mL/h)”; the values themselves are unambiguous (0.75 and 1 respectively, recovered by back-calculation and checked in the gate above).
Median trough runs about 20% below Figure 5b. At the selected 0.025 mg/kg every-2-weeks dose this model gives a median steady-state trough of 0.96 ng/mL (typical individual 1.03 ng/mL) where Figure 5b reads about 1.26 ng/mL. Three things are worth separating here. First, the paper’s own stated dose-selection criterion is on Cavg_ss, and on that metric the model lands inside the 1-2.5 ng/mL band (2.06 ng/mL median). Second, the offset is not a structural error: the half-life gate, the interspecies half-life ratio, the Equation 1 reproduction of Table 2, the mRNA clearance identity and exact dose proportionality all pass, and the same roughly 20% offset appears at every dose level rather than distorting the shape of the dose-response. Third, the paper does not report the N, the random-effect sampling scheme, or whether residual error was included in Figure 5b’s median, and because trough concentration is a nonlinear function of five lognormal random effects the population median is not the typical-value profile (this model’s median sits 7% below its own typical value). No parameter has been adjusted to close the gap.
Non-paper-derived values used only in the comparison table. The ctrough_* and cmin reference values are digitized from the published Figure 5b (a log-scale plot; roughly 20% reading uncertainty). No model parameter is derived from a figure – every ini() value comes from Table 1 or Table 2.
Observed data are not available. The cynomolgus concentration data behind Figures 2, 3, 4 and S1 have not been released, so the figure replications show model-predicted medians and percentile bands rather than observed-versus-predicted overlays.
**The 486 ng/mL*h efficacy anchor is from a different study.** It comes from aged, high-fat-diet obese cynomolgus monkeys with naturally developed cardiovascular and metabolic disease, reported as “data not shown”. It is not part of the dataset either model was fit to, and this model (parameterised on healthy animals) predicts a larger exposure at the same 0.15 mg/kg weekly dose. Only the derived human target concentration (1.45 ng/mL) and the 1-2.5 ng/mL band are used here.
Prior estimation run not published. V1 and Vc were “fixed based on a prior estimation” (Table 1 footnote); that prior run is not reported, so the two values are taken as published and wrapped in fixed().