Osteoarthritis (OA) is the most prevalent form of arthritis and a leading cause of musculoskeletal disability worldwide, affecting approximately 595 million people globally [1–3]. The knee is the most commonly affected joint, accounting for over 56% of cases [3, 4]. Current management remains primarily palliative, and no disease-modifying osteoarthritis drug (DMOAD) has received FDA or EMA approval [5]. This substantial unmet need, compounded by projections of 75%–95% increases in prevalence by 2050 [1, 3], underscores the urgency of developing mechanistic frameworks to guide therapeutic optimization.

Despite decades of research, multiple mechanistically promising DMOAD candidates have failed to achieve consistent structural modification in clinical trials [6–8]. Key explanations include OA heterogeneity across molecular endotypes [9–11], the insensitivity of radiographic joint-space width (JSW) as a structural endpoint [12, 13], and the historical treatment of OA as a homogeneous disease in clinical trial design [10, 11]. These challenges argue for mechanistic, systems-level approaches such as Quantitative Systems Pharmacology (QSP) to integrate biology, pharmacokinetics, and clinical endpoints.

Glucosamine, administered as glucosamine sulfate or glucosamine hydrochloride, is widely used for OA management [14, 15]. As an amino monosaccharide, it serves as an essential substrate for glycosaminoglycan (GAG) biosynthesis [16, 17]. Pharmacologically, glucosamine primarily exerts anti-catabolic effects by inhibiting IL-1β–induced inflammatory signaling, suppressing NF-κB activation, and reducing the expression of matrix metalloproteinases (MMP-1, MMP-3, MMP-13) and pro-inflammatory mediators (COX-2, PGE2, TNF-α) in human osteoarthritic chondrocytes [18–20]. Importantly, glucosamine does not reliably enhance anabolic matrix synthesis at therapeutic concentrations; comparative studies confirm that its chondroprotective effects derive primarily from limiting cytokine-driven matrix degradation rather than stimulating new matrix production [21–23].

Glucosamine exhibits rapid absorption, extensive first-pass metabolism, and low oral bioavailability (estimated at ∼19% from preclinical data; absolute human bioavailability has not been determined owing to the absence of an intravenous reference [24–26]), attributed primarily to gut-wall presystemic loss rather than hepatic extraction. Absorption is primarily mediated by facilitated transport through glucose transporters (GLUTs) rather than by passive diffusion, indicating that membrane permeability is not the rate-limiting step in oral absorption [24, 25]. Following a 1,500 mg once-daily dose, plasma concentrations increase ∼30-fold from baseline, reaching a Cmax of ∼10 μM at 3 h [26]. Glucosamine is subsequently eliminated in a multi-exponential fashion with substantial distribution into extravascular compartments; plasma concentrations remain measurable and above baseline for up to 48 h, and the terminal half-life was tentatively estimated at ∼15 h [26]. This sustained systemic exposure, together with synovial fluid concentrations that remain elevated at 12 h post-dose [27], supports once-daily dosing and the use of a steady-state synovial concentration as the pharmacodynamic driver (developed further in Module A: Pharmacokinetics). The glucosamine controversy has been characterized as fundamentally a pharmacokinetic issue [28], with preclinical studies demonstrating dose/concentration-dependent anti-inflammatory effects in adjuvant arthritis models [29, 30]. However, the assumption that crystalline glucosamine sulfate (cGS) has superior bioavailability compared with other glucosamine salt forms has recently been challenged. A randomized, double-blind, crossover pharmacokinetic study comparing cGS with regular glucosamine sulfate (rGS) in healthy volunteers found no significant difference in AUC_0–24, Cmax, or Tmax between formulations, with rGS showing numerically higher systemic exposure (AUC_0–24 GMR = 1.69, favoring rGS, p = 0.136) [31].

Furthermore, crystallographic analyses have demonstrated that commercially available “glucosamine sulfate” products, including stabilized “crystalline” forms, are physical mixtures of glucosamine chloride and potassium sulfate (or sodium chloride) rather than true sulfate salts [32]. Genuine glucosamine sulfate is highly hygroscopic and chemically unstable, precluding formulation into stable commercial dosage forms [32]. Since glucosamine sulfate dissociates in gastric fluid to yield the same glucosamine cation as glucosamine hydrochloride, the pharmacologically active moiety is identical regardless of the starting salt form [31, 32]. The USP reference standard for glucosamine is based on the hydrochloride salt, further underscoring that distinctions between commercial “sulfate” and “hydrochloride” products may reflect differences in excipients, manufacturing quality control, and labeling rather than fundamental differences in the active ingredient [33].

Despite these formulation controversies, glucosamine achieves therapeutically relevant synovial concentrations, with selective cartilage accumulation and prolonged synovial retention relative to plasma [27, 34–36]. This biphasic profile, with rapid initial plasma distribution but sustained synovial exposure, motivates the use of synovial fluid concentration as the driver of cartilage effects in the QSP model; the supporting data are detailed in the Methods (Module A: Pharmacokinetics).

The clinical efficacy of glucosamine remains contentious. Two pivotal European trials, both funded by the manufacturer of a patented crystalline glucosamine sulfate product, demonstrated significant structure-modifying effects over 3 years, with minimal joint space narrowing (−0.06 mm) compared with placebo (−0.19 to −0.31 mm), as well as symptomatic improvement [37, 38]. These structural benefits have not been replicated in independent trials. In contrast, the NIH-funded GAIT trial using glucosamine hydrochloride found no significant overall benefit over placebo, though patients with moderate-to-severe pain showed benefit with combination therapy (glucosamine hydrochloride 1,500 mg/day plus chondroitin sulfate 1,200 mg/day) [39]. Notably, 78% of GAIT participants had mild baseline pain, and the placebo response rate was 60.1%, substantially higher than typically observed in OA trials enrolling more symptomatic patients [39]. The controversy has historically been attributed to differences in formulation [28], but recent evidence challenges this interpretation. Pharmacokinetic comparisons have found no significant bioavailability advantage for crystalline glucosamine sulfate over regular glucosamine sulfate [31], and crystallographic evidence suggests both products deliver the same active moiety (glucosamine chloride) after gastric dissolution [32]. Therefore, the inconsistent trial results likely reflect a combination of factors: (1) differences in trial design, patient selection, and outcome assessment methodology; (2) the absence of dose-optimization studies, as every trial has used the same 1,500 mg dose established decades ago, without evidence of optimality [40]; (3) heterogeneity in patient populations and placebo response rates [41]; and (4) the reliance of positive structural outcomes on industry-sponsored trials that have not been independently replicated. A further consideration is that the glucosamine content of commercial formulations has been shown to vary substantially from label claims [42], and the trials reviewed here administered products based on labeled doses without independent verification of actual potency. The nominal doses reported across studies, including those tabulated in Table 1, therefore reflect label claims rather than confirmed delivered potency, introducing an additional source of variability that is independent of formulation bioavailability.

The glucosamine controversy presents an ideal case study demonstrating where traditional empirical approaches fall short and mechanistic modeling becomes essential. The disconnect between robust preclinical anti-inflammatory effects and inconsistent clinical outcomes suggests that simple dose-response relationships cannot capture the complexity of glucosamine’s action in OA. Although empirical pharmacokinetic and pharmacodynamic approaches are developed to adequately describe the data, they are generally not designed to provide quantitative insights into specific underlying mechanisms or to facilitate the use of these mechanistic insights to extrapolate to new conditions [43].

Key modeling implications include: (1) the need to integrate pharmacokinetics (absorption variability, synovial distribution) with pharmacodynamics (IL-1β/NF-κB pathway inhibition, cartilage turnover); (2) accounting for disease heterogeneity and patient stratification based on OA severity; (3) linking molecular-level mechanisms to clinically measurable endpoints such as WOMAC scores and joint space width; and (4) exploring dose optimization beyond the empirically established 1,500 mg daily regimen. A QSP framework can reconcile the apparent contradiction between mechanistic evidence and clinical trial results by explicitly modeling the conditions under which therapeutic concentrations are achieved at target tissues.

The disconnect between robust preclinical anti-inflammatory effects and inconsistent clinical outcomes suggests that empirical dose-response relationships cannot capture the complexity of glucosamine’s action. Quantitative Systems Pharmacology (QSP) integrates computational modeling across biological and temporal scales and species [44], making it uniquely suited to address this challenge. Glucosamine’s mechanism spans absorption kinetics, plasma-synovial fluid partitioning, intracellular uptake, and downstream NF-κB inhibition, affecting cartilage homeostasis; processes that require mechanistic integration rather than empirical fitting.

The regulatory landscape increasingly supports such approaches. The FDA Modernization Act 2.0 (December 2022) removed the blanket mandate for animal testing, allowing computational and other non-animal methods as alternatives in drug development [45]. QSP submissions to the FDA have increased steadily, exceeding 80 submissions in 2023 alone [46].

Despite extensive investigation, many disease-modifying osteoarthritis drug (DMOAD) programs have failed to demonstrate consistent clinical benefit [8, 47]. Increasingly, these failures are recognized not solely as biological shortcomings but as challenges in quantitative translation between drug exposure, tissue pharmacology, and clinical outcomes. Model-Informed Drug Development (MIDD) frameworks integrate pharmacokinetics, pharmacodynamics, and systems biology to reduce uncertainty in dose selection and clinical trial design [48, 49]. Within this context, QSP models provide a mechanistic framework that links drug exposure to disease progression, enabling hypothesis testing, evaluation of variability drivers, and simulation of clinical scenarios prior to costly trials [43, 50].

Despite extensive literature, no comprehensive QSP model integrates glucosamine pharmacokinetics, mechanisms, and clinical outcomes in OA. Such a model would: (1) quantitatively reconcile contradictory trial results through a mechanistic understanding of PK differences; (2) predict optimal dosing strategies and formulations; (3) identify patient populations most likely to benefit; and (4) exemplify model-informed drug development (MIDD) aligned with the FDA Modernization Act 2.0 principles.

Glucosamine is an ideal QSP candidate given extensive human trial data on 3-year structural outcomes; well-characterized pharmacokinetics, including synovial fluid concentrations; elucidated mechanisms (NF-κB inhibition); and validated endpoints (joint space width, WOMAC) with established MCIDs.

The objectives of this study were to: (1) develop a QSP model integrating glucosamine PK, cartilage homeostasis, and clinical outcomes; (2) calibrate against landmark European trial data; (3) validate using independent datasets; and (4) apply translational simulations exploring dose-response, bioavailability enhancement, treatment duration, and patient stratification by disease severity. Beyond characterizing glucosamine pharmacology, this work evaluates whether a mechanistic QSP framework can inform development-relevant questions, including exposure requirements for structural modification, variability in clinical response, and optimization of dosing and formulation strategies.

A structured data extraction template was developed in Microsoft Excel with separate worksheets for: study characteristics (design, population, intervention), WOMAC pain outcomes (mean, SD, change from baseline by timepoint), joint space width measurements (absolute values, change, measurement method), cartilage biomarkers (CTX-II, COMP concentrations), pharmacokinetic parameters (absorption rate, bioavailability, clearance), placebo response dynamics, and responder rates. The template included standardized fields for data source documentation (table number, figure, text page) and quality flags.

Data extraction was performed independently by two reviewers using the extracted PDF text files. Discrepancies were resolved by consensus or reference to original publications. For outcomes reported only in figures, data were extracted using WebPlotDigitizer with inter-rater reliability assessment. When standard deviations were not reported, they were estimated from standard errors, 95% confidence intervals, or imputed from studies with similar populations following Cochrane Handbook recommendations.

Quality assurance included: (1) verification of extracted values against original publications; (2) unit standardization (WOMAC scores normalized to 0–100 scale, JSW in mm); (3) documentation of scale transformations; (4) flagging of imputed or estimated values. All extraction discrepancies and corrections were documented in the verification log.

The quantitative systems pharmacology model integrates three interconnected modules representing glucosamine pharmacokinetics, cartilage matrix dynamics, and clinical outcomes. The model was implemented as a system of ordinary differential equations that explicitly represents the mechanistic pathway from drug exposure to clinical benefit. The complete set of model equations and parameter details is provided in the Supplementary Appendix.

Five parameters were estimated by weighted least-squares: cartilage turnover rates (k_syn, k_deg), glucosamine anti-catabolic effect (I_max,deg), and pain dynamics (Pain₀, P_max). Pharmacokinetic parameters (Fsulfate = 0.22, FHCl = 0.11, CL = 12 L/h, Ceff,ref = 0.30 μg/mL) and the potency anchor (IC50,deg = 3.0 μg/mL) were fixed rather than estimated. We note that absolute oral bioavailability and systemic clearance have not been directly determined for glucosamine in humans; Persiani et al. [26] characterized plasma pharmacokinetics but could not estimate absolute bioavailability (no intravenous arm), and the available F estimates derive from animal studies (≈3–6% in horses, 12% in dogs, and 21% in rats [24], using glucosamine hydrochloride). The values adopted here should therefore be regarded as plausible assumed inputs. Because F and CL enter the steady-state exposure only through the ratio F/CL, uncertainty in clearance is equivalent to a proportional change in bioavailability and is captured by the bioavailability sensitivity analysis (Supplementary Figure S4). We further note that the published steady-state plasma AUC over the dosing interval (∼14,560 ng·h/mL at 1,500 mg [26]) implies an apparent oral clearance (CL/F) of ∼100 L/h; the nominal CL and F adopted here are therefore illustrative, literature-informed inputs rather than independently fitted quantities. The target-site driver in the implemented model is fixed by the effective synovial reference (∼0.30 μg/mL, Module A: Pharmacokinetics), which is set conservatively relative to the measured synovial data [35] and does not depend on the individual values of CL or F; the associated exposure uncertainty is bracketed by Supplementary Figure S4.

Observation Model: An explicit statistical observation model was assumed for each endpoint. For JSW measurements, observations were assumed to follow a normal distribution around model predictions with endpoint-specific measurement error: J S W obs   ∼   N J S W pred , σ JSW 2 where σ_JSW = 0.08 mm, reflecting the standard error of JSW measurement reported in the calibration trials. For WOMAC pain scores: P a i n obs   ∼   N P a i n pred , σ Pain 2 where σ_Pain = 5.0 WOMAC points, based on typical within-group standard deviations in OA trials. These observation-error estimates were used to compute weighted residuals for diagnostic purposes and to scale the relative contributions of each endpoint to the objective function.

Each optimization enforced box constraints and the structural constraint k_deg > k_syn (Supplementary Table S6). The resulting calibrated values are specified in the MATLAB pipeline (R2023b) for all downstream simulation, figure generation, and uncertainty analysis; bootstrap confidence intervals and the wide-range re-calibration re-estimate the free parameters using fmincon (interior-point algorithm).

Objective Function: The objective function was the weighted sum of squared errors (WSSE), equivalent to maximum likelihood estimation under the Gaussian observation model: WSSE = ∑ i J S W obs , i − J S W pred , i σ JSW 2 + ∑ j P a i n obs , j − P a i n pred , j σ Pain 2

This formulation ensures that JSW and WOMAC endpoints contribute to measurement precision in proportion to their precision, avoiding dominance by either endpoint due to scale differences. In the implemented objective, the WOMAC pain residuals were additionally down-weighted relative to the structural endpoint to prioritize the joint-space-width data that drive the disease-modifying predictions. Calibration targets included JSW change at Years 1 and 3 from Reginster 2001 and Pavelka 2002 (n = 8 observations) and WOMAC pain at baseline, Weeks 4, 8, 16, and 24 from GAIT 2006 for both placebo and GH arms (n = 10 observations), yielding 18 observations for estimating 5 parameters.

Data Sparsity Considerations: With 18 observations and 5 estimated parameters, the effective degrees of freedom (df = 13) are modest, limiting the precision of the parameter estimates. This constraint motivated our reliance on literature-derived values for pharmacokinetic and potency parameters rather than attempting to estimate them from the available data. Uncertainty quantification through bootstrap resampling accounts for the influence of individual data points on parameter estimates, and the sensitivity analyses in the Results characterize the robustness of predictions to parameter uncertainty.

Parameter estimates were constrained to physiologically meaningful ranges (E_max ≤ 1, I_max ≤ 1, k_deg > k_syn). Two structural identifiability considerations merit discussion. First, for inhibitory E_max models with fixed potency (IC₅₀), only the product I_max × (C/IC₅₀) is fully identifiable from efficacy data; consequently, uncertainty in IC₅₀ propagates inversely to I_max estimates. Second, the cartilage turnover parameters k_syn and k_deg enter the model in a manner that renders only their net balance (k_net = k_deg − k_syn) identifiable from JSW progression data; the individual rates are structurally correlated. We therefore report both the individual parameter estimates (for mechanistic interpretability) and the derived net catabolic rate (for inferential validity). The constraint k_deg > k_syn was imposed based on the biological requirement that OA cartilage exhibits net matrix loss; this assumption is consistent with the progressive JSW narrowing observed in placebo arms of all calibration trials. Sensitivity analyses examining these relationships are presented in Results.

Model validation followed a multi-level strategy addressing structural, mechanistic, and predictive validity.

The validated model was used to conduct four translational simulation analyses to address clinically relevant questions regarding the optimization of glucosamine therapy. Each simulation used a virtual population of 500 subjects with interindividual variability sampled from expanded parameter distributions reflecting real-world population heterogeneity (coefficients of variation approximately 1.5× those used for calibration; see Supplementary Table S1).

To address recent evidence questioning the assumed bioavailability difference between glucosamine salt forms [31, 32], a post hoc sensitivity analysis was conducted comparing model predictions under four bioavailability scenarios: the base case (F_sulfate = 0.22, F_HCl = 0.11), equal bioavailability at the lower value (both F = 0.11), at the higher value (both F = 0.22), and at an intermediate value (both F = 0.165). For each scenario, the 3-year JSW treatment effect, 24-week WOMAC pain difference (GAIT design), and minimum effective dose were calculated to quantify the dependence of model predictions on the formulation-related bioavailability assumption.

All analyses were conducted in MATLAB R2023b (MathWorks, Natick, MA). The model implementation used forward Euler integration with weekly time steps for the numerical solution of the GAG dynamics. Primary parameter estimates were obtained by weighted least-squares optimization and are specified in the calibration script (M04_set_parameters.m). The bootstrap confidence intervals (M07_bootstrap_ci.m) and the wide-range IC_50,deg re-calibration (M15_IC50_recalibration.m) re-estimate the free parameters using fmincon (interior-point algorithm, Optimization Toolbox). Visualizations were created using MATLAB’s graphics functions.

The systematic review was conducted following PRISMA guidelines. The QSP model development followed good practices for model-informed drug development as recommended by the International Society of Pharmacometrics (ISoP) and the International Consortium for Innovation and Quality in Pharmaceutical Development (IQ) guidelines [59, 60].

The complete analysis pipeline, including all MATLAB scripts for model calibration, validation, and simulation, is publicly available at https://github.com/hamedgilzadkohan/glucosamine-qsp. The repository includes: (1) QSP model implementation (simulate_glucosamine.m); (2) the calibration script that specifies the parameter values, with documented bounds (M04_set_parameters.m); (3) bootstrap and profile likelihood analysis scripts (M07_bootstrap_ci.m, M06_profile_likelihood.m); (4) virtual population simulation code (M09_virtual_population.m); (5) translational simulation code (M11_translational_simulations.m); (6) visualization scripts for all figures (M12_figures.m, M13_supplementary_figures.m); (7) equal-bioavailability sensitivity analysis (M14_equal_bioavailability_sensitivity.m); and (8) a master script (run_all.m) that reproduces all analyses. A random seed of rng (42) was used for all stochastic analyses (bootstrap resampling, virtual population generation, and translational simulations) to ensure reproducibility. For computational efficiency, the steady-state synovial concentration was implemented using an equivalent reference-scaling approach, where the effective concentration at standard dosing (Fsulfate = 0.22, Dose = 1,500 mg) yields a synovial concentration of approximately 0.30 μg/mL, and other formulations and doses are scaled proportionally.

Local sensitivity analysis identified parameter subsets with differential influence on model outputs (Figure 3). For the primary efficacy endpoint (3-year JSW treatment effect), the most influential parameters were baseline GAG content (GAG₀), degradation rate (k_deg), and synthesis rate (k_syn), reflecting the dominant role of cartilage turnover balance. Among modifiable drug-related parameters, bioavailability (F_sulfate) and anti-catabolic potency (I_max,deg) were most influential, indicating that the magnitude of the treatment effect is primarily driven by drug exposure and its interaction with the degradation pathway.

For disease progression in the placebo arm (3-year JSW change), the dominant parameters were baseline GAG content (GAG₀) and the cartilage turnover rates (k_deg, k_syn), with sensitivity coefficients >2.0. This confirms that the balance between synthesis and degradation, independent of drug effects, governs the natural progression of disease.

For symptomatic outcomes (Week 24 placebo pain), baseline pain (Pain₀), GAG content, and placebo effect magnitude (P_max) were most influential. The placebo rate constant (k_pl) showed high sensitivity within the physiological range but exhibited diminishing influence at higher values, consistent with parameter sloppiness beyond k_pl ≈ 0.05 days^-1 (Supplementary Figure S5). Time-varying sensitivity analysis showed that parameter influences remained stable over the 3-year simulation period, supporting the model’s structural validity. The bioavailability analysis revealed a near-linear relationship between F and treatment effect, suggesting that formulation improvements could proportionally enhance efficacy (Supplementary Figure S4).

A virtual population of 500 patients was generated by sampling parameters from distributions reflecting interindividual variability (CVs of 15–50%, depending on parameter). The population simulations demonstrated robust model performance with appropriate uncertainty quantification (Figure 4).

For JSW outcomes over 3 years, the 90% prediction intervals captured the observed data from both calibration studies. The placebo arm showed a median JSW change of −0.446 mm (90% PI: −1.668 to +0.216 mm), while the glucosamine sulfate arm showed a median JSW change of −0.247 mm (90% PI: −1.494 to +0.492 mm). The median treatment effect was 0.186 mm (mean 0.208 mm), consistent with the calibration point estimate. Notably, the distribution of individual treatment effects under calibration conditions revealed substantial benefit across the virtual population. Under these conditions (crystalline glucosamine sulfate, calibration trial populations, CV 15–50%), 98.2% of virtual patients showed positive benefit (treatment effect >0), and the majority exceeded the MCID threshold of 0.1 mm JSW preservation. However, this high responder rate reflects the relatively low interindividual variability assumed for controlled-trial populations; the Translational Simulations in the Results, which employ expanded variability (CV 25–65%) to represent real-world heterogeneity, yield a more conservative 83% responder rate at the standard 1,500 mg/day dose. This difference reflects both the broader population variability (CV 25–65% vs. 15–50%) and the stricter MCID criterion (>0.1 mm) applied in the translational analysis versus the positive-benefit criterion (>0) reported here. The distribution was approximately normal. These high responder rates reflect the calibration scenario representing ideal conditions; the Translational Simulations in the Results explore responder rates under varied assumptions, including broader population heterogeneity.

External validation results are summarized in Table 6 and Figure 5. The GUIDE 2007 trial [51] provided primary placebo-controlled validation using the same crystalline glucosamine sulfate formulation as the calibration studies. The model accurately predicted trajectories for both the placebo and glucosamine sulfate arms over the 6-month treatment period (Figure 5A). All observed data points fell within the 90% prediction intervals; because the interval width (driven by the assumed interindividual variability) is wide, this coverage indicates consistency rather than a stringent test of accuracy. At Month 6, predicted pain scores were 24.9 and 23.5 for placebo and glucosamine sulfate arms, respectively, compared to observed values of 30.5 and 25.5 (Table 6).

The predicted treatment effect (glucosamine minus placebo) at 6 months was −1.4 points, whereas the observed effect was −5.0 points, resulting in an underestimation of the between-arm treatment difference by 3.6 points (72% relative error). This is distinct from the systematic over-prediction of pain reduction noted above, which affected both arms similarly. Mean absolute error (MAE) was 3.8 WOMAC points, and root mean square error (RMSE) was 4.2 points, both meeting the pre-specified “excellent” criterion of <5 points (Figure 5D; Table 6).

The validated model was applied to four translational simulation analyses. These simulations employed a broader virtual population (N = 500 per scenario) with expanded distributions of interindividual variability (coefficients of variation approximately 1.5× those used for calibration) to represent real-world heterogeneity beyond the controlled trial populations used for calibration. Consequently, predicted responder rates in these simulations are more conservative than those observed in the internal validation (83% vs. 98.2% at 1,500 mg/day), reflecting the expected attenuation when extrapolating from ideal trial conditions to diverse clinical populations.

To assess the dependence of model predictions on the assumed bioavailability difference between formulations, four scenarios were compared (Table 7). Equalizing bioavailability at the lower value (F = 0.11) reduced the predicted 3-year JSW treatment effect by 48% (from 0.211 to 0.109 mm) and shifted the minimum effective dose from 750 to 1,500 mg/day. At an intermediate value (F = 0.165), the treatment effect decreased by 24% (0.161 mm) with a minimum effective dose of 1,000 mg/day. Equalizing at the higher value (F = 0.22) preserved the JSW treatment effect but nearly doubled the predicted GAIT pain difference (1.0–1.9 WOMAC points). These results demonstrate that the model’s structural and symptomatic predictions are sensitive to the bioavailability assumptions.

The principal findings of the translational simulations and their clinical implications are summarized in Table 8.

This study presents a mechanistic quantitative systems pharmacology (QSP) framework for evaluating the disease-modifying potential of glucosamine sulfate in knee osteoarthritis. The model integrates pharmacokinetics, cartilage homeostasis dynamics, and clinical pain outcomes to provide a system-level understanding of glucosamine’s disease-modifying potential. The results provide quantitative evidence supporting several development-relevant conclusions. First, the model was calibrated to the landmark trial results from Reginster et al. [38] and Pavelka et al. [37], reproducing a 3-year joint space width (JSW) treatment effect of 0.211 mm against the observed 0.240 mm (12% error); this agreement reflects a constrained fit under the assumed potency anchor rather than an independent structural prediction, a point developed in “Re-calibration across the biologically realistic potency range” and Strengths and Limitations. Second, external validation against the independent GUIDE 2007 trial [51] demonstrated acceptable predictive performance for glucosamine sulfate monotherapy (MAE = 3.8 points on the WOMAC scale), with all observations falling within the model’s 90% prediction intervals and a treatment effect prediction error of 3.6 points. Exploratory comparison against the MOVES 2015 trial (combination therapy without placebo arm) [52] demonstrated that the glucosamine-only model appropriately under-predicted combination efficacy by approximately 15–20 WOMAC points, providing a quantitative estimate of chondroitin sulfate’s contribution that is consistent with its known pharmacology. Third, virtual population simulations under calibration conditions showed that 98.2% of patients benefited, with a median treatment effect of 0.186 mm; however, as detailed below and in “Re-calibration across the biologically realistic potency range”, this predicted benefit is strongly contingent on the assumed anti-catabolic potency and should not be interpreted as established efficacy. Fourth, translational simulations using expanded population heterogeneity identified actionable insights including a minimum effective dose of ∼750 mg/day, mean treatment effect of 0.22 mm at the standard 1,500 mg/day dose (with 83% responder rate), potential for 73% efficacy enhancement through bioavailability improvement to 40%, and consistent treatment effects (NNT = 2) across all disease severity stages from very early to severe OA.

From a drug development perspective, the present model demonstrates how mechanistic integration of pharmacokinetics and disease biology can clarify sources of variability that may obscure therapeutic effects in clinical trials. The simulations suggest that insufficient or heterogeneous joint exposure may contribute substantially to the historically inconsistent clinical outcomes reported for glucosamine. Clinically, this implies that variability in patient response may reflect pharmacokinetic limitations rather than the absence of biological activity, helping reconcile previously inconsistent clinical trial outcomes without requiring alternative biological mechanisms. Within a development program, such insights could support exposure-guided trial design, formulation optimization, and enrichment of patient populations predicted to achieve therapeutic tissue concentrations. Collectively, these findings illustrate how mechanistic QSP modeling can reduce uncertainty in efficacy interpretation and support quantitative decision-making regarding dose selection, formulation development, and clinical trial design. Such quantitative frameworks align with emerging model-informed drug development (MIDD) practices increasingly encouraged by regulatory agencies to contextualize efficacy outcomes and support dose justification [64]. It must be emphasized, however, that these interpretations are conditional on the model’s potency and exposure assumptions. As shown in “Re-calibration across the biologically realistic potency range”, when the anti-catabolic potency is set to values consistent with the concentrations at which glucosamine acts in vitro, the predicted structural benefit largely disappears at achievable synovial exposures, indicating that the historically inconsistent trial outcomes may reflect a genuine absence of sufficient target engagement rather than merely heterogeneous exposure.

It is important to note that treatment effect estimates differ across analysis contexts because of deliberate differences in assumptions about population variability. The calibration point estimate (0.211 mm) represents the expected effect under fixed parameter values that match the trial conditions. Internal validation with calibration-level variability (CV 15–50%) yielded a median effect of 0.186 mm (mean 0.208 mm). Translational simulations employed expanded variability (CV 25–65%) to reflect real-world heterogeneity beyond controlled-trial populations, yielding a mean treatment effect of 0.22 mm with a lower responder rate (83% vs. 98.2%). These differences are not contradictory; instead, they demonstrate how population heterogeneity attenuates responder rates while maintaining consistent underlying pharmacology and comparable mean treatment effects. All three estimates are consistent with the observed clinical trial treatment effect of 0.24 mm when appropriate uncertainty is accounted for.

This work exemplifies the transformative potential of computational approaches in pharmaceutical development, directly aligned with recent regulatory evolution toward New Approach Methodologies (NAMs) and Model-Informed Drug Development (MIDD). On December 29, 2022, President Biden signed the FDA Modernization Act 2.0 into law, fundamentally revising the Federal Food, Drug, and Cosmetics Act of 1938 by removing the blanket mandate for animal testing and explicitly authorizing “certain alternatives to animal testing, including cell-based assays and computer models” to support investigational new drug applications [45, 59]. This legislative milestone reflects growing recognition that over 90% of drugs that appear safe and effective in animal models fail to receive FDA approval in humans, primarily due to safety and efficacy concerns [60].

Our QSP modeling approach embodies the principles championed by this regulatory evolution. The ICH M15 draft guidance on General Principles for Model-Informed Drug Development, endorsed in November 2024, articulates that “appropriate use of MIDD can enable greater efficiency in drug development, while harmonized approaches to assessment can promote consistent and transparent evaluation of MIDD evidence to inform regulatory decision making” [48]. The MIDD Paired Meeting Program under PDUFA VII (2023–2027) prioritizes explicit dose selection, clinical trial simulation, and predictive mechanistic safety evaluation, all of which are addressed by our glucosamine QSP model [48, 49].

The FDA’s commitment to NAMs is further evidenced by the October 2024 Science Board recommendations to create a uniform qualification framework, establish centralized NAM databases, and develop transparent review processes for NAM-based applications [65]. Recent FDA/CDER/OND perspectives emphasize the agency’s “steadfast commitment to the 3Rs” (replace, reduce, refine) and encourage submission of NAM-based nonclinical data packages [66]. Our QSP model, calibrated to human clinical trial data without requiring additional animal studies, represents precisely the type of human-relevant in silico methodology envisioned by these initiatives.

QSP modeling has emerged as an increasingly accepted methodology in pharmaceutical development. A landscape analysis of FDA submissions through December 2020 demonstrated notable growth in QSP applications across drug development phases and therapeutic areas [71], with 2020 submissions alone accounting for approximately 4% of the agency’s annual IND submissions [72]. QSP approaches have been applied across therapeutic areas, including metabolism, autoimmunity, oncology, and neuroscience [44, 72]. Our glucosamine model extends this established methodology to osteoarthritis, a disease area with no approved DMOADs, where mechanistic modeling could inform both development strategy and regulatory decision-making.

The model architecture follows the six-stage workflow advocated by Gadkar et al. [50]: scope definition, data assembly, model specification, calibration, validation, and application. Our use of multi-start optimization, sensitivity analysis, virtual population generation, and external validation aligns with emerging best practices for QSP model assessment [71, 73]. The explicit representation of interindividual variability through parameter distributions (CV 15–50%) enables population-level predictions with appropriate uncertainty quantification, a critical requirement for regulatory applications. The virtual population framework implemented here represents an early disease-specific digital twin paradigm, in which mechanistic models simulate individualized disease trajectories under alternative therapeutic scenarios [74]. Such approaches enable in silico exploration of treatment strategies, allowing evaluation of dosing, formulation, and patient variability prior to clinical testing. By integrating mechanistic biology with patient-level variability, digital-twin–like models may provide a scalable framework for hypothesis testing and trial optimization in osteoarthritis drug development.

Strengths: This study has several notable strengths. The systematic literature review and data extraction pipeline ensured comprehensive capture of available evidence, progressing from 177 initial records through rigorous screening to 5 core calibration studies representing over 2,900 patients. The multi-module model structure, incorporating PK, cartilage homeostasis, and clinical outcomes, provides mechanistic interpretability beyond empirical curve-fitting. Calibration against both structural (JSW) and symptomatic (WOMAC pain) endpoints ensures relevance to the dual goals of disease modification and symptom relief. We used two entirely independent trials for external validation: GUIDE 2007 [51] for placebo-controlled validation, and MOVES 2015 [52] as an active comparator. The GUIDE trial was particularly valuable because it used the same crystalline glucosamine sulfate formulation as the calibration studies and included an accurate placebo arm, enabling direct validation of treatment-effect predictions. Finally, the translational simulation framework generates clinically actionable hypotheses regarding dose optimization, formulation enhancement, treatment duration, and patient selection.

Limitations: Several limitations merit consideration. First, the model was calibrated to a limited dataset of 18 aggregate observations (8 JSW, 10 WOMAC) from published trials, rather than individual patient-level data; this constrains parameter precision and precludes direct estimation of between-subject variability from the calibration data. Access to individual patient data from the calibration trials would enable more rigorous mixed-effects modeling and improved uncertainty quantification.

Second, the model was calibrated primarily to data from European trials using crystalline glucosamine sulfate; generalization to dietary supplement formulations with variable quality and bioavailability should be interpreted cautiously.

An important extension of this limitation concerns the validity of the calibration data themselves. The structural outcome calibration relies on the Reginster (2001) [38] and Pavelka (2002) [37] trials, both funded by the manufacturer of the patented crystalline glucosamine sulfate product tested. These positive structural results have not been replicated in independent, non-industry-funded trials. The model-predicted minimum effective dose of approximately 750 mg/day, the NNT estimates, and the bioavailability enhancement projections are all directly dependent on the effect sizes from the calibration trial. If the European trial results overestimate the true treatment effect, due to trial conduct factors, patient selection, or industry influence, all downstream model predictions would be proportionally optimistic. The apparent paradox that the model predicts efficacy at 750 mg/day while independent trials show minimal benefit at 1,500 mg/day likely reflects this calibration dependence rather than a genuine dose-response relationship. Relatedly, the external validation presented here did not incorporate independent trials that reported null results [75–77]. These trials were not suitable for formal treatment-effect validation against our calibration endpoints because they employed a discontinuation design [77], studied hip rather than knee osteoarthritis [76], or assessed structure by MRI using glucosamine hydrochloride rather than radiographic joint space width with crystalline sulfate [75]. Their existence nonetheless underscores that the present model is validated only against the subset of trials reporting positive structural outcomes, and that its predictions should be interpreted as conditional on those data. Recalibration against pooled data from both positive and null trials represents an important avenue for future work. Furthermore, recent pharmacokinetic evidence challenges the bioavailability assumptions underlying the model’s formulation comparisons. The base model assumes F = 0.22 for glucosamine sulfate and F = 0.11 for glucosamine hydrochloride, based on published estimates [26]. However, a randomized crossover study found no significant pharmacokinetic difference between crystalline and regular glucosamine sulfate [31], and crystallographic evidence demonstrates that commercially marketed “glucosamine sulfate” is chemically a mixture of glucosamine chloride and an alkali sulfate salt [32]. If the true bioavailability difference between formulations is negligible, the model’s explanation of the GAIT vs. European trial discrepancy on pharmacokinetic grounds alone would be insufficient, and alternative factors (trial design, patient selection, placebo response heterogeneity, industry sponsorship effects) would need to be invoked. Future model iterations should explore calibration to pooled data from both industry-sponsored and independent trials, and sensitivity analyses comparing model behavior under equal-bioavailability assumptions.

As shown in the Results (Table 7), the model’s key predictions depend substantially on the assumed two-fold bioavailability difference between formulations. When this difference is removed by equalizing bioavailability at the lower value, the predicted minimum effective dose shifts upward to 1,500 mg/day, the dose at which independent trials have failed to demonstrate significant structural benefit, and the pharmacokinetic explanation for the discrepancy between European and GAIT outcomes is correspondingly weakened. Because this assumption is challenged by recent evidence [31, 32], the formulation-based account of the trial discrepancy should be regarded as provisional.

Third, the pharmacokinetic module uses steady-state average synovial concentrations rather than dynamic plasma profiles, and therefore does not capture peak-to-trough fluctuations within the dosing interval. Two features mitigate this: glucosamine’s relatively long terminal half-life (∼15 h [26]) limits the magnitude of intra-interval fluctuation, and the cartilage turnover processes driving the modeled endpoints operate on a timescale of weeks to months, far longer than the dosing interval. This simplification is therefore unlikely to materially affect predictions of structural outcomes.

Fourth, although the model includes a minor direct analgesic component (β_drug = 5.0, contributing approximately 0.45 WOMAC points at therapeutic concentrations), glucosamine’s symptomatic effects are predominantly mediated through the structural preservation (GAG) pathway. Future work could explore whether direct anti-inflammatory mechanisms warrant a larger role in the pain model.

Fifth, the cartilage homeostasis module employs a simplified representation of GAG dynamics that does not capture the full complexity of matrix turnover, including collagen remodeling, chondrocyte senescence, and inflammatory cytokine cascades.

Sixth, while external validation against the GUIDE 2007 trial [51] was encouraging for glucosamine sulfate monotherapy, the MOVES 2015 trial [52] could not serve as a proper validation of treatment effect due to its use of combination therapy and absence of a placebo arm. Validation against additional independent cohorts, including those with different demographic characteristics, comorbidities, or concurrent medications, would strengthen confidence in the model’s robustness.

Seventh, the current model does not incorporate chondroitin sulfate pharmacology, limiting its applicability to predicting combination therapy.

Eighth, external validation revealed systematic over-prediction of pain reduction in the GUIDE 2007 trial, with the model under-predicting endpoint pain scores by 2–6 WOMAC points. This bias likely reflects heterogeneity in placebo-response dynamics between the calibration (GAIT: mild OA, 60% placebo response) and validation (GUIDE: moderate-severe OA) populations. The placebo response module, calibrated on 24-week data, may not reliably extrapolate to other patient populations or longer durations. Future iterations should incorporate population-specific placebo parameters or mixed-effects approaches to capture between-study variability.

Ninth, the model does not integrate structural heterogeneity within the joint (e.g., focal vs. diffuse cartilage loss) that may influence treatment response.

Tenth, although the core calibration studies were assessed as having low overall risk of bias (Supplementary Table S7), clinical trial data inevitably carry uncertainties in outcome measurement, dropout patterns, and adherence that propagate through calibration. The “some concerns” ratings for GAIT 2006 (high placebo response, differential dropout) and MOVES 2015 (open-label design) were taken into account when interpreting model performance on these datasets.

The anti-catabolic effect parameter (I_max,deg) was estimated at its physiological upper bound of 1.0. Several considerations inform the interpretation of this finding: i) the bootstrap 95% CI [0.6348, 1.0000] includes the point estimate, with 94.6% of replicates concentrated at the upper bound (≥0.95); the wide lower tail reflects bootstrap samples with less informative data configurations, which is expected behavior for boundary-constrained parameters. ii) sensitivity analyses demonstrated that model predictions are robust across the plausible I_max,deg range (0.8–1.0), with treatment effects of 0.18–0.22 mm all exceeding the MCID threshold (Supplementary Table S4). iii) only the composite quantity I_max,deg × (C_syn/IC_50,deg) is fully identifiable from clinical outcome data; the individual parameters are partially confounded. This limitation is inherent to calibrating mechanistic models against aggregate clinical endpoints, given the lack of direct measurements of target engagement. An important caveat tempers the interpretation of this near-complete inhibition. The in vitro studies demonstrating suppression of IL-1β–induced catabolic signaling did so at concentrations of approximately 10–1,000 μg/mL [53], which are one to three orders of magnitude higher than the steady-state synovial concentrations achievable after oral dosing (≈0.78 μg/mL; 4.34 μM [35]). The model reproduces the observed structural benefit only because the assumed potency anchor (IC50,deg = 3.0 μg/mL) places the achievable concentration on the steep part of the inhibition curve. When the model is re-calibrated using IC50 values consistent with the in vitro effective range, the predicted structural benefit falls below the minimal clinically important difference and the calibration fit degrades substantially (“Re-calibration across the biologically realistic potency range”). The predicted disease-modifying effect should therefore be interpreted as conditional on this optimistic potency assumption rather than as an established property of glucosamine.

Additionally, bootstrap analysis revealed a high correlation between the synthesis and degradation rate constants (k_syn and k_deg; r = 0.993), indicating that these parameters are not independently identifiable from the available data. This structural non-identifiability implies that, while the individual turnover rates carry substantial uncertainty (CV ≈ 17%), the net catabolic rate, knet = kdeg − ksyn, that drives disease progression is well-constrained (CV = 12.5%; Supplementary Table S5). The constraint k_deg > k_syn was imposed a priori based on the requirement that OA cartilage exhibits net matrix loss; relaxing this constraint would allow solutions with k_syn > k_deg that are biologically implausible for progressive OA. Future studies incorporating direct measurements of cartilage turnover biomarkers (e.g., CTX-II for degradation, CPII for synthesis) could help resolve this correlation by providing independent constraints on synthesis and degradation rates.