Population PK/PD analyses were conducted using pooled clinical data from a Phase I comparative PK study (SB16-1001) in healthy male subjects and a Phase III comparative efficacy and safety study (SB16-3001) in patients with PMO (Table 1). In the Phase I study, 168 male healthy volunteers (HV) were randomized to receive a single subcutaneous (SC) dose of 60 mg SB16, EU-DEN, or US-DEN (Lee et al., 2024). In the Phase III study, 456 PMO patients were randomized to receive either SB16 or reference EU-DEN 60 mg SC at month 0 and month 6. At month 12, patients initially treated with EU-DEN were randomized again to either continue on EU-DEN or were transitioned to SB16. Patients in the SB16 treatment group continued to receive SB16 at month 12 (Langdahl et al., 2024). The key demographic and clinical characteristics of the two study populations are summarized in Table 2. The difference in age, gender, height, and body weight among the subjects reflected the inherent distinctions between the study population (I.e., male HV vs. female PMO patients). Most subjects in both studies were Caucasian. These subject characteristics were considered during the model development process.

Both clinical trials were conducted in accordance with the Declaration of Helsinki (1996), the International Council for Harmonisation E6 (R2) Good Clinical Practice guidelines, and all applicable local regulatory requirements. Study protocols were reviewed and approved by the relevant Institutional Review Boards or Independent Ethics Committees at each study center. Written informed consent was obtained from all participants prior to enrollment.

Overall, 6,583 serum denosumab concentration data points from 615 subjects across the Phase I and III studies, collected at prespecified timepoints in each study, along with 1,716 lumbar spine bone mineral density (BMD) measurements from 456 patients in the Phase III study, were included in the for PK/PD modeling (Table 1). In addition, a comparison of subject characteristics between DEN and SB16 in the Phase III study (SB16-3001) is provided in Supplementary Table S1.

All samples were analyzed at certified laboratories using validated analytical methods. Serum denosumab concentrations were measured employing a validated electrochemiluminescence immunoassay (ECLIA) specifically developed for the detection and quantification of denosumab. The lower limit of quantification (LLOQ) for this assay was 20 ng/mL. Among these, 615 pre-dose samples (9.34%) with concentrations below the LLOQ were set to zero for model initialization.

Of the 4,262 post-dose serum concentration samples, 3,133 samples (73.51%) were quantifiable and included in the model estimation, while 1,129 samples (26.49%) were below the LLOQ and thus treated as missing values (MDV = 1). Although excluded from the likelihood calculation, they were retained in the dataset to preserve the time structure and dosing history, ensuring accurate model fitting and simulation.

BMD of the lumbar spine (L1-4) vertebrae were measured during screening and prior to dosing at months 6, 12, and 18 using dual-energy X-ray absorptiometry (DXA) scanners from either GE Lunar [GE Healthcare, Chicago, IL, United States] or Hologic [Marlborough, MA, United States]. All scanners were certified, and DXA scan images were analyzed by the central reading center Calyx [Massachusetts, MA, United States] accordance with established practice guidelines. All patients’ subsequent DXA scans were acquired on the same scanner on which they had their baseline scan in screening. Investigator sites conducted quality control activities for DXA systems prior to the first patient scan (EMA, 2025).

An overview of the full modeling and simulation workflow is provided in Figure 1. In summary, serum concentration data pooled from the Phase I and III studies were analyzed using a nonlinear mixed effects population modeling approach. Covariate screening was subsequently performed during PK modeling to identify significant predictors of PK parameters and to explore potential differences in PK profiles between the biosimilar, SB16 and the reference product, DEN.

PK/PD modeling was conducted in a sequential Population PK parameters and Data (PPP&D) method approach. The population PK parameters estimated from the PK model were fixed in the PK-PD model, and both PK and PD data were used simultaneously to estimate PD parameters and their IIV, as well as the individual PK parameters (Lunn et al., 2009; Lacroix et al., 2012; Zhang et al., 2003).

Covariate screening was also conducted during PD modeling to identify significant predictors of the PD outcome and to assess treatment effects (SB16 vs. DEN). In the final PK/PD model, treatment effects on clearance (CL) were intentionally included irrespective of statistical significance to enable comparative simulations between SB16 and DEN.

All PK/PD models were implemented using the stochastic approximation expectation-maximization (SAEM) method in Monolix Suite™ 2024R1 (Lixoft SAS, Antony, France), and data processing was performed using R software v4.3.1 (R Foundation for Statistical Computing, Vienna, Austria).

Consistent with the known PK characteristics of denosumab, exploratory analysis of serum concentration-time profiles following administration of SB16 and EU/US-DEN exhibited target-mediated drug disposition (TMDD) (Sutjandra et al., 2011; Gabrielsson et al., 2016). Accordingly, one and two-compartmental TMDD models with several approximations to drug-target receptor binding kinetics were tested (Dua et al., 2015; Gibiansky et al., 2008), and a two-compartment, quasi-steady-state (QSS) approximation of the TMDD model best described the PK data (Figure 2). The amount of drug administered subcutaneously ( A s c ) into the depot, the total (free and bound) concentration of the drug ( C t o t ) in the central compartment, the amount of free drug in the peripheral (A _p ) compartment, and total concentration of free and bound target RANKL ( R t o t ), as defined by the QSS model, are described by respective differential Equations 1–4: d A SC dt = − k a · A SC (1) d C tot dt = k a · A SC V C − CL + Q · C − K int · R tot · C K ss + C + Q   · A p V p (2) d A p dt = Q   · C − A p V p (3) d R tot dt = k syn − k deg · R tot − k int − k deg · R tot · C K SS + C (4)

In the absence of drug target level data, it was assumed that the free drug, target, and drug-target complex are in a QSS, where the binding rate is balanced by the sum of dissociation and internalization rates. This is described by steady-state constant K s s = k int + k o f f k o n (Dua et al., 2015; Gibiansky et al., 2008).

Under QSS conditions, the concentration of free drug (C) and drug-target complexes (RC) are described using Equation 5–6: RC = R tot · C K ss + C (5) C = 1 2 C tot − R tot − K ss + C tot − R tot − K ss 2 + 4 · K ss · C tot (6)

Prior to drug administration, it was assumed that the turnover rate of free target receptors ( R 0 ) was at steady-state, such that k s y n = R 0 × k deg (Equation 7), and no drug-bound target was present in the system at baseline. The baseline receptor concentration (R ₀ ), QSS constant (K _SS ), and internalization rate constant (k _int ) were estimated from the observed PK data under QSS assumptions. R 0 = k s y n k deg (7)

For the PD model, an indirect response model based on turnover kinetics was employed to describe the changes in BMD following drug administration. Based on the mechanism of action of denosumab which reduces bone resorption by preventing formation and activity of osteoclast the inhibition of the elimination rate of lumbar spine BMD ( k o u t was selected from among various indirect response models evaluated during PD model development (Dayneka et al., 1993). The inhibition effect was modeled using a sigmoid maximum effect (I _max) function, and the hill coefficient (hill) was estimated as a model parameter to characterize the steepness of the concentration-response relationship, as described in Equation 8 (Gabrielsson et al., 2007): dBMD dt = k in − k out · BMD ·   1 − I max · C hill IC 50 hill + C hill (8)

Prior to drug administration, it was assumed that the BMD turnover system was at steady state, such that the BMD production rate constant (kin) equals the product of the BMD elimination rate constant (kout) and the baseline BMD (BMD₀), as shown in Equation 9: k i n = k o u t   · B M D 0 (9)

Under this assumption, k _in was not directly estimated but calculated accordingly.

To confine the estimated I _max parameter within the biologically plausible range of 0 and 1, and to implement interindividual variability on I _max, a logit transformation was applied as shown in Equation 10 (Adedokun et al., 2020): I max = e I max F 1 + e I max F (10) where I_maxF is an unbounded parameter representing the logit-transformed value of I_max.

Interindividual variability (IIV) in most PK/PD parameters was modeled using an exponential distribution (Equation 11). The IIV for ith individual, η _i , is a random variable assumed to be independently selected from normal distribution with a mean of zero and variance of ω². P i = P T V × exp η i (11) where P T V : typical value (population value) of the fixed effect parameter

The residual error model describing the difference between observed and the model-predicted values was selected based on the likelihood ratio test (LRT) and inspection of goodness-of-fit (GOF) plots.

Based on the analysis, a combined error model (Equation 12) was selected to describe the residual error in serum denosumab concentration, and an additive error model (Equation 13) was selected for lumbar spine BMD. The residual variability (RV) for the jth observed value in the ith individual, ε_add,ij, ε_prop,ij, represent additive and proportional components, respectively. These errors are assumed to be independent and normal distribution, each with a mean of zero and a variance of σ². y i j = I P R E D i j × 1 + ε p r o p , i j + ε a d d , i j (12) y i j = I P R E D i j + ε a d d , i j (13) where y_ij: jth observed value in the ith individual; IPRED_ij: jth model predicted value in the ith individual; ε_add,ij additive error; ε_prop,ij proportional error.

Model discrimination and validation were performed using statistical and graphical approaches. For statistical methods, the LRT was used to differentiate between hierarchical models at p-value <0.05, considering that the difference of −2 log-likelihood (-2LL) of the model follows an approximate chi-square distribution. A p-value of 0.05, corresponding to a decrease in objective function value (OFV) of 3.84 points, was considered as statistically significant. For the Wald test, a 95% CI for each parameter estimate was constructed using the point estimate and standard error derived from the model outcomes and checked if 95% CI excluded 0.

Basic GOF plots and visual predicted checks (VPCs) were applied to compare the observed values and model-predicted values, and model prediction was evaluated in terms of central tendency and variability. VPCs were performed by simulating 1000 replicates, and the resulting simulations were overlaid and compared with the original PK/PD data in plots generated using Monolix.

Covariate selection was performed using the conditional sampling use for the stepwise approach based on the correlation (COSSAC) algorithm, as implemented in Monolix. Study population (HV vs. PMO patients), body weight, age, race, and treatment group (biosimilar SB16 vs. DEN) were tested on relevant PK/PD parameters. At each iteration for covariate screening, covariates were checked for addition at p-value of 0.01 corresponding to a 6.635 points difference in log-likelihood (LL) then removal at p-value of 0.001 corresponding to 10.828 points difference in LL. The procedure alternated between forward inclusion and backward elimination steps, and the final model was selected based on the LRT (Ayral et al., 2021).

Continuous covariates were incorporated into the model using a normalized power function centered at the median value, and categorical covariates were modeled using an exponential function as presented in Equations 14,15, respectively: T V P i = T V P r e f · C o v i C o v r e f θ (14) T V P i = T V P r e f · e θ · C a t i (15) where TVP_i is the typical value for ith individual, Cov_i is the ith individual covariate value, Cov_ref is the reference (median) value, θ is the estimated covariate effect, and Cat_i is a binary indicator variable representing the categorical covariate (1 for tested variable, 0 for reference).

The final covariate model was automatically selected by the COSSAC algorithm as the model with the lowest −2 log-likelihood among the candidate models explored through stepwise forward and backward procedures. After selection, the numerical stability of parameter estimates (e.g., standard errors and relative standard errors) was evaluated, and covariates associated with unstable estimates were removed.

Monte-Carlo simulations of serum denosumab concentration-time and lumbar spine BMD-time profile were performed by the treatment group and covariate subgroups using Simulx (Monolix Suite™ 2024R1). The simulations were based on the final population PK/PD models. For each of SB16 and EU/US-DEN treatment groups, simulation data were generated for 1,000 virtual participants receiving three successive SC administrations of 60 mg of drug at 6-month interval. The IIV and RV were incorporated to reflect population variability. Covariates were fixed based on predefined subgroup scenarios rather than sampled from distributions.

The fifth, 50th (median), and 95th percentiles of the PK metrics were calculated and compared across subgroups identified as significant predictors of PK during covariate assessment. Similarly, the percentage change from baseline in lumbar spine BMD was calculated and summarized using the median and 90% prediction interval (PI) for each subgroup.

The treatment effects of the biosimilar SB16 and reference DEN, pooled from EU/US-DEN data, on the CL and I_max were assessed during covariate analysis. Furthermore, empirical Bayes estimates (EBEs) for key PK and PD parameters were visually compared between the two treatments using boxplots.

For the simulations used to compare PK and PD profiles of SB16 and DEN, the treatment group was intentionally retained as a covariate in the final PK/PD model, irrespective of its statistical significance in the covariate analysis. For each subject, PK metrics such as C_max, time to peak drug concentration (T_max) and area under the concentration-time curve within a dosing interval (τ) at steady state (AUC_τ,ss) and change from baseline in lumbar spine BMD were calculated from the simulation data. The median and 90% PI were then summarized and compared between the two treatment groups.

The PK of denosumab was well-described by a two-compartment TMDD model with QSS approximation, adequately characterizing data from the Phase I study in healthy male volunteers and the Phase III study in female PMO patients (Table 3). IIV was estimated for all PK parameters. Statistically significant correlation between IIVs of CL/F and V_P/F was identified after testing the correlations among all the IIVs. The PD model was described using an indirect response model with inhibitory effect on the elimination process (Table 4). Including IIV on the baseline BMD (BMD₀) and the IC₅₀ significantly improved the model fit.

All the fixed-effect PK/PD parameters and random effect parameters for IIV were estimated with acceptable precision, as indicted by relative standard errors (RSEs) < 25%. Although IIV of inter-compartmental clearance (Q) in healthy subjects was estimated to be high (coefficient of variation, CV: 295.99%), the estimate of the IIV was also obtained with high precision (RSE: 6.97%), suggesting that the observed variability is supported by the data. No significant signs of parameter instability or over-parameterization were observed.

Covariate analysis revealed that the study population (HV vs. PMO patients) statistically significantly affected k_a (ratio of HV/PMO: 0.55), R₀ (ratio: 15.49), and Q (ratio: 0.18). Body weight influenced V_C/F, V_P/F, and CL/F through a power model, with exponents of 1.5, 0.93, and 0.52, respectively. Race influenced CL/F, with values 1.14 times higher in Black individuals and 1.22 times higher in Asians compared to those of Caucasians. No patient-related variables, including treatment groups, were identified as covariates for any PD parameter.

GOF plots for the PK/PD model demonstrated good agreement between observed and predicted values (Figures 3, 4). The individual prediction plots (IPRED) showed data points closely clustered around the line of identity, indicating that the model adequately captured IIV (Figures 3B, 4B). Residual diagnostics further supported model adequacy, as individual weighted residuals (IWRES) were randomly distributed around zero across both time and predicted concentration axes. No patterns or evidence of heteroscedasticity were observed (Figures 3C,D, 4C,D).

The VPC for the final PK/PD model showed that model predictions closely aligned with the observed data for serum denosumab concentrations and BMD (Figure 5). Most observed values fell within the 90% PI, indicating strong predictive performance.

Simulated denosumab exposure (steady-state AUC) across significant covariate subgroups are shown in Figure 6. The simulated AUC was approximately 4% lower in Phase III PMO patients compared to that in Phase I healthy subjects (Figure 6A). When stratified by body weight–represented by 45 kg (low), 64 kg (median), and 90 kg (high) the AUC was approximately 45% higher in the 45 kg group and 39% lower in the 90 kg group, relative to the 64 kg reference (Figure 6B). Additionally, simulated AUC values were approximately 11% and 19% lower in Black and Asian subjects, respectively, compared to Caucasians.

Model-based simulations were performed to evaluate how difference in denosumab exposure across covariates influences PD responses, measured by changes from baseline in lumbar spine BMD. The percentage change from baseline in lumbar spine BMD following 6-month treatment intervals over 18 months was slightly higher in the Phase I healthy subjects (6.65%) than in the Phase III PMO patients (6.45%) (Figure 7A). Simulated changes from baseline in lumbar spine BMD by weight groups were 7.11%, 6.66%, and 5.54% for the 45 kg, 64 kg and 90 kg groups, respectively (Figure 7B). Similarly, simulated changes from baseline in lumbar spine BMD were 6.45%, 6.15%, and 5.93% for Caucasian, Black, and Asian subjects, respectively (Figure 7C).

Treatment (SB16 vs. DEN) was not identified as a statistically significant covariate since decrease in OFV was observed during covariate screening. Consequently, it was not included in the final PK-PD model. Boxplots of the EBEs for key PK and PD parameters showed that the distributions of ETAs were highly comparable between the two treatment groups, with no significant difference observed (Figure 8).

For treatment-specific simulation, the treatment group was intentionally included as a covariate in the final PK-PD model. The implemented CL/F ratio of SB16 to DEN was 0.9982, indicating a negligible difference between treatments. As shown in Figure 9, the time–concentration profile, time-lumbar spine BMD profile, and percentage change from baseline in lumbar spine BMD over time were nearly identical between the biosimilar and the reference product. A comparison of secondary PK and PD parameters demonstrated highly similar results between the two treatment groups (Table 5). To further confirm the robustness of this comparability, an additional simulation was conducted using individual PK and PD parameters estimated specifically in postmenopausal women with osteoporosis. The resulting profiles remained consistent with the primary findings and are provided in Supplementary Figure S1.