Skip to contents

Benzathine benzylpenicillin G, intramuscular in children (Hand 2019)

Source: vignettes/articles/Hand_2019_benzathine_benzylpenicillin_g.Rmd
Hand_2019_benzathine_benzylpenicillin_g.Rmd

Model and source

  • Citation: Hand RM, Salman S, Newall N, Vine J, Page-Sharp M, Bowen AC, Gray K, Baker A, Kado J, Joseph J, Marsh J, Ramsay J, Sika-Paotonu D, Batty KT, Manning L, Carapetis J. A population pharmacokinetic study of benzathine benzylpenicillin G administration in children and adolescents with rheumatic heart disease: new insights for improved secondary prophylaxis strategies. J Antimicrob Chemother. 2019;74(7):1984-1991. doi:10.1093/jac/dkz076. The fixed-elimination, absorption-limited structural approach is shared with the same group’s later analyses; see modellib(‘Kado_2020_benzathine_benzylpenicillin_g’) and modellib(‘Kado_2023_benzathine_benzylpenicillin_g’).
  • Description: One-compartment population PK model for benzylpenicillin released from an intramuscular depot of benzathine benzylpenicillin G (Bicillin L-A) with two sequential first-order absorption stages (t1/2,abs-1 ~ 0.46 day for dissolution of the crystal suspension, t1/2,abs-2 ~ 8.9 day for the rate-limiting release into plasma) and a fixed, allometrically scaled elimination rate constant. Fat-free mass scales the apparent volume of distribution, and a body mass index at or above 25 kg/m^2 increases the slow absorption half-life by 86.5%, producing flip-flop kinetics in which absorption, not clearance, determines the observed terminal half-life. Developed from 256 dried-blood-spot benzylpenicillin concentrations collected over six monthly injection cycles in 18 children and adolescents receiving secondary prophylaxis for rheumatic heart disease (Hand 2019).
  • Article: https://doi.org/10.1093/jac/dkz076

Monthly intramuscular benzathine benzylpenicillin G (BPG) has been the backbone of rheumatic heart disease secondary prophylaxis since the 1950s, on the premise that it holds plasma benzylpenicillin above 0.02 mg/L – the group A Streptococcus susceptibility breakpoint – for most of the interval between injections. Hand 2019 is the first population pharmacokinetic study of that regimen in the population it is actually prescribed to: children and adolescents with established rheumatic heart disease, sampled by dried blood spot across six consecutive monthly injection cycles.

The paper’s two structural findings are what make it worth packaging:

  • Elimination is not identifiable and was fixed. Benzylpenicillin’s own elimination half-life is 20-60 minutes, so nothing observed over a 28-day cycle reflects clearance. Fitting kel freely returned elimination half-lives “much longer than previously reported”; the authors therefore fixed kel at 1.32 h^-1 per 70 kg with theory-based allometric exponents and let the absorption model carry the shape of the curve. Everything observed is flip-flop kinetics.
  • Two sequential first-order absorption stages, parameterised by their half-lives (Figure 1): a fast stage (t1/2,abs-1 ~ 0.46 day), interpreted as dissolution of the BPG crystal suspension in the depot, and a slow stage (t1/2,abs-2 ~ 8.9 day), interpreted as hydrolysis of BPG to benzathine and penicillin and release into plasma. The slow stage is rate-limiting and therefore is the observed terminal half-life.

The single retained covariate on absorption is a body mass index at or above 25 kg/m^2, which increases t1/2,abs-2 by 86.5% – nearly doubling it. The authors’ mechanistic reading is that in heavier participants the “intramuscular” injection is often not intramuscular at all but intra-adipose or subcutaneous, and that the resulting slower release flattens the concentration-time profile. Fat-free mass is the size descriptor for the apparent volume of distribution.

Population

Eighteen children and adolescents contributed 256 benzylpenicillin concentrations (Table 1; 22 were enrolled and 4 withdrew before contributing samples). Median age was 14.1 years (range 7.9-17.7), median weight 62.9 kg (range 29.9-149), median height 1.61 m (range 1.36-2.05) and median body mass index 23.6 kg/m^2 (range 16.2-44.4). Eight (44%) were male and eight (44%) had a BMI at or above 25 kg/m^2, which is the covariate stratification used throughout the paper (n = 10 lower BMI, n = 8 higher BMI). Fourteen participants identified as Aboriginal, two as Maori and two as Samoan. All received the full 900 mg (1.2 MIU) Bicillin L-A dose by intramuscular injection into the upper outer gluteal quadrant every 28 days, alternating sides.

Samples were dried blood spots assayed by LC-MS/MS (LLOQ 0.0025 mg/L, LOD 0.001 mg/L); 25 of 256 concentrations (9.7%) were below the limit of quantification and retained by the M3 likelihood method. Creatinine was not measured, but no participant had known renal disease. Estimation was by NONMEM 7.2.0 with LAPLACIAN INTER on natural-log-transformed concentrations.

The same information is available programmatically via the model’s population metadata (readModelDb("Hand_2019_benzathine_benzylpenicillin_g")()$population).

Source trace

Every ini() entry in inst/modeldb/specificDrugs/Hand_2019_benzathine_benzylpenicillin_g.R carries an in-file comment naming its source location. They are collected here for review.

Equation / parameter Value Source location
lkel (fixed) 1.32 h^-1 per 70 kg = 31.68 day^-1 Table 3, “kel (h^-1 70 kg^-1)”; bootstrap column reads “fixed”
lvc 72.2 L per 70 kg Table 3, “V (L 70 kg^-1)”; bootstrap 72.0 (64.0-84.2)
e_ffm_vc (fixed) 1 Results: kel fixed with an exponent of -1/4, “equivalent to an exponential of 3/4 for CL and 1 for V”
e_ffm_kel (fixed) -0.25 Results: “the elimination rate constant was fixed with allometric scaling with an exponential of -1/4”
lthalf_abs1 0.455 day Table 3, “t1/2, abs-1 (days)”; bootstrap 0.461 (0.174-0.948)
lthalf_abs2 8.88 day Table 3, “t1/2, abs-2 (days)”; bootstrap 8.79 (5.71-12.5)
e_bmi_thalf_abs2 0.865 Table 3, “increase in t1/2, abs-2 with BMI >= 25 kg/m^2 (%)” = 86.5; bootstrap 86.8 (33.4-198)
etalvc omega^2 = 0.26^2 = 0.0676 Table 3, “IIV in V” = 26% (shrinkage 9%); footnote: values are 100 x sqrt(variance)
etalthalf_abs1 omega^2 = 0.78^2 = 0.6084 Table 3, “IIV in t1/2, abs-1” = 78% (shrinkage 12%)
etalthalf_abs2 omega^2 = 0.63^2 = 0.3969 Table 3, “IIV in t1/2, abs-2” = 63% (shrinkage 12%)
etaiov_lthalf_abs2_1 omega^2 = 0.30^2 = 0.09 Table 3, “IOV in t1/2, abs-2” = 30% (shrinkage 46%)
corr(t1/2,abs-1, t1/2,abs-2) -1, fixed -> encoded as -0.99 Table 3, “r (t1/2, abs-1, t1/2, abs-2)” = -1, bootstrap “fixed”; see Assumptions
corr(t1/2,abs-2, V) -0.746 Table 3, “r (t1/2, abs-2, V)”; bootstrap -0.808 (-1.00 to -0.316)
corr(t1/2,abs-1, V) +0.746 (derived) Forced by r = -1 above; see Assumptions
propSd 0.35 Table 3, “RV (%)” = 35 (shrinkage 13%); bootstrap 34 (30-38). Additive on the natural-log scale = proportional on the linear scale
d/dt(depot), d/dt(transit1), d/dt(central) n/a Figure 1: Bolus ->k_a-1-> Absorption ->k_a-2-> V ->k_el->
FFM as the size descriptor reference 70 kg Results: “Fat-free mass was the best size parameter for allometric scaling on V”; Table 3 column headings “per 70 kg”
BMI >= 25 kg/m^2 threshold 25 kg/m^2 Results: “BMI as a categorical variable, with a threshold of >= 25 kg/m^2, resulted in the best fit” 
mod <- readModelDb("Hand_2019_benzathine_benzylpenicillin_g")
mod
#> function() {
#>   description <- "One-compartment population PK model for benzylpenicillin released from an intramuscular depot of benzathine benzylpenicillin G (Bicillin L-A) with two sequential first-order absorption stages (t1/2,abs-1 ~ 0.46 day for dissolution of the crystal suspension, t1/2,abs-2 ~ 8.9 day for the rate-limiting release into plasma) and a fixed, allometrically scaled elimination rate constant. Fat-free mass scales the apparent volume of distribution, and a body mass index at or above 25 kg/m^2 increases the slow absorption half-life by 86.5%, producing flip-flop kinetics in which absorption, not clearance, determines the observed terminal half-life. Developed from 256 dried-blood-spot benzylpenicillin concentrations collected over six monthly injection cycles in 18 children and adolescents receiving secondary prophylaxis for rheumatic heart disease (Hand 2019)."
#>   reference <- paste(
#>     "Hand RM, Salman S, Newall N, Vine J, Page-Sharp M, Bowen AC, Gray K,",
#>     "Baker A, Kado J, Joseph J, Marsh J, Ramsay J, Sika-Paotonu D, Batty KT,",
#>     "Manning L, Carapetis J. A population pharmacokinetic study of benzathine",
#>     "benzylpenicillin G administration in children and adolescents with",
#>     "rheumatic heart disease: new insights for improved secondary prophylaxis",
#>     "strategies. J Antimicrob Chemother. 2019;74(7):1984-1991.",
#>     "doi:10.1093/jac/dkz076.",
#>     "The fixed-elimination, absorption-limited structural approach is shared",
#>     "with the same group's later analyses; see",
#>     "modellib('Kado_2020_benzathine_benzylpenicillin_g') and",
#>     "modellib('Kado_2023_benzathine_benzylpenicillin_g').",
#>     sep = " "
#>   )
#>   vignette <- "Hand_2019_benzathine_benzylpenicillin_g"
#>   units <- list(time = "day", dosing = "mg", concentration = "mg/L")
#> 
#>   # Doses are expressed as milligrams of the benzathine benzylpenicillin G salt
#>   # (Methods, 'Clinical study procedures': "Bicillin L-A ... 2.3 mL containing
#>   # 900 mg (1.2 MIU) of BPG"), and the measured analyte is benzylpenicillin
#>   # (Methods, 'Measuring penicillin from DBS'). V is therefore an apparent
#>   # volume V/F that absorbs both the salt-to-penicillin mass conversion and the
#>   # unknown bioavailable fraction; no separate F is estimated. Concentrations
#>   # were assayed in dried blood spots but the paper analysed and reports them as
#>   # plasma concentrations (Methods, 'Pharmacokinetic modelling and
#>   # simulations': "Log_e plasma concentration-time datasets for
#>   # benzylpenicillin"), so `central` is annotated as plasma.
#>   compartmentData <- list(
#>     depot = list(
#>       analyte = "benzathine benzylpenicillin g",
#>       units = "mg",
#>       specimen = "administration site",
#>       verified = TRUE
#>     ),
#>     transit1 = list(
#>       analyte = "benzathine benzylpenicillin g",
#>       units = "mg",
#>       specimen = "administration site",
#>       verified = TRUE
#>     ),
#>     central = list(
#>       analyte = "benzylpenicillin",
#>       units = "mg",
#>       specimen = "plasma",
#>       verified = TRUE
#>     )
#>   )
#> 
#>   covariateData <- list(
#>     FFM = list(
#>       description = "Fat-free mass, the allometric size descriptor for the apparent volume of distribution and for the fixed elimination rate constant.",
#>       units = "kg",
#>       type = "continuous",
#>       reference_category = NULL,
#>       notes = paste(
#>         "Reference 70 kg: Table 3 reports kel as 'h^-1 . 70 kg^-1' and V as",
#>         "'L . 70 kg^-1'. Results: 'Fat-free mass was the best size parameter for",
#>         "allometric scaling on V.' The same size descriptor also carries the",
#>         "elimination rate constant. The paper does not restate that explicitly,",
#>         "but it is forced by its own equivalence claim: kel was fixed 'with",
#>         "allometric scaling with an exponential of -1/4 (equivalent to an",
#>         "exponential of 3/4 for CL and 1 for V)'. Because CL = kel * V, the",
#>         "exponents -1/4 and 1 sum to 3/4 only if kel and V are scaled by the SAME",
#>         "size descriptor; scaling V by fat-free mass and kel by total body weight",
#>         "would not produce a clean 3/4-power clearance and would contradict the",
#>         "paper's stated equivalence. Both are therefore scaled by FFM against a",
#>         "70 kg reference. Absorption parameters are NOT allometrically scaled.",
#>         "The paper derived FFM from body weight and body mass index using the",
#>         "published model of Anderson BJ, Holford NH, Drug Metab Pharmacokinet",
#>         "2009;24:25-36 (reference 20; Methods, 'Pharmacokinetic modelling and",
#>         "simulations': 'Fat-free mass was estimated from weight and BMI from a",
#>         "published model in children'). That reference reparameterises the",
#>         "Janmahasatian et al. equation already registered for this column;",
#>         "FFM = WHSmax * HT^2 * WT / (WHS50 * HT^2 + WT) is algebraically",
#>         "identical to WHSmax * WT / (WHS50 + BMI), with WHSmax = 42.92 and",
#>         "WHS50 = 30.93 for males and WHSmax = 37.99 and WHS50 = 35.98 for",
#>         "females. Hand 2019 itself prints none of those constants; a downstream",
#>         "user must either supply FFM directly or compute it from WT, HT and SEXF.",
#>         "Assumed time-fixed at baseline.",
#>         sep = " "
#>       ),
#>       source_name = "FFM"
#>     ),
#>     BMI = list(
#>       description = "Body mass index at baseline; dichotomised inside model() at the 25 kg/m^2 threshold the paper selected, and applied to the slow absorption half-life t1/2,abs-2.",
#>       units = "kg/m^2",
#>       type = "continuous",
#>       reference_category = "BMI < 25 kg/m^2 (the reference stratum, n = 10 of 18)",
#>       notes = paste(
#>         "The model consumes BMI as a CONTINUOUS column and forms the binary",
#>         "indicator internally, because the paper's covariate is explicitly a",
#>         "thresholded version of a continuous measurement: 'Although many body",
#>         "composition covariates were correlated with k_a-2, BMI as a categorical",
#>         "variable, with a threshold of >=25 kg/m^2, resulted in the best fit and",
#>         "was associated with an 86.5% increase in t1/2,abs-2' (Results). Applied",
#>         "as the linear multiplier t1/2,abs-2 = theta * (1 + 0.865 * [BMI >= 25]).",
#>         "BMI is also the second input to the fat-free-mass derivation described",
#>         "under FFM, so a downstream dataset that carries BMI and WT can produce",
#>         "both covariate columns. The threshold is the conventional",
#>         "overweight cut-off and is the same 25 kg/m^2 boundary the group used to",
#>         "stratify enrolment in its later phase 1 study",
#>         "(modellib('Kado_2023_benzathine_benzylpenicillin_g')). Assumed",
#>         "time-fixed at baseline. No other covariate relationship reached",
#>         "significance (Results: 'No other significant covariate relationships",
#>         "were identified').",
#>         sep = " "
#>       ),
#>       source_name = "BMI"
#>     )
#>   )
#> 
#>   population <- list(
#>     species = "human",
#>     n_subjects = 18L,
#>     n_studies = 1L,
#>     age_range = "7.9-17.7 years",
#>     age_median = "14.1 years",
#>     weight_range = "29.9-149 kg",
#>     weight_median = "62.9 kg",
#>     height_range = "1.36-2.05 m",
#>     height_median = "1.61 m",
#>     bmi_range = "16.2-44.4 kg/m^2",
#>     bmi_median = "23.6 kg/m^2",
#>     sex_female_pct = 56,
#>     race_ethnicity = c(Aboriginal = 77.8, Maori = 11.1, Samoan = 11.1),
#>     disease_state = "children and adolescents with a history of acute rheumatic fever or established rheumatic heart disease, receiving monthly benzathine benzylpenicillin G as secondary prophylaxis; none had known established renal disease",
#>     dose_range = "900 mg (1.2 MIU) benzathine benzylpenicillin G by intramuscular injection into the upper outer gluteal quadrant, alternating sides, once every 28 days; all 18 analysed participants received the full 900 mg dose",
#>     regions = "Australia (metropolitan Perth, Western Australia)",
#>     notes = paste(
#>       "Longitudinal observational study conducted March-November 2017 through the",
#>       "Princess Margaret Hospital ambulatory care service. 22 participants were",
#>       "enrolled and 4 withdrew before contributing samples; 18 contributed 256",
#>       "benzylpenicillin concentrations, 16 (89%) with full datasets across all six",
#>       "monthly injection cycles. Samples were dried blood spots assayed by a",
#>       "validated LC-MS/MS method (LLOQ 0.0025 mg/L, LOD 0.001 mg/L); 25 (9.7%)",
#>       "concentrations were below the limit of quantification and were retained",
#>       "using the M3 likelihood method of Beal 2001. Intensive sampling on days 1,",
#>       "3, 6, 12 and 21 after the injection was performed in two of the six cycles,",
#>       "with an additional trough before each injection and unscheduled samples",
#>       "triggered by sore throat. Creatinine was not measured. Estimation was by",
#>       "NONMEM 7.2.0 LAPLACIAN with INTER on log_e-transformed concentrations.",
#>       "Baseline demographics: Table 1. Final estimates and bootstrap: Table 3.",
#>       "Individual post hoc absorption half-lives and exposure metrics by BMI",
#>       "stratum: Table 2. Model schematic: Figure 1. Goodness of fit: Figure 2.",
#>       "Prediction-corrected VPC stratified by BMI: Figure 3. Growth-chart-based",
#>       "dosing simulations: Figure 4.",
#>       sep = " "
#>     )
#>   )
#> 
#>   ini({
#>     # ==========================================================================
#>     # Disposition
#>     # ==========================================================================
#>     # kel was FIXED, not estimated: Table 3 gives the bootstrap column for kel as
#>     # "fixed", and Results explains why -- "Initial analysis using standard
#>     # compartmental modelling with various absorption models resulted in
#>     # estimates of elimination t1/2 that were much longer than previously
#>     # reported for benzylpenicillin. Therefore, the elimination rate constant was
#>     # fixed with allometric scaling with an exponential of -1/4 ... based on
#>     # previously published data in children receiving intravenous
#>     # benzylpenicillin." Table 3 reports it per hour; converted to 1/day because
#>     # every absorption parameter in this model is expressed in days. The same
#>     # 1.32 h^-1 . 70 kg^-1 value is carried forward by the group's later
#>     # analysis, modellib('Kado_2023_benzathine_benzylpenicillin_g').
#>     lkel <- fixed(log(1.32 * 24))
#>     label("Elimination rate constant at 70 kg fat-free mass (1/day)") # Table 3: kel = 1.32 h^-1 . 70 kg^-1, bootstrap 'fixed' -> 31.68 1/day
#> 
#>     lvc <- log(72.2)
#>     label("Apparent central volume of distribution V/F at 70 kg fat-free mass (L)") # Table 3: V = 72.2 L . 70 kg^-1 (bootstrap 72.0, 95% CI 64.0-84.2)
#> 
#>     # Allometric exponents were imposed a priori, not estimated (Results: the
#>     # exponent of -1/4 on kel is "equivalent to an exponential of 3/4 for CL and
#>     # 1 for V"), so both are fixed.
#>     e_ffm_vc <- fixed(1)
#>     label("Allometric exponent of fat-free mass on V/F (unitless), applied a priori") # Results: 'equivalent to an exponential of 3/4 for CL and 1 for V'
#>     e_ffm_kel <- fixed(-0.25)
#>     label("Allometric exponent of fat-free mass on kel (unitless), applied a priori") # Results: kel fixed 'with allometric scaling with an exponential of -1/4'
#> 
#>     # ==========================================================================
#>     # Absorption
#>     # ==========================================================================
#>     # Figure 1: Bolus --k_a-1--> Absorption --k_a-2--> V --k_el-->. Two
#>     # sequential FIRST-ORDER stages; Results states "First-order absorption for
#>     # both these stages performed better than models with zero-order process"
#>     # and "The addition of peripheral compartment(s) did not improve the model."
#>     # Both are parameterised by half-life, which model() converts to a rate
#>     # constant as k = log(2) / t1/2.
#>     lthalf_abs1 <- log(0.455)
#>     label("Fast absorption half-life from the injection depot, t1/2,abs-1 (day)") # Table 3: t1/2,abs-1 = 0.455 days (bootstrap 0.461, 95% CI 0.174-0.948)
#>     lthalf_abs2 <- log(8.88)
#>     label("Slow absorption half-life into the central compartment, t1/2,abs-2 (day)") # Table 3: t1/2,abs-2 = 8.88 days (bootstrap 8.79, 95% CI 5.71-12.5)
#> 
#>     # Linear fractional increase in t1/2,abs-2 for the higher-BMI stratum.
#>     # Estimated (the bootstrap column reports a median and a CI), not fixed.
#>     e_bmi_thalf_abs2 <- 0.865
#>     label("Fractional increase in t1/2,abs-2 when BMI is at or above 25 kg/m^2 (unitless)") # Table 3: 'increase in t1/2,abs-2 with BMI >=25 kg/m2 (%)' = 86.5 (bootstrap 86.8, 95% CI 33.4-198)
#> 
#>     # ==========================================================================
#>     # Inter-individual variability
#>     # ==========================================================================
#>     # Table 3 footnote: 'IIV, IOV and RV are presented as 100% x sqrt(variability
#>     # estimate)'. The tabulated percentage is therefore 100 x omega, so
#>     # omega^2 = (percentage / 100)^2 directly -- NOT log(1 + CV^2). Results:
#>     # 'A full covariance matrix model was used', so all three etas form one
#>     # block.
#>     #
#>     # Off-diagonals. The paper prints two of the three correlations:
#>     #   r(t1/2,abs-1, t1/2,abs-2) = -1     (Table 3, bootstrap column 'fixed')
#>     #   r(t1/2,abs-2, V)          = -0.746 (bootstrap -0.808, 95% CI -1.00 to -0.316)
#>     # The third is forced by the first: with r12 exactly -1 the two absorption
#>     # etas are perfectly collinear, eta1 = -(sd1/sd2) * eta2, so
#>     #   r(t1/2,abs-1, V) = -r(t1/2,abs-2, V) = +0.746.
#>     #
#>     # A correlation of exactly -1 makes the 3x3 block singular and rxode2's
#>     # Cholesky sampler cannot decompose it. Following the repository convention
#>     # for published perfect correlations, the OFF-DIAGONALS are scaled by 0.99
#>     # (so r12 becomes -0.99); the diagonal variances -- the published IIV values
#>     # -- are untouched. The resulting matrix is positive definite (smallest
#>     # eigenvalue 4.8e-03).
#>     #
#>     #   var(t1/2,abs-1) = 0.78^2 = 0.6084
#>     #   var(t1/2,abs-2) = 0.63^2 = 0.3969
#>     #   var(V)          = 0.26^2 = 0.0676
#>     #   cov(abs-1, abs-2) = 0.99 * (-1.000) * 0.78 * 0.63 = -0.486486
#>     #   cov(abs-1, V)     = 0.99 * (+0.746) * 0.78 * 0.26 = +0.149776
#>     #   cov(abs-2, V)     = 0.99 * (-0.746) * 0.63 * 0.26 = -0.120973
#>     etalthalf_abs1 + etalthalf_abs2 + etalvc ~ c(
#>       0.6084,
#>       -0.486486, 0.3969,
#>       0.149776, -0.120973, 0.0676
#>     )
#> 
#>     # ==========================================================================
#>     # Inter-occasion variability
#>     # ==========================================================================
#>     # Table 3: 'IOV in t1/2,abs-2' = 30% (shrinkage 46%; bootstrap 31, 95% CI
#>     # 20-48); omega^2 = 0.30^2 = 0.09. An occasion is one monthly injection
#>     # cycle, and the source data span six of them. rxode2 parses but cannot
#>     # simulate the native 'eta ~ var | occ' multi-level syntax, so this is
#>     # encoded as a single occasion-indexed eta (the repository's registered
#>     # etaiov_<param>_<occasion> form) that is drawn once per subject per solve.
#>     # A single-cycle simulation reproduces the source exactly; a multi-cycle
#>     # simulation reuses the occasion-1 draw across every cycle and therefore
#>     # under-represents within-subject cycle-to-cycle variation. See the
#>     # vignette's Assumptions and deviations section.
#>     etaiov_lthalf_abs2_1 ~ 0.09
#> 
#>     # ==========================================================================
#>     # Residual variability
#>     # ==========================================================================
#>     # Concentrations were analysed on the natural-log scale (Methods,
#>     # 'Pharmacokinetic modelling and simulations': 'Log_e plasma
#>     # concentration-time datasets ... were analysed'), so the additive
#>     # log-scale residual is a proportional error model on the linear scale.
#>     propSd <- 0.35
#>     label("Proportional residual error (fraction)") # Table 3: RV = 35% (shrinkage 13%; bootstrap 34, 95% CI 30-38)
#>   })
#> 
#>   model({
#>     # ---- Reference values ----------------------------------------------------
#>     ffmRef <- 70 # kg      -- Table 3 reports kel and V per 70 kg
#>     bmiCut <- 25 # kg/m^2  -- Results: 'BMI as a categorical variable, with a threshold of >=25 kg/m2'
#> 
#>     # ---- Covariate dichotomisation ------------------------------------------
#>     # The paper's covariate is the indicator [BMI >= 25 kg/m^2]; the model takes
#>     # continuous BMI and applies the published threshold here.
#>     bmiHigh <- 0
#>     if (BMI >= bmiCut) {
#>       bmiHigh <- 1
#>     }
#> 
#>     # ---- Disposition ---------------------------------------------------------
#>     # Volume carries the estimated typical value, the a priori exponent of 1 on
#>     # fat-free mass, and the only disposition IIV. kel is a fixed rate constant
#>     # carrying the a priori exponent of -1/4 on the same size descriptor and no
#>     # IIV, so clearance kel * V scales with the 3/4 power of fat-free mass, which
#>     # is exactly the equivalence the Results section asserts.
#>     vc <- exp(lvc + etalvc) * (FFM / ffmRef)^e_ffm_vc
#>     kel <- exp(lkel) * (FFM / ffmRef)^e_ffm_kel
#> 
#>     # ---- Absorption ----------------------------------------------------------
#>     # t1/2,abs-1 carries IIV only. t1/2,abs-2 carries IIV, IOV, and the linear
#>     # BMI effect applied to the typical value.
#>     # The IOV term is applied on a second line so that the mu-referenced line
#>     # carries exactly one subject-level random effect, as rxode2 requires.
#>     thalf_abs1 <- exp(lthalf_abs1 + etalthalf_abs1)
#>     thalf_abs2Base <- exp(lthalf_abs2 + etalthalf_abs2)
#>     thalf_abs2 <- thalf_abs2Base * exp(etaiov_lthalf_abs2_1) *
#>       (1 + e_bmi_thalf_abs2 * bmiHigh)
#> 
#>     ka1 <- log(2) / thalf_abs1
#>     ka2 <- log(2) / thalf_abs2
#> 
#>     # ---- Structure (Figure 1) ------------------------------------------------
#>     # depot ("Bolus") --ka1--> transit1 ("Absorption") --ka2--> central (V) --kel-->
#>     d/dt(depot) <- -ka1 * depot
#>     d/dt(transit1) <- ka1 * depot - ka2 * transit1
#>     d/dt(central) <- ka2 * transit1 - kel * central
#> 
#>     Cc <- central / vc
#>     Cc ~ prop(propSd)
#>   })
#> }
#> <environment: 0x55f5e8c14e10>

Structural checks against the typical-value model

Before simulating a cohort, three properties that follow deterministically from the published parameters are checked against a typical-value solve (rxode2::zeroRe()). These are numerical identities, not statistical comparisons, so they are asserted tightly: any of them failing means a parameter, an exponent or an equation was mis-transcribed.

dose_mg <- 900 # Methods: Bicillin L-A, 900 mg (1.2 MIU) per injection
ffm_ref <- 70 # Table 3: kel and V are reported per 70 kg
ffm_typ <- 45 # kg; a representative fat-free mass, see the cohort section below

mod_typ <- rxode2::zeroRe(mod)

typical_profile <- function(bmi) {
  ev <- rxode2::et(amt = dose_mg, cmt = "depot") |>
    rxode2::et(seq(0, 200, by = 0.05), cmt = "central")
  dat <- as.data.frame(ev)
  dat$FFM <- ffm_typ
  dat$BMI <- bmi
  out <- as.data.frame(rxode2::rxSolve(mod_typ, dat, returnType = "data.frame"))
  out[out$time > 0, c("time", "Cc")]
}

typ_low <- typical_profile(20) # BMI < 25 stratum
#> ℹ omega/sigma items treated as zero: 'etalthalf_abs1', 'etalthalf_abs2', 'etalvc', 'etaiov_lthalf_abs2_1'
typ_high <- typical_profile(30) # BMI >= 25 stratum
#> ℹ omega/sigma items treated as zero: 'etalthalf_abs1', 'etalthalf_abs2', 'etalvc', 'etaiov_lthalf_abs2_1'

# Terminal slope over the last 60 days, where the profile is log-linear.
terminal_half_life <- function(prof) {
  tail_dat <- prof[prof$time > 140 & prof$Cc > 0, ]
  log(2) / -stats::coef(stats::lm(log(Cc) ~ time, data = tail_dat))[["time"]]
}

# Trapezoidal AUC over the full 200-day window (Clast is < 1e-9 of Cmax, so
# this is AUC0-inf to well within the assertion tolerance).
trapz_auc <- function(prof) {
  sum(diff(prof$time) * (utils::head(prof$Cc, -1) + utils::tail(prof$Cc, -1)) / 2)
}

kel_typ <- 1.32 * 24 * (ffm_typ / ffm_ref)^-0.25
vc_typ <- 72.2 * (ffm_typ / ffm_ref)^1
auc_expected <- dose_mg / (kel_typ * vc_typ)

structural <- data.frame(
  check = c(
    "Terminal t1/2, BMI < 25 (day)",
    "Terminal t1/2, BMI >= 25 (day)",
    "Ratio of terminal t1/2 (high / low BMI)",
    "AUC0-inf, BMI < 25 (mg*day/L)",
    "AUC0-inf, BMI >= 25 (mg*day/L)"
  ),
  expected = c(8.88, 8.88 * 1.865, 1.865, auc_expected, auc_expected),
  observed = c(
    terminal_half_life(typ_low),
    terminal_half_life(typ_high),
    terminal_half_life(typ_high) / terminal_half_life(typ_low),
    trapz_auc(typ_low),
    trapz_auc(typ_high)
  )
)
structural$pct_diff <- 100 * (structural$observed - structural$expected) /
  structural$expected
knitr::kable(structural, digits = 4)
check expected observed pct_diff
Terminal t1/2, BMI < 25 (day) 8.8800 8.8800 0.0000
Terminal t1/2, BMI >= 25 (day) 16.5612 16.5612 0.0000
Ratio of terminal t1/2 (high / low BMI) 1.8650 1.8650 0.0000
AUC0-inf, BMI < 25 (mg*day/L) 0.5481 0.5480 -0.0078
AUC0-inf, BMI >= 25 (mg*day/L) 0.5481 0.5479 -0.0280

The first three rows confirm the flip-flop claim: the terminal half-life of the packaged model is not log(2)/kel (about 32 minutes) but t1/2,abs-2 itself, and the higher-BMI stratum reproduces the published 86.5% increase exactly. The last two rows confirm mass balance: total exposure equals Dose / (kel . V) and is untouched by the absorption covariate, which only redistributes that exposure over time.

stopifnot(
  # Numerical identities from the published parameters; no cohort sampling and
  # no random draws are involved, so these are asserted tightly.
  all(abs(structural$pct_diff) < 0.5)
)

Virtual cohort

Original observed data are not publicly available. The cohort below matches the Table 1 demographics, split into the paper’s two analysis strata with 200 simulated participants each.

Fat-free mass is the model’s size covariate. Hand 2019 derived it “from weight and BMI from a published model in children” citing Anderson & Holford (Drug Metab Pharmacokinet 2009;24:25-36), which is a reparameterisation of the Janmahasatian et al. equation registered for the FFM column in inst/references/covariate-columns.md:

FFM = WHSmax . HT^2 . WT / (WHS50 . HT^2 + WT), with WHSmax = 42.92 and WHS50 = 30.93 for males and WHSmax = 37.99 and WHS50 = 35.98 for females. Hand 2019 prints none of those constants itself; see Assumptions and deviations.

rxode2::rxSetSeed(20190401)
set.seed(20190401)

n_per_arm <- 200 # 200 per arm is the repository cap and is ample for a VPC

# Truncated sampler so every covariate stays inside the Table 1 observed range.
rtrunc_norm <- function(n, mean, sd, lo, hi) {
  out <- stats::rnorm(n, mean, sd)
  while (any(bad <- out < lo | out > hi)) {
    out[bad] <- stats::rnorm(sum(bad), mean, sd)
  }
  out
}

make_arm <- function(arm, bmi_lo, bmi_hi, bmi_mean, bmi_sd, id_offset) {
  data.frame(
    id = id_offset + seq_len(n_per_arm),
    arm = arm,
    # Table 1: 8 of 18 (44%) male, so SEXF = 1 with probability 0.56.
    SEXF = stats::rbinom(n_per_arm, 1, 0.56),
    # Table 1: median height 1.61 m, range 1.36-2.05.
    HT = rtrunc_norm(n_per_arm, 1.61, 0.16, 1.36, 2.05),
    BMI = rtrunc_norm(n_per_arm, bmi_mean, bmi_sd, bmi_lo, bmi_hi)
  )
}

cohort <- bind_rows(
  # Table 1: BMI median 23.6, range 16.2-44.4; 10 of 18 below the 25 threshold.
  make_arm("BMI < 25", 16.2, 24.99, 21.0, 2.8, 0L),
  make_arm("BMI >= 25", 25.0, 44.4, 30.5, 5.5, 1000L)
) |>
  mutate(
    WT = BMI * HT^2,
    # Janmahasatian / Anderson-Holford fat-free mass; see the text above.
    WHSMAX = ifelse(SEXF == 1, 37.99, 42.92),
    WHS50 = ifelse(SEXF == 1, 35.98, 30.93),
    FFM = WHSMAX * HT^2 * WT / (WHS50 * HT^2 + WT)
  ) |>
  select(id, arm, SEXF, HT, BMI, WT, FFM)

cohort |>
  group_by(arm) |>
  summarise(
    n = n(),
    across(c(HT, BMI, WT, FFM), ~ stats::median(.x)),
    .groups = "drop"
  ) |>
  dplyr::rename(
    "Stratum" = arm,
    "N" = n,
    "Height (m)" = HT,
    "BMI (kg/m2)" = BMI,
    "Weight (kg)" = WT,
    "FFM (kg)" = FFM
  ) |>
  knitr::kable(digits = 1, caption = "Median simulated covariates by stratum.")
Median simulated covariates by stratum.
Stratum N Height (m) BMI (kg/m2) Weight (kg) FFM (kg)
BMI < 25 200 1.6 20.9 55.4 39.8
BMI >= 25 200 1.6 32.1 83.8 50.5

Simulated weights (median about 55 kg in the lower-BMI arm and 83 kg in the higher-BMI arm) sit inside the published 29.9-149 kg range, and the higher-BMI arm is heavier – matching the paper’s observation that the higher-BMI participants “also had lower weights” reversed, i.e. that the lower-BMI stratum was lighter.

Which dosing interval does Table 2 describe?

Table 2 of Hand 2019 reports each participant’s Cmin, Cmax, Tmax and time above target, but the Methods do not say whether those per-participant metrics cover a single injection cycle in isolation or a cycle at steady state after months of monthly dosing. The distinction matters here: the model’s slow absorption half-life is 8.88 days in the lower-BMI stratum but 16.56 days in the higher-BMI stratum against a 28-day interval, so accumulation is modest in one stratum and substantial in the other. The two readings are therefore distinguishable, and the paper’s own numbers settle it.

Both are solved below at the typical-value parameters, with fat-free mass taken from the Table 1 medians rather than from a sampled cohort, so that the only thing varying between the two candidate answers is the number of preceding doses.

# Fat-free mass implied by the Table 1 medians (weight 62.9 kg, BMI 23.6
# kg/m^2), sex-weighted at the Table 1 ratio of 8 of 18 male.
ffm_male <- 42.92 * 62.9 / (30.93 + 23.6)
ffm_female <- 37.99 * 62.9 / (35.98 + 23.6)
ffm_table1 <- 0.44 * ffm_male + 0.56 * ffm_female

tau <- 28 # day; Methods: monthly injections

# Solve a 28-day cycle preceded by `n_cycles - 1` earlier monthly injections,
# for one lower-BMI and one higher-BMI typical individual. The BMI values only
# have to fall on either side of the published 25 kg/m^2 threshold; the model
# uses the indicator, not the magnitude.
typical_cycle <- function(n_cycles) {
  t0 <- (n_cycles - 1L) * tau
  ev <- rxode2::et(
    amt = dose_mg,
    time = seq(0, by = tau, length.out = n_cycles),
    cmt = "depot",
    id = 1:2
  ) |>
    rxode2::et(seq(t0, t0 + tau, by = 0.25), cmt = "central")
  dat <- as.data.frame(ev)
  dat$FFM <- ffm_table1
  dat$BMI <- ifelse(dat$id == 1, 23.6, 26.0)
  out <- as.data.frame(rxode2::rxSolve(mod_typ, dat, returnType = "data.frame"))
  out |>
    filter(time >= t0) |>
    mutate(
      tad = time - t0,
      arm = ifelse(id == 1, "BMI < 25", "BMI >= 25")
    ) |>
    select(arm, tad, Cc)
}

cycle_metrics <- function(prof) {
  prof |>
    group_by(arm) |>
    summarise(
      `Cmax (ug/L)` = max(Cc) * 1000,
      `Cmin (ug/L)` = Cc[which.max(tad)] * 1000,
      `Tmax (h)` = tad[which.max(Cc)] * 24,
      `Time > 0.02 mg/L (days)` = sum(Cc > 0.02) * 0.25,
      `Time > 0.01 mg/L (days)` = sum(Cc > 0.01) * 0.25,
      .groups = "drop"
    )
}

single_dose <- typical_cycle(1L)
#> ℹ omega/sigma items treated as zero: 'etalthalf_abs1', 'etalthalf_abs2', 'etalvc', 'etaiov_lthalf_abs2_1'
#> Warning: multi-subject simulation without without 'omega'
steady_state <- typical_cycle(18L)
#> ℹ omega/sigma items treated as zero: 'etalthalf_abs1', 'etalthalf_abs2', 'etalvc', 'etaiov_lthalf_abs2_1'
#> Warning: multi-subject simulation without without 'omega'

published <- data.frame(
  arm = c("BMI < 25", "BMI >= 25"),
  `Cmax (ug/L)` = c(34.8, 19.8),
  `Cmin (ug/L)` = c(5.64, 7.15),
  `Tmax (h)` = c(45.6, 43.0),
  `Time > 0.02 mg/L (days)` = c(9.75, 0),
  `Time > 0.01 mg/L (days)` = c(19.0, 18.5),
  check.names = FALSE
)

bind_rows(
  published |> mutate(source = "Hand 2019 Table 2 (median)"),
  cycle_metrics(single_dose) |> mutate(source = "Model, single dose"),
  cycle_metrics(steady_state) |> mutate(source = "Model, 18th monthly dose")
) |>
  relocate(source) |>
  arrange(arm, source) |>
  dplyr::rename("Source" = source, "Stratum" = arm) |>
  knitr::kable(digits = 2)
Source Stratum Cmax (ug/L) Cmin (ug/L) Tmax (h) Time > 0.02 mg/L (days) Time > 0.01 mg/L (days)
Hand 2019 Table 2 (median) BMI < 25 34.80 5.64 45.6 9.75 19.00
Model, 18th monthly dose BMI < 25 41.85 5.80 48.0 11.75 21.00
Model, single dose BMI < 25 36.89 5.14 48.0 10.25 19.25
Hand 2019 Table 2 (median) BMI >= 25 19.80 7.15 43.0 0.00 18.50
Model, 18th monthly dose BMI >= 25 30.73 10.73 54.0 12.75 28.25
Model, single dose BMI >= 25 20.99 7.41 60.0 2.50 20.50

The single-dose column reproduces every Table 2 entry to within about 11% except Tmax, which is flat-topped and grid-limited. The eighteenth-dose column does not: it overshoots the higher-BMI Cmax by 55% and the higher-BMI Cmin by 50%, while overshooting the lower-BMI stratum by far less. The discriminating quantity is the ratio between the strata, which depends only on the covariate effect and not at all on the fat-free-mass anchor:

ratio_of <- function(tbl) {
  tbl <- as.data.frame(tbl)
  tbl[["Cmax (ug/L)"]][tbl$arm == "BMI >= 25"] /
    tbl[["Cmax (ug/L)"]][tbl$arm == "BMI < 25"]
}
cmax_ratio <- data.frame(
  source = c(
    "Hand 2019 Table 2",
    "Model, single dose",
    "Model, 18th monthly dose"
  ),
  cmax_ratio_high_over_low = c(
    ratio_of(published),
    ratio_of(cycle_metrics(single_dose)),
    ratio_of(cycle_metrics(steady_state))
  )
)
cmax_ratio$pct_diff_vs_published <- 100 *
  (cmax_ratio$cmax_ratio_high_over_low -
    cmax_ratio$cmax_ratio_high_over_low[1]) /
  cmax_ratio$cmax_ratio_high_over_low[1]
cmax_ratio |>
  dplyr::rename(
    "Source" = source,
    "Cmax ratio (high / low BMI)" = cmax_ratio_high_over_low,
    "% diff vs published" = pct_diff_vs_published
  ) |>
  knitr::kable(digits = 3)
Source Cmax ratio (high / low BMI) % diff vs published
Hand 2019 Table 2 0.569 0.000
Model, single dose 0.569 0.000
Model, 18th monthly dose 0.734 29.057 
stopifnot(
  # The single-dose Cmax ratio reproduces the published ratio; the
  # steady-state one does not. This is a pure covariate-effect quantity --
  # deterministic, with no cohort sampling and no dependence on the fat-free-
  # mass anchor -- so it is asserted tightly.
  abs(cmax_ratio$pct_diff_vs_published[2]) < 2,
  abs(cmax_ratio$pct_diff_vs_published[3]) > 20
)

The per-participant metrics in Table 2 therefore describe a single 900 mg injection cycle with no carry-over from preceding injections, and the rest of this vignette simulates that interval. The accumulation the model does predict across repeated monthly dosing is reported at the end for completeness.

Simulation

Each participant receives one 900 mg injection and is observed every 6 h over the following 28 days – the sampling grid Hand 2019 used for its own simulations (Methods, “Pharmacokinetic modelling and simulations”: “A plasma benzylpenicillin concentration was simulated every 6 h, between doses of a 28 day dosing period”).

ev <- rxode2::et(amt = dose_mg, time = 0, cmt = "depot", id = cohort$id) |>
  # Observations go on the ODE state `central`; rxode2 returns the algebraic
  # observable Cc as a column at those records.
  rxode2::et(seq(0, tau, by = 0.25), cmt = "central")

ev_df <- as.data.frame(ev) |>
  left_join(cohort, by = "id")

sim_raw <- as.data.frame(
  rxode2::rxSolve(mod, ev_df, returnType = "data.frame")
)
# rxode2 echoes the input covariate columns back in the solve; drop them before
# joining so the cohort frame stays the single authoritative source and dplyr
# does not create BMI.x / BMI.y suffixed duplicates.
sim <- sim_raw[, setdiff(names(sim_raw), setdiff(names(cohort), "id"))] |>
  left_join(cohort, by = "id") |>
  dplyr::rename(tad = time)

Replicating the published figures

Figure 3 – concentration-time profile by BMI stratum

Figure 3 of the paper is a prediction-corrected VPC of plasma benzylpenicillin over the 28 days from the last dose, stratified by BMI, on a log10 scale with the BLQ fraction marked. The equivalent simulated percentile bands are shown below.

bands <- sim |>
  filter(tad > 0) |>
  group_by(arm, tad) |>
  summarise(
    p10 = stats::quantile(Cc, 0.10),
    p50 = stats::quantile(Cc, 0.50),
    p90 = stats::quantile(Cc, 0.90),
    .groups = "drop"
  )

ggplot(bands, aes(x = tad)) +
  geom_ribbon(aes(ymin = p10, ymax = p90), alpha = 0.25, fill = "steelblue") +
  geom_line(aes(y = p50), linewidth = 0.8) +
  geom_hline(yintercept = 0.02, linetype = "dashed", colour = "firebrick") +
  geom_hline(yintercept = 0.0025, linetype = "dotted") +
  scale_y_log10() +
  facet_wrap(~arm) +
  labs(
    x = "Time from last dose (days)",
    y = "Benzylpenicillin (mg/L, log10 scale)"
  ) +
  theme_bw()
Simulated benzylpenicillin concentrations over a 28-day injection cycle, by BMI stratum. Replicates the layout of Figure 3 of Hand 2019. The dashed line is the 0.02 mg/L target; the dotted line is the assay LLOQ of 0.0025 mg/L.

Simulated benzylpenicillin concentrations over a 28-day injection cycle, by BMI stratum. Replicates the layout of Figure 3 of Hand 2019. The dashed line is the 0.02 mg/L target; the dotted line is the assay LLOQ of 0.0025 mg/L.

The paper’s Figure 3 reports the fraction of observations below the limit of quantification alongside the concentration percentiles. The simulated equivalent is shown below; the assay LLOQ was 0.0025 mg/L.

sim |>
  filter(tad > 0) |>
  group_by(arm, tad) |>
  summarise(frac_blq = mean(Cc < 0.0025), .groups = "drop") |>
  ggplot(aes(x = tad, y = frac_blq, colour = arm)) +
  geom_line(linewidth = 0.8) +
  labs(x = "Time from last dose (days)", y = "Fraction BLQ", colour = "Stratum") +
  theme_bw()
Simulated fraction of concentrations below the 0.0025 mg/L assay limit of quantification, by BMI stratum. Companion to the lower panels of Figure 3 of Hand 2019.

Simulated fraction of concentrations below the 0.0025 mg/L assay limit of quantification, by BMI stratum. Companion to the lower panels of Figure 3 of Hand 2019.

Figure 4 – exposure versus weight

Figure 4 plots the percentage of the dosing interval above 0.02 mg/L, peak concentrations and trough concentrations against body weight, for each BMI stratum. Hand 2019 built its weight axis from WHO and CDC growth references; the cohort here is sampled from the Table 1 demographics instead, so the comparison is of shape and separation rather than of the exact weight grid.

per_subject <- sim |>
  group_by(id, arm, WT, BMI, FFM) |>
  summarise(
    cmax = max(Cc),
    tmax = tad[which.max(Cc)],
    # Under a single-dose 28-day interval the paper's Cmin -- the trough
    # immediately before the next injection -- is the last concentration.
    cmin = Cc[which.max(tad)],
    # Time above target, from the 6-hourly grid the paper itself simulated on.
    days_above_target = sum(Cc > 0.02) * 0.25,
    pct_above_target = 100 * mean(Cc > 0.02),
    days_above_low = sum(Cc > 0.01) * 0.25,
    pct_above_low = 100 * mean(Cc > 0.01),
    .groups = "drop"
  )
per_subject |>
  select(arm, WT, pct_above_target, cmax, cmin) |>
  tidyr::pivot_longer(
    c(pct_above_target, cmax, cmin),
    names_to = "metric",
    values_to = "value"
  ) |>
  mutate(
    metric = factor(
      metric,
      levels = c("pct_above_target", "cmax", "cmin"),
      labels = c(
        "(a) % of interval > 0.02 mg/L",
        "(b) Cmax (mg/L)",
        "(c) Cmin (mg/L)"
      )
    )
  ) |>
  ggplot(aes(x = WT, y = value, colour = arm)) +
  geom_point(alpha = 0.35, size = 0.9) +
  geom_smooth(method = "loess", formula = y ~ x, se = FALSE, linewidth = 0.9) +
  facet_wrap(~metric, ncol = 1, scales = "free_y") +
  labs(x = "Weight (kg)", y = NULL, colour = "Stratum") +
  theme_bw()
Simulated exposure metrics over one injection cycle versus body weight, by BMI stratum. Replicates Figure 4 of Hand 2019 (a: percentage of the interval above 0.02 mg/L; b: peak concentration; c: trough concentration).

Simulated exposure metrics over one injection cycle versus body weight, by BMI stratum. Replicates Figure 4 of Hand 2019 (a: percentage of the interval above 0.02 mg/L; b: peak concentration; c: trough concentration).

Both of the paper’s qualitative conclusions reproduce: exposure falls monotonically as weight rises (because the 900 mg dose is flat above 20 kg, so the mg/kg dose falls), and the higher-BMI stratum has markedly lower peaks and a smaller fraction of the interval above target.

per_subject |>
  group_by(arm) |>
  summarise(
    `Time > 0.02 mg/L (days)` = stats::median(days_above_target),
    `Time > 0.02 mg/L (%)` = stats::median(pct_above_target),
    `Time > 0.01 mg/L (days)` = stats::median(days_above_low),
    `Time > 0.01 mg/L (%)` = stats::median(pct_above_low),
    .groups = "drop"
  ) |>
  dplyr::rename("Stratum" = arm) |>
  knitr::kable(
    digits = 1,
    caption = "Simulated median target attainment; compare Table 2 of Hand 2019 (9.75 days / 35% and 19.0 days / 65% for BMI < 25; 0 days / 0% and 18.5 days / 69% for BMI >= 25)."
  )
Simulated median target attainment; compare Table 2 of Hand 2019 (9.75 days / 35% and 19.0 days / 65% for BMI < 25; 0 days / 0% and 18.5 days / 69% for BMI >= 25).
Stratum Time > 0.02 mg/L (days) Time > 0.02 mg/L (%) Time > 0.01 mg/L (days) Time > 0.01 mg/L (%)
BMI < 25 9.5 33.6 19.1 67.7
BMI >= 25 0.0 0.0 15.5 54.9 
med_target <- per_subject |>
  group_by(arm) |>
  summarise(
    days02 = stats::median(days_above_target),
    days01 = stats::median(days_above_low),
    .groups = "drop"
  )

stopifnot(
  # The paper's headline finding: no child spends the whole interval above the
  # 0.02 mg/L target, and the higher-BMI stratum spends less of it above target
  # than the lower-BMI stratum. Asserted on stratum medians, never on an
  # extreme of a random cohort.
  all(med_target$days02 < tau),
  med_target$days02[med_target$arm == "BMI >= 25"] <
    med_target$days02[med_target$arm == "BMI < 25"],
  # Against the lower 0.01 mg/L threshold the paper reports the two strata as
  # nearly equal (19.0 vs 18.5 days) -- the flatter higher-BMI profile loses at
  # the high threshold but not at the low one. That inversion of the gap is a
  # structural consequence of the covariate and is reproduced here.
  abs(diff(med_target$days01)) < abs(diff(med_target$days02))
)

PKNCA validation

Non-compartmental analysis of the 28-day injection cycle gives the exposure metrics the paper tabulates. half.life is the flip-flop terminal half-life (that is, t1/2,abs-2), and clast.obs at 28 days is the paper’s Cmin, the trough immediately before the next injection.

# Only `!is.na(Cc)`: adding `time > 0` or `Cc > 0` would drop the time-zero row
# PKNCA needs to anchor the AUC interval.
sim_nca <- sim |>
  filter(!is.na(Cc)) |>
  select(id, arm, time = tad, Cc)

# Guarantee a time = 0 record per (id, arm) even if the grid ever changes.
sim_nca <- bind_rows(
  sim_nca,
  sim_nca |> distinct(id, arm) |> mutate(time = 0, Cc = 0)
) |>
  distinct(id, arm, time, .keep_all = TRUE) |>
  arrange(id, arm, time)

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

dose_df <- cohort |>
  mutate(time = 0, amt = dose_mg) |>
  select(id, arm, time, amt)
dose_obj <- PKNCA::PKNCAdose(dose_df, amt ~ time | arm + id)

intervals <- data.frame(
  start = 0,
  end = tau,
  cmax = TRUE,
  tmax = TRUE,
  clast.obs = TRUE,
  auclast = TRUE,
  aucinf.obs = TRUE,
  half.life = TRUE
)

nca_res <- PKNCA::pk.nca(
  PKNCA::PKNCAdata(conc_obj, dose_obj, intervals = intervals)
)
nca_long <- as.data.frame(nca_res) |>
  filter(
    PPTESTCD %in%
      c("cmax", "tmax", "clast.obs", "auclast", "aucinf.obs", "half.life")
  ) |>
  select(arm, id, PPTESTCD, PPORRES)

nca_long |>
  group_by(arm, PPTESTCD) |>
  summarise(
    median = stats::median(PPORRES, na.rm = TRUE),
    q25 = stats::quantile(PPORRES, 0.25, na.rm = TRUE),
    q75 = stats::quantile(PPORRES, 0.75, na.rm = TRUE),
    .groups = "drop"
  ) |>
  mutate(PPTESTCD = nlmixr2lib::ncaParamLabel(PPTESTCD)) |>
  dplyr::rename(
    "Stratum" = arm,
    "NCA parameter" = PPTESTCD,
    "Median" = median,
    "Q1" = q25,
    "Q3" = q75
  ) |>
  knitr::kable(digits = 4)
Stratum NCA parameter Median Q1 Q3
BMI < 25 AUC0-∞ (obs) 0.5948 0.4984 0.7258
BMI < 25 AUClast 0.4900 0.4058 0.5828
BMI < 25 Clast 0.0054 0.0024 0.0083
BMI < 25 Cmax 0.0372 0.0272 0.0482
BMI < 25 t½ 8.8906 5.7965 13.7716
BMI < 25 Tmax 2.0000 1.5000 2.5625
BMI >= 25 AUC0-∞ (obs) 0.4993 0.4083 0.5943
BMI >= 25 AUClast 0.3263 0.2625 0.3999
BMI >= 25 Clast 0.0061 0.0042 0.0072
BMI >= 25 Cmax 0.0191 0.0130 0.0251
BMI >= 25 t½ 16.2997 10.4121 25.4101
BMI >= 25 Tmax 2.5000 1.9375 3.2500

Comparison against the published NCA

Table 2 of Hand 2019 reports individual t1/2,abs-2, Cmin, Cmax and Tmax by BMI stratum as median (IQR) (range); those medians are the reference below. Because Table 2 summarises the paper’s 18 real participants at their own (unpublished) fat-free masses, the simulated side is the typical-value profile at the fat-free mass implied by the Table 1 medians, which is the like-for-like comparison: population typical value against observed median, with no dependence on how a virtual cohort’s weights were sampled.

Two unit conversions make the two sides commensurable: the paper’s concentrations are in ug/L and the model’s in mg/L (divide by 1000), and the paper’s Tmax is in hours while the model’s time unit is days (divide by 24). The paper’s t1/2,abs-2 is the reference terminal half-life, which is what it is under flip-flop kinetics.

typ_nca_in <- single_dose |>
  dplyr::rename(time = tad)

typ_conc <- PKNCA::PKNCAconc(typ_nca_in, Cc ~ time | arm)
typ_dose <- PKNCA::PKNCAdose(
  data.frame(arm = unique(single_dose$arm), time = 0, amt = dose_mg),
  amt ~ time | arm
)
typ_nca <- PKNCA::pk.nca(
  PKNCA::PKNCAdata(typ_conc, typ_dose, intervals = intervals)
)

typ_long <- as.data.frame(typ_nca) |>
  select(arm, PPTESTCD, PPORRES)
# Hand 2019 Table 2, median column, converted to the model's units.
reference_nca <- data.frame(
  arm = c("BMI < 25", "BMI >= 25"),
  cmax = c(34.8, 19.8) / 1000, # ug/L -> mg/L
  clast.obs = c(5.64, 7.15) / 1000, # the paper's Cmin, ug/L -> mg/L
  tmax = c(45.6, 43.0) / 24, # h -> day
  half.life = c(9.8, 20.3) # t1/2,abs-2 in days = the flip-flop terminal t1/2
)

nca_tbl <- nlmixr2lib::ncaComparisonTable(
  simulated = typ_long,
  reference = reference_nca,
  by = "arm",
  params = c("cmax", "clast.obs", "tmax", "half.life"),
  units = c(
    cmax = "mg/L",
    clast.obs = "mg/L",
    tmax = "day",
    half.life = "day"
  )
)
knitr::kable(nca_tbl, digits = 4)
NCA parameter arm Reference Simulated % diff
Cmax (mg/L) BMI < 25 0.0348 0.0369 +6.0%
Cmax (mg/L) BMI >= 25 0.0198 0.021 +6.0%
Tmax (day) BMI < 25 1.9 2 +5.3%
Tmax (day) BMI >= 25 1.79 2.5 +39.5%*
Clast (mg/L) BMI < 25 0.00564 0.00514 -8.8%
Clast (mg/L) BMI >= 25 0.00715 0.00741 +3.6%
t½ (day) BMI < 25 9.8 8.91 -9.1%
t½ (day) BMI >= 25 20.3 16.6 -18.2% 
attr(nca_tbl, "footnote")
#> [1] "* differs from reference by more than ±20%."
# ncaComparisonTable() returns the "% diff" column as formatted CHARACTER
# ("+6.0%", "-18.2%", with a trailing "*" on flagged rows), so strip the
# decoration before comparing.
as_pct <- function(x) as.numeric(gsub("[%*+ ]", "", x))
pct <- as_pct(nca_tbl[["% diff"]])
stopifnot(!anyNA(pct))
stopifnot(
  # A mis-transcribed dose, volume, half-life or unit conversion would move a
  # whole stratum by tens of percent. Asserted on the centre of the comparison
  # and on a robust envelope -- never on the single worst row.
  stats::median(abs(pct)) < 15,
  stats::quantile(abs(pct), 0.75) < 30
)

The one row that exceeds the 20% tolerance is Tmax in the higher-BMI stratum. That is expected and benign: the slow absorption stage produces a broad, nearly flat peak, so Tmax is poorly determined. On the 6-hourly grid the paper itself used, the higher-BMI peak spans several adjacent grid points of nearly identical concentration, and moving between them changes Tmax by hours while changing Cmax by well under 1%.

The simulated terminal half-life sits below the published post-hoc medians in both strata (8.9 versus 9.8 days, 16.6 versus 20.3 days) because the model’s typical values are being compared with the median of 18 participants’ empirical Bayes estimates, which is a different quantity and need not coincide.

The same comparison over the simulated cohort

For completeness, the cohort simulated above is compared against the same reference. This is a weaker check than the typical-value comparison because it additionally depends on how the virtual cohort’s weights and heights were sampled within each BMI stratum – something Hand 2019 does not report – and because the median of a lognormally distributed peak exceeds the peak of the typical-value profile.

cohort_tbl <- nlmixr2lib::ncaComparisonTable(
  simulated = nca_long,
  reference = reference_nca,
  by = "arm",
  params = c("cmax", "clast.obs", "tmax", "half.life"),
  units = c(
    cmax = "mg/L",
    clast.obs = "mg/L",
    tmax = "day",
    half.life = "day"
  )
)
knitr::kable(cohort_tbl, digits = 4)
NCA parameter arm Reference Simulated % diff
Cmax (mg/L) BMI < 25 0.0348 0.0372 +6.8%
Cmax (mg/L) BMI >= 25 0.0198 0.0191 -3.6%
Tmax (day) BMI < 25 1.9 2 +5.3%
Tmax (day) BMI >= 25 1.79 2.5 +39.5%*
Clast (mg/L) BMI < 25 0.00564 0.00537 -4.7%
Clast (mg/L) BMI >= 25 0.00715 0.00605 -15.4%
t½ (day) BMI < 25 9.8 8.89 -9.3%
t½ (day) BMI >= 25 20.3 16.3 -19.7% 
cohort_pct <- as_pct(cohort_tbl[["% diff"]])
stopifnot(!anyNA(cohort_pct))
stopifnot(
  # Loose by design; see the paragraph above. Still tight enough to catch a
  # factor-level transcription error, which would move a row by 100% or more.
  stats::median(abs(cohort_pct)) < 40,
  stats::quantile(abs(cohort_pct), 0.75) < 60
)

Accumulation across repeated monthly dosing

The model predicts real, stratum-dependent accumulation across repeated monthly injections even though Table 2 describes a single cycle. Because the slow absorption half-life is 8.88 days in the lower-BMI stratum and 16.56 days in the higher-BMI stratum against a 28-day interval, the higher-BMI stratum accumulates considerably more.

accum <- inner_join(
  cycle_metrics(single_dose) |> select(arm, cmax_sd = `Cmax (ug/L)`),
  cycle_metrics(steady_state) |> select(arm, cmax_ss = `Cmax (ug/L)`),
  by = "arm"
) |>
  mutate(
    accumulation_ratio = cmax_ss / cmax_sd,
    theoretical = 1 /
      (1 - 2^(-tau / c(8.88, 8.88 * 1.865)[match(arm, c("BMI < 25", "BMI >= 25"))]))
  )

accum |>
  dplyr::rename(
    "Stratum" = arm,
    "Cmax, single dose (ug/L)" = cmax_sd,
    "Cmax, 18th dose (ug/L)" = cmax_ss,
    "Accumulation ratio" = accumulation_ratio,
    "1 / (1 - 2^(-tau/t1/2,abs-2))" = theoretical
  ) |>
  knitr::kable(digits = 3)
Stratum Cmax, single dose (ug/L) Cmax, 18th dose (ug/L) Accumulation ratio 1 / (1 - 2^(-tau/t1/2,abs-2))
BMI < 25 36.886 41.845 1.134 1.127
BMI >= 25 20.987 30.726 1.464 1.449 

stopifnot(
  # The simulated accumulation ratio agrees with the closed-form ratio for a
  # first-order process whose rate constant is the slow absorption step --
  # another expression of the flip-flop claim. Deterministic, so tight.
  all(abs(accum$accumulation_ratio / accum$theoretical - 1) < 0.05)
)

Both strata accumulate, but the higher-BMI stratum roughly 1.5-fold against roughly 1.1-fold for the lower-BMI stratum, so repeated dosing narrows the gap between them that a single cycle shows. This is a prediction of the packaged model, not a published result of Hand 2019.

Assumptions and deviations

  • Fat-free mass is the size descriptor for kel as well as for V. The paper states only that “fat-free mass was the best size parameter for allometric scaling on V”. Applying it to kel too is forced by the paper’s own equivalence claim: kel was fixed “with allometric scaling with an exponential of -1/4 (equivalent to an exponential of 3/4 for CL and 1 for V)”, and because CL = kel . V, the exponents -1/4 and 1 sum to 3/4 only if both are scaled by the same descriptor. Scaling V by fat-free mass and kel by total body weight would contradict the stated equivalence.
  • The -1 correlation is encoded as -0.99. Table 3 fixes r(t1/2,abs-1, t1/2,abs-2) to exactly -1, which makes the 3 x 3 OMEGA block singular; rxode2’s Cholesky-based sampler cannot decompose it. Following the repository convention for published perfect correlations, every off-diagonal is scaled by 0.99 while the diagonal variances – the published IIV values – are left untouched. The resulting matrix is positive definite.
  • The third correlation is derived, not published. The paper reports a full covariance matrix but prints only two of the three correlations. With r(t1/2,abs-1, t1/2,abs-2) = -1 the two absorption etas are perfectly collinear, so r(t1/2,abs-1, V) = -r(t1/2,abs-2, V) = +0.746 is forced by arithmetic rather than assumed.
  • Inter-occasion variability is drawn once per solve. rxode2 parses but cannot simulate the native eta ~ var | occ multi-level syntax, so the 30% IOV on t1/2,abs-2 is encoded as a single occasion-indexed eta (etaiov_lthalf_abs2_1). The cohort simulation in this vignette is a single injection cycle, so that draw is exact; a user who simulates many cycles will reuse one draw across all of them and therefore under-represent within-subject cycle-to-cycle variation. The paper’s own reported IOV shrinkage is 46%, so this term is weakly informed in any case.
  • Table 2 is read as a single dosing cycle, not as a steady-state cycle. The paper does not state which interval its per-participant Cmin, Cmax, Tmax and time-above-target values cover. The model discriminates the two readings because the accumulation ratio differs sharply between the strata (about 1.1 in the lower-BMI stratum against about 1.5 in the higher-BMI one), and the published between-stratum Cmax ratio of 0.569 matches the single-dose prediction to three significant figures while the steady-state prediction is 29% away from it. The “Which dosing interval does Table 2 describe?” section above shows the arithmetic. Note that this reading applies only to Table 2; the Figure 4 dosing simulations are explicitly described in the Methods as being at steady state, and the “Accumulation across repeated monthly dosing” section reports what the packaged model predicts there.
  • The fat-free-mass equation constants are not in Hand 2019. The paper derives FFM from weight and BMI citing Anderson & Holford 2009 (reference 20) but prints none of that equation’s constants. The vignette uses the Janmahasatian / Anderson-Holford form already registered for the FFM column in inst/references/covariate-columns.md; the two are algebraically identical (WHSmax . HT^2 . WT / (WHS50 . HT^2 + WT) equals WHSmax . WT / (WHS50 + BMI)). The model file contains none of these constants – it consumes FFM as an input column – so this affects only the virtual cohort. The Anderson & Holford paediatric FFM model predates the 2015 Al-Sallami paediatric age correction, which is therefore not applied.
  • Sex enters only through the fat-free-mass derivation. Hand 2019 did not retain sex as a covariate; it is sampled here at the Table 1 ratio (44% male) solely because the FFM equation is sex-specific.
  • Body mass index is consumed as a continuous column. The model forms the >= 25 kg/m^2 indicator inside model() from a continuous BMI value, because the paper’s covariate is explicitly a thresholded continuous measurement and because BMI is also the second input to the fat-free-mass derivation.
  • The dose is the benzathine salt, and V is apparent. Doses are milligrams of benzathine benzylpenicillin G (900 mg = 1.2 MIU Bicillin L-A) while the measured analyte is benzylpenicillin, so V is a V/F that absorbs both the salt-to-penicillin mass conversion and the unknown bioavailable fraction. No separate bioavailability term is estimated, and absolute V should not be compared with intravenous benzylpenicillin volumes.
  • The lower 450 mg dose band is not exercised here. Australian guidelines give 450 mg (0.6 MIU) below 20 kg, and Hand 2019 simulated that band, but all 18 participants who contributed data weighed above 20 kg (Discussion, limitations: “all participants weighed > 20 kg, therefore simulation values < 20 kg were unable to be verified”). The cohort here stays inside the observed weight range.
  • Concentrations below the assay limit are not censored. The paper retained BLQ observations by the M3 likelihood method; the simulations here report continuous concentrations without applying the 0.0025 mg/L LLOQ, so simulated troughs are lower than any that could have been measured.
  • No supplement was used. Supplementary data for this article (Figure S1 and assay / NONMEM detail) is referenced by the paper but is not required for the model: Table 3 reports the complete final parameter set. No erratum for this article was located.