Skip to contents

mRNA vaccine QSP (Selvaggio 2021)

Source: vignettes/articles/Selvaggio_2021_mRNA_vaccine_QSP.Rmd
Selvaggio_2021_mRNA_vaccine_QSP.Rmd

Model and source

  • Citation: Selvaggio G, Leonardelli L, Lofano G, Fresnay S, Parolo S, Medini D, Siena E, Marchetti L. A quantitative systems pharmacology approach to support mRNA vaccine development and optimization. CPT Pharmacometrics Syst Pharmacol. 2021;10(12):1448-1451. doi:10.1002/psp4.12721. Adaptive-immune backbone from Chen X, Hickling T, Vicini P. A mechanistic, multiscale mathematical model of immunogenicity for therapeutic proteins: Part 1-theoretical model. CPT Pharmacometrics Syst Pharmacol. 2014;3:e133. doi:10.1038/psp.2014.30. Innate-cell parameters fitted to macaque data from Liang F, et al. Efficient targeting and activation of antigen-presenting cells in vivo after modified mRNA vaccine administration in rhesus macaques. Mol Ther. 2017;25:2635-2647. doi:10.1016/j.ymthe.2017.08.006.
  • Article: https://doi.org/10.1002/psp4.12721
  • Description: QSP. Preclinical (rhesus macaque, calibration). Multicompartmental quantitative systems pharmacology model of mRNA vaccine immunogenicity: kinetics of a lipid nanoparticle (LNP) carrying mRNA at the injection site, uptake and antigen expression by dendritic cells (myeloid + plasmacytoid), monocytes and neutrophils, migration to the draining lymph node, T-cell and B-cell adaptive immune response, and antibody production. Adapted and extended from the Chen 2014 immunogenicity backbone (adaptive immune response retained unvaried; Selvaggio 2021 supplement ‘Model equations’). Innate-cell recruitment / uptake / migration parameters fitted to Liang 2017 rhesus macaque data (fluorescent-protein-encoding modified mRNA-LNP). Time in days; LNP in number of molecules; cell states in cell counts; antibody in pmole. Delivered as a mean-field extraction with a single mean-affinity B-cell subgroup (Chen 2014 encodes 17 affinity-subclones for the BCR distribution; the single-subgroup extraction with mean Ka is a documented deviation, see vignette Assumptions) and an explicit low-antibody closure for the receptor-competition equation (paper Eq for total antigen requires a Newton solve on Ag_f per timestep; not natively expressible in rxode2).

Selvaggio et al. (2021) present a Commentary in CPT: Pharmacometrics and Systems Pharmacology introducing a mechanistic quantitative systems pharmacology (QSP) framework for mRNA vaccine development. The framework adapts and extends the multiscale therapeutic-protein immunogenicity model of Chen, Hickling, and Vicini (2014) to represent the early events after injection of a lipid nanoparticle (LNP) carrying an mRNA cargo. The main text and Figure 1 present the two-compartment (injection site + draining lymph node) topology; the online Supplementary Material lists every ODE and every parameter value (Table S1).

The paper is a QSP application, not a population-PK / -PD analysis: there is no fitted between-subject variability, no residual error, and no subject-level fit table. The full parameter set (Table S1) is delivered as a single deterministic mean-field simulation calibrated to the rhesus macaque migration data of Liang et al. (2017); the paper’s principal result (Table 1) is a sensitivity analysis of the antibody titre at day 30 under 50 ug mRNA-LNP dosing.

Population

The QSP framework is deterministic and describes a “typical” innate + adaptive immune response to a single intramuscular injection of an mRNA-LNP vaccine. The innate-cell layer (mDC, pDC, monocytes, neutrophils) is fit to rhesus macaque cell-count data from Liang et al. (2017), where macaques received a fluorescent-protein-encoding modified mRNA-LNP and cell counts at the injection site and draining lymph node were quantified over time. The adaptive-immune layer (T cells, B cells, antibody) inherits parameter values from Chen et al. (2014), a therapeutic-protein immunogenicity backbone developed against human data. Selvaggio 2021 explicitly maintains the adaptive-immune layer “unvaried” while introducing new mRNA-specific early-events dynamics.

The model’s population metadata is available programmatically via readModelDb("Selvaggio_2021_mRNA_vaccine_QSP")()$population.

Source trace

Every parameter used by the model traces back to Selvaggio 2021 Supplementary Material Table S1. The provenance is recorded inline in inst/modeldb/therapeuticArea/Selvaggio_2021_mRNA_vaccine_QSP.R. Below is a consolidated summary organised by biological subsystem.

Subsystem Parameters Source
LNP kinetics k_deg_lnp Table S1 (fitted, Liang 2017)
Myeloid DC mdc0, k_dt_mdc, k_rc_mdc, k_up_mdc, k_exp_mdc, k_is2ln_mdc_lnp, k_is2ln_mdc_ag, k_ln2bl_mdc_ag, k_is2bl_mdc Table S1 (mixed: initial numbers ref [2], death rates ref [1], all others fitted)
Plasmacytoid DC k_dt_pdc, k_rc_pdc, k_up_pdc, k_exp_pdc, k_is2ln_pdc_lnp, k_is2ln_pdc_ag, k_ln2bl_pdc_ag, k_is2bl_pdc Table S1 (mixed; initial number = 0 hard-coded in model())
DC MHC / transport mhc0, k_tr Table S1, ref [1]
Monocytes mc0, k_dt_mc, k_rc_mc, k_up_mc, k_exp_mc, k_is2ln_mc_lnp, k_is2ln_mc_ag, k_ln2bl_mc_ag, k_is2bl_mc Table S1 (mixed)
Neutrophils k_dt_np, k_rc_np, k_up_np, k_exp_np, k_is2ln_np_lnp, k_is2ln_np_ag, k_ln2bl_np_ag, k_is2bl_np Table S1 (mixed; initial number = 0 hard-coded)
T cells nt0, k_dt_nt, k_dt_at, k_dt_mt, k_dt_ft, k_act_nt, k_act_mt, k_prol_at, f1_t Table S1, ref [1]
B cells nb0, brn, ka_mean, kr_bcr, cc_n, cc_m, k_prol_anb, k_prol_amb, k_dt_ab_b, k_dt_mb, k_dt_lp, k_dt_sp, g1_b, g2_b, k_act_nb, k_act_mb Table S1, ref [1] (mean-affinity subgroup extraction; see Assumptions)
Antibody k_prod_ab, k_deg_ab Table S1, ref [1]
Lymph-node volume v_ln Not in paper; canonical 0.5 mL default. See Assumptions.

The ODEs themselves are printed page-by-page in the Selvaggio 2021 Supplementary Material section “Model equations” (pp. 1-5 of the supplement PDF). Each of the 45 compartments in the packaged model corresponds to a state-variable equation in the supplement, and every term in every ODE line is a one-to-one transcription except the two documented deviations below.

Simulation

The model has 45 ODE states (LNP + 39 innate-cell states + 5 T-cell states + 7 B-cell / antibody states, with a single mean-affinity B-cell subgroup). Dosing is via amt into the lnp compartment, in number of LNP molecules. A 50 ug mRNA-LNP dose (as used in the Selvaggio 2021 sensitivity analysis, Table 1) converts to ~5.6e16 molecules for the 535 g/mol LNP MW listed in Table S1:

lnp_mw   <- 535           # g/mol, Table S1
n_avog   <- 6.022e23      # #molecules/mol
dose_ug  <- 50            # single 50 ug mRNA-LNP (Selvaggio 2021 Table 1)
dose_lnp <- dose_ug * 1e-6 / lnp_mw * n_avog
sprintf("50 ug mRNA-LNP -> %.3g LNP molecules", dose_lnp)
#> [1] "50 ug mRNA-LNP -> 5.63e+16 LNP molecules"
mod <- readModelDb("Selvaggio_2021_mRNA_vaccine_QSP")

Baseline (no-dose) hold

Before injecting the vaccine, the model should sit at its physiological baseline: naive mDC = 4288, naive T = 1445, naive B = 5200, and every other state (pDC, NP, LNP, antigen-processing intermediates, active / memory / functional cells, antibody) starting at zero. Simulating with amt = 0 for 30 days should keep every state at its initial value.

ev0 <- rxode2::et(amt = 0, cmt = "lnp", time = 0) |>
  rxode2::et(seq(0, 30, by = 1))
sim0 <- rxode2::rxSolve(mod, ev0)
message("Range of mdc_is over 30 days: [",
        format(min(sim0$mdc_is), digits = 4), ", ",
        format(max(sim0$mdc_is), digits = 4), "]")
#> Range of mdc_is over 30 days: [4288, 4288]
message("Range of Ab over 30 days: [",
        format(min(sim0$Ab), digits = 4), ", ",
        format(max(sim0$Ab), digits = 4), "]")
#> Range of Ab over 30 days: [0, 0]

At physiological baseline, mdc_is, nt, nb hold at their initial values within numerical precision, and the antibody state stays at zero. This confirms the steady-state closure k_br^X = X0 * (k_dt^X + k_IS2BL^X) (Selvaggio 2021 Table S1 note) is encoded consistently.

50 ug single-dose mRNA-LNP simulation

The following simulation reproduces the model’s response to a single 50 ug mRNA-LNP dose over 60 days:

ev50 <- rxode2::et(amt = dose_lnp, cmt = "lnp", time = 0) |>
  rxode2::et(seq(0, 60, by = 0.5))
sim50 <- rxode2::rxSolve(mod, ev50)

Innate cells at the injection site

The four innate-cell populations mDC / pDC / MC / NP are recruited to the injection site by LNP-driven chemokine signalling (k_rc^X), take up LNP (k_up^X), and either migrate to the draining lymph node carrying LNP or antigen or clear peripherally. The panel below shows the total mDC / pDC / MC / NP counts at the injection site as a function of time, aggregated across their naive / LNP-carrying / antigen-expressing / antigen-presenting sub-states. This is the model analogue of the paper’s supplementary Figure S1, upper row (“total cell recruitment at the IS”).

sim50 |>
  dplyr::mutate(
    mDC = mdc_is + mdc_is_lnp + mdc_is_ag + mdc_is_agl + mdc_is_agm + mdc_is_agh,
    pDC = pdc_is + pdc_is_lnp + pdc_is_ag + pdc_is_agl + pdc_is_agm + pdc_is_agh,
    MC  = mc_is  + mc_is_lnp  + mc_is_ag,
    NP  = np_is  + np_is_lnp  + np_is_ag
  ) |>
  dplyr::select(time, mDC, pDC, MC, NP) |>
  tidyr::pivot_longer(-time, names_to = "population", values_to = "cells") |>
  dplyr::mutate(population = factor(population, levels = c("mDC", "pDC", "MC", "NP"))) |>
  ggplot(aes(time, cells)) +
  geom_line() +
  facet_wrap(~population, scales = "free_y") +
  labs(x = "Time (days)", y = "Cell count at injection site",
       title = "Innate cells at the injection site",
       caption = "Model analogue of Selvaggio 2021 supplement Figure S1 (top row).")

Lymph-node antigen presentation

Following uptake and maturation, DCs migrate to the draining lymph node and present antigen at low / medium / high membrane-loading levels. Total antigen in the lymph node (weighted 0.1 / 0.5 / 0.9 across AgL / AgM / AgH) drives B-cell activation.

sim50 |>
  dplyr::mutate(
    Ag_total = 95.1e3 * (0.1 * (mdc_ln_agl + pdc_ln_agl) +
                          0.5 * (mdc_ln_agm + pdc_ln_agm) +
                          0.9 * (mdc_ln_agh + pdc_ln_agh))
  ) |>
  ggplot(aes(time, Ag_total)) +
  geom_line() +
  labs(x = "Time (days)", y = "Total antigen presented (#molecules)",
       title = "Total presented antigen in the draining lymph node",
       caption = "Ag = MHC0 * [0.1 * AgL + 0.5 * AgM + 0.9 * AgH] over mDC+pDC (Selvaggio 2021 supp).")

Adaptive immune response and antibody

The lymph-node antigen drives T-cell activation (naive T -> active T -> memory T + functional T) and B-cell activation (naive B -> active B -> memory B + short-/long-lived plasma cells -> antibody).

sim50 |>
  dplyr::mutate(
    Functional_T   = ft,
    Memory_T       = mt,
    Short_plasma   = sp,
    Long_plasma    = lp,
    Memory_B       = mb,
    Antibody_pmole = ab
  ) |>
  dplyr::select(time, Functional_T, Memory_T, Short_plasma, Long_plasma,
                Memory_B, Antibody_pmole) |>
  tidyr::pivot_longer(-time, names_to = "state", values_to = "value") |>
  ggplot(aes(time, value)) +
  geom_line() +
  facet_wrap(~state, scales = "free_y") +
  labs(x = "Time (days)", y = "Cell count (or pmole for antibody)",
       title = "Adaptive-immune response after 50 ug mRNA-LNP",
       caption = "Chen 2014 T/B/Ab kinetics driven by Selvaggio 2021 innate-cell antigen presentation.")

Dose-response sensitivity

Selvaggio 2021 Table 1 reports a sensitivity analysis of antibody titre at day 30 with respect to k_tr, k_up, k_rc and k_exp under a 50 ug injection. The table shows that increasing k_up (vector uptake) and k_rc (adjuvant-mediated recruitment) has the strongest impact on the antibody titre, with a “minor synergy beyond additive” between k_up and k_rc. As a proxy check on the model’s monotonic dose-response, we run three dose levels (10 / 50 / 250 ug):

doses_ug <- c(10, 50, 250)
ev_dr <- purrr::map2_dfr(seq_along(doses_ug), doses_ug, function(id, du) {
  d_lnp <- du * 1e-6 / lnp_mw * n_avog
  rxode2::et(amt = d_lnp, cmt = "lnp", time = 0, id = id) |>
    rxode2::et(seq(0, 60, by = 0.5)) |>
    as.data.frame() |>
    dplyr::mutate(id = id, dose_ug = du)
})
sim_dr <- rxode2::rxSolve(mod, ev_dr, keep = "dose_ug")
#> Warning: multi-subject simulation without without 'omega'

sim_dr |>
  ggplot(aes(time, Ab, group = factor(dose_ug), colour = factor(dose_ug))) +
  geom_line() +
  labs(x = "Time (days)", y = "Antibody titre (pmole)", colour = "Dose (ug)",
       title = "Antibody titre vs. time by mRNA-LNP dose",
       caption = "Selvaggio 2021 Table 1 reports the day-30 antibody titre as the sensitivity endpoint.") +
  scale_x_continuous(breaks = seq(0, 60, 10))

The model produces a monotonic dose-response in antibody titre over the 10-250 ug range, consistent with the paper’s sensitivity finding that higher LNP-driven recruitment / uptake amplifies the antibody response.

Day-30 antibody titre by dose 

sim_dr |>
  dplyr::filter(abs(time - 30) < 1e-6) |>
  dplyr::select(dose_ug, Ab) |>
  dplyr::arrange(dose_ug) |>
  dplyr::rename("Dose (ug mRNA-LNP)" = dose_ug,
                "Antibody at day 30 (pmole)" = Ab) |>
  knitr::kable(digits = 3,
               caption = "Day-30 antibody titre by dose (Selvaggio 2021 Table 1 sensitivity endpoint).")
Day-30 antibody titre by dose (Selvaggio 2021 Table 1 sensitivity endpoint).
Dose (ug mRNA-LNP) Antibody at day 30 (pmole)
10 0
50 0
250 0

Assumptions and deviations

The extraction packages the Selvaggio 2021 QSP model with the following documented deviations from a bit-for-bit reproduction. Each is noted inline in the model file and expanded here so downstream users are aware of the scope of the packaged deterministic simulator.

  1. Single-subgroup B-cell layer. Selvaggio 2021 (inheriting from Chen 2014) encodes 17 B-cell affinity-subclones (“J = 17”) with Ka_j spanning 3.91e-1 to 2.56e-4 pM^-1. Each subclone carries its own NB_i / ANB_i / AMB_i / MB_i / SP_i / LP_i / Ab_i state. The packaged model uses a single mean-affinity subgroup with ka_mean = 0.010005 /pM (geometric mean of the endpoints). Reason: 17 subclones times 7 states plus the transcendental receptor- competition equation for Ag_f are not natively expressible in rxode2 (see item 2), and the practical simulation cost of enumerating 119 additional ODEs is high enough to break the 5-minute vignette render budget. Downstream users who need the full affinity distribution can extend the shipped model by hand-authoring the 17-clone block; the mean-field version is sufficient for the sensitivity-style analyses that motivated the original paper.
  2. Explicit Ag_f closure. The paper’s receptor-competition equation Ag = Ag_f * [1 + 2 * sum_i(K_a_i * Ab_i / (1 + K_a_i * Ag_f)) + sum_i(K_a_i * BCR_i / (1 + K_a_i * Ag_f))] is implicit in Ag_f and requires a Newton solve per timestep. The packaged model uses the low-antibody closure Ag_f = Ag / (1 + 2 * K_a * Ab_conc + K_a * BCR_conc). This is exact in the limit K_a * Ag_f << 1 and progressively over-estimates Ag_f (and therefore ro, F_i, G_i, H_i) as the antibody / BCR loading approaches saturation. The approximation is documented here so users interpreting late-time antibody dynamics understand the direction of the systematic bias.
  3. LNP mass-balance interpretation. The Selvaggio 2021 supplement writes the LNP loss equation as dLNP/dt = -k_up^X * (sum of all IS states of cell type X) * ... - k_deg^LNP * LNP (summing across every subtype of every cell), which conflicts with the individual-cell equations that source LNP uptake only from the naive form (k_up^X * LNP * X_IS_naive). The packaged model takes the mass-balance-consistent interpretation: naive cells are the only sink for LNP, matching the mDC_IS -> mDC_IS^LNP transition rate in the individual-cell equations. This resolves an inconsistency in the published supplement.
  4. Lymph-node volume default. The receptor-competition closure converts per-cell BCR / Ab counts to lymph-node concentrations using a lymph-node volume v_ln. Neither Selvaggio 2021 nor Chen 2014 tabulates this constant; the packaged model uses a canonical 0.5 mL reference volume as a fixed literature default. Users who need to match a specific anatomy can override v_ln at simulation time.
  5. Immature-DC death rate not carried. Table S1 lists k_dt^ImDC = 0.0924 /day and k_dt^IpDC = 0.0924 /day for immature dendritic cells; the paper’s ODEs use only the mature-cell death rate k_dt^mDC = k_dt^pDC = 0.02310 /day throughout. The packaged model correspondingly carries only the mature rates, matching the ODEs.
  6. Nominal residual placeholder. The nlmixr2 UI observation declaration requires an estimable residual-error parameter; the packaged model uses addSd_Ab = 0.01 (pmole) as a placeholder. The paper reports no residual-error model (the QSP is deterministic). For typical-value simulations, suppress the residual via mod |> rxode2::zeroRe().
  7. amt/mta rename. The paper’s active-memory-T state (AMT) is called mta in the packaged model because amt is a reserved rxode2 dose-amount keyword. All other state names transliterate directly.

Errata

  • Selvaggio 2021 supplement Table S1 lists k_IS2BL^pDC = 02605 /day; the missing decimal point reads as 0.02605 /day in context (matching the per-day order of the neighbouring k_IS2BL^mDC = 0.6603 /day and k_IS2BL^MC = 2.5318e-5 /day) and is entered as such.
  • The B-cell IIV distribution “NB_0 is assumed to be normally distributed within the 17 subgroups” is stated in Table S1’s leading note but only the group mean (5200 cells) is published; the per-subgroup variance is not reported. This does not affect the mean-field extraction (see Assumption 1).