As bone diseases continue to impact global health with serious effects on quality of life, mathematical modeling offers insights into the complex cellular and molecular mechanisms, which control bone remodeling. The objective of the current work is to develop a multi-scale computational model of bone remodeling which can predict bone diseases related to dysfunctional parathyroid gland activity.

The parathyroid gland is responsible for the production of parathyroid hormone (PTH), one of the major hormones in vertebrates besides calcitriol and calcitonin for regulation of calcium homeostasis and bone health (Hernández-Castellano et al., 2020). The pulsatile nature of PTH secretion is a fundamental characteristic shared by many hormones, where pulsatility is believed to modulate target organ responsiveness and shows deviations in disease states as well as circadian and seasonal fluctuations (Chiavistelli et al., 2015). The PTH secretion pattern is characterized by tonic (i.e., constant) and pulsatile components. In healthy subjects, the tonic part of PTH secretion constitutes the majority (70%), whereas approximately 30% is secreted in low-amplitude and high-frequency bursts occurring every 10–20 min, superimposed on the tonic secretion (Chiavistelli et al., 2015). While the exact mechanisms underlying pulsatile hormone secretion are complex and not yet fully understood, the biological importance of PTH pulsatility to bone metabolism is supported by experimental evidence showing that intermittent PTH drug administration produces anabolic effects, whereas continuous administration results in catabolic outcomes (Locklin et al., 2003; Capriani et al., 2012). To understand this complex regulatory mechanism, it is necessary to examine PTH effects at the cellular level.

PTH plays a central role in maintaining calcium homeostasis through a feedback loop: parathyroid cells express calcium-sensing receptors that detect changes in serum calcium levels, with low calcium triggering increased PTH pulse amplitude and frequency, while high calcium results in the opposite effect (Chiavistelli et al., 2015). This creates a regulatory cycle where low serum calcium stimulates PTH secretion, which enhances bone remodeling activity and calcium release, leading to increased serum calcium levels. As calcium levels rise, PTH secretion is reduced and the thyroid releases calcitonin, which lowers the blood calcium level (Sun et al., 2020). Beyond the parathyroid-bone axis, calcium homeostasis involves complex interactions between multiple organ systems, including bone, intra- and extracellular fluid compartments, gut (oral intake), kidney, and the parathyroid gland (Peterson and Riggs, 2010; Sun et al., 2020). Dynamic PTH secretion patterns are fundamental to calcium homeostasis regulation and the pathogenesis of bone diseases.

PTH targets the parathyroid hormone/parathyroid hormone-related protein receptor (PTH/PTHrP type 1 receptor), also commonly known as PTH1R (Cheloha et al., 2015). PTH1R is a G-protein-coupled receptor that regulates skeletal development, bone turnover and mineral homeostasis. PTH1R transduces stimuli from PTH and PTH-related protein (PTHrP) into the interior of target cells (i.e., cells of the osteoblastic lineage) to promote several divergent signaling cascades (Cheloha et al., 2015). Changes in the PTH secretion pattern have been associated with various diseases including primary and secondary osteoporosis, and hyperparathyroidism (Bilezikian et al., 2018; Schaefer, 2000; Harms et al., 1994b; a; Bonadonna et al., 2005), ultimately leading to an imbalanced bone remodeling activity and distorted calcium homeostasis. Simulating bone diseases related to PTH glandular secretion and the corresponding bone cellular activities is the first step in exploring efficient drug treatments.

Osteoporosis (OP) is one of the most frequent musculoskeletal diseases affecting people worldwide (Aibar-Almazán et al., 2022; Salari et al., 2021). OP is characterized by low bone mass and altered bone quality, which ultimately leads to bone fractures. This disease is characterized by imbalanced bone remodeling–the fundamental process that regulates bone homeostasis. In bone remodeling, osteoclastic cells resorb the existing bone matrix, while osteoblastic cells replace the bone matrix by initially forming osteoid, which subsequently gets mineralized. In OP the balance between bone resorption and formation is biased towards resorption with diminished bone formation. Depending on the underlying causes, several types can be distinguished including postmenopausal and senile osteoporosis (Aibar-Almazán et al., 2022). Over the last few decades, a large variety of drugs have been developed, which help combat osteoporosis. The rapid increase in bone biology knowledge has led to the development of mechanobiological pharmacokinetic-pharmacodynamic (PK-PD) models of osteoporosis treatments. As reviewed in Pivonka et al. (2024), these in silico models allow predictions beyond bone mineral density (BMD), i.e., bone microdamage and degree of mineralization. Hence, in silico trials may serve as complementary tools to experimental studies, potentially contributing to our understanding of drug dosing and combinational treatments, though extensive validation of bone remodeling models across multiple scales remains essential before any clinical application.

Two primary formulations exist for modeling bone cell population dynamics (Cook et al., 2024). The first approach by Komarova et al. (2003) uses power laws where exponent terms represent the accumulated effects of signaling molecules governing both self- and externally-regulated cellular pathways. This results in a relatively small parameter space but creates inherent limitations in interpretability and extensibility. This approach has been further developed for spatio-temporal dynamics (Ryser et al., 2009; 2010). The accumulation of signaling effects makes it difficult to isolate the contributions of individual biomolecular factors on specific cell types. The second approach by Lemaire et al. (2004) employs mass kinetics formulations with explicit Michaelis-Menten and Hill equations for enzyme and ligand binding kinetics. This enables direct identification of how specific signaling factors affect osteoblast and osteoclast concentrations, providing clear biomolecular targets for drug treatments and disease modeling. The explicit modeling of receptor-ligand interactions results in a larger parameter space, but enables studying hormone signaling patterns in healthy and disease states. Thus, we adopt the mass kinetics framework in the present study due to its explicit incorporation of PTH1R receptors expressed on osteoblastic cells and PTH-PTH1R binding kinetics.

Currently no bone disease progression models with links to PTH glandular secretion patterns exist. Endocrine diseases are typically defined by comparing serum levels of endocrine factors with the “normal (or reference) range”. This reference range is used to discern hyper- and hypofunction of respective glands. Dynamic, time-dependent diseases may evolve within the normal range and are characterized by increased or decreased secretory dynamics (Chiavistelli et al., 2015; Schmitt et al., 2005). PTH glandular secretion governs osteoblastic cellular responses in bone and a disturbed function of the parathyroid gland can lead to development of progressive bone diseases. Consequently, the existing knowledge of healthy and pathological PTH glandular secretion patterns need to be incorporated into bone cell population models to accurately describe disease progression.

Schappacher-Tilp et al. (2019) developed a comprehensive mathematical model, which includes the major adaptive mechanisms governing the production, secretion, and degradation of PTH in patients with chronic kidney disease on hemodialysis. This work aimed to investigate the efficacy of parathyroid drugs. The model focused on simulating hemodialysis patients with secondary hyperparathyroidism, but has the potential to be extended and applied to other diseases such as primary hyperparathyroidism or hypo- and hypercalcemia. The ionized calcium concentration–which regulates the parathyroid gland response via the calcium sensing receptor–was provided as an input parameter for this model. However, (ionized) calcium concentration depends directly on bone turnover. Consequently, consideration of the bone remodeling process is essential.

Most computational models of bone remodeling are formulated as bone cell population models and include the action of PTH on bone cells in a simplistic manner: constant concentrations of PTH in the central- and/or bone fluid compartment in combination with one-state receptor models (Lemaire et al., 2004; Pivonka et al., 2008; Pivonka and Komarova, 2010; Scheiner et al., 2013; Trichilo et al., 2019). These models simulate catabolic action of PTH on osteoclasts by using a constant PTH concentration as input parameter to mimic both healthy state and particular bone diseases. PTH binds to its receptor PTH1R, expressed on osteoblasts (Datta and Abou-Samra, 2009); the receptor-ligand binding reaction is described by a one-state receptor model. The binding of PTH to its receptor is much faster than a cellular response such as differentiation, proliferation and/or production of ligands of osteoblasts, with binding occurring within minutes (Okazaki et al., 2008) compared to cellular processes that take tens of hours to days (Roberts et al., 1982). Hence, a steady-state assumption is used for solving these binding reactions resulting in a Hill-type equation for the activator/repressor functions, consistent with previous mathematical modeling approaches (Lemaire et al., 2004; Pivonka et al., 2008; Scheiner et al., 2013). The activator/repressor functions are based on receptor occupancy as a function of PTH concentration and total number of receptors expressed on osteoblasts (Datta and Abou-Samra, 2009). In the bone cell population models of Lemaire et al. (2004) and Pivonka et al. (2008), Pivonka and Komarova (2010), the PTH activator and repressor functions influence an intracellular communication pathway, which results in increased osteoclast activity and consequently catabolic bone resorption. While this type of approach is practical and simple to apply for creating catabolic bone remodeling states, it fails to address the link between pathological PTH release patterns and the different observed bone diseases.

Potter et al. (2005) proposed the first computational model to analyze PTH1R kinetics, focusing on the response to constant vs pulsatile dosing patterns of PTH. They introduced a measure of sensitization with values between 1 (highly sensitized) and 0 (desensitized). The study investigated clinically prescribed PTH drug dosing patterns and found a sensitization measure of around 0.9. However, they found a value of 0.89 for healthy glandular PTH secretion patterns. This proximity indicates that the proposed measure of sensitization is not meaningful to use as an activator function for osteoblast response as it is not able to distinguish between anabolic and catabolic actions of PTH.

Recently, Pivonka and co-workers applied the two-state receptor model to PTH1R to analyze the effects of PTH glandular and external dosing patterns on bone cell activity (Martonová et al., 2023). The work focused on clinically observed catabolic bone diseases related to perturbations of PTH glandular secretion. Following the approach proposed by Li and Goldbeter (1989), a cellular osteoblast activity function was developed, which can distinguish various aspects of the stimulation signal including peak dose, time of ligand exposure, and exposure period. Using this formulation, the potential of pharmacological manipulation of the diseased glandular secretion to restore healthy bone, cellular responsiveness via clinical approved external PTH injections was investigated. While it was mentioned that the so derived cellular activity function could be linked with a cell population model of bone remodeling, no description on how this could be accomplished was provided.

In this paper we develop a multi-scale bone cell population model based on a two-state receptor model of PTH1R accounting for dynamic PTH secretion patterns. This model is based on our previous work on a (osteoblastic) cellular responsiveness function, which can distinguish different PTH dosing patterns in health and disease (Martonová et al., 2023). We propose a calibration strategy that compares cellular activity values from the two-state receptor model with receptor occupancy values from the one-state model. We solve this problem in a semi-coupled way based on fulfillment of separation of time-scale condition: the two-state receptor model of PTH1R has a characteristic timescale ranging from tens of seconds to tens or hundreds of minutes, while the bone cell population model operates on timescales of hours to tens of days. After model calibration, we validate the model for several glandular disease states and analyze the effect of pulse characteristics on the bone cell response.

The bone cell population model (BCPM) describes the temporal behavior of osteoblasts and osteoclasts in various states including their regulation by receptor-ligand interactions. In this section, we present the bone cell population model proposed by Lemaire et al. (2004), and further developed and refined by Pivonka and co-workers (Pivonka et al., 2008; Pivonka and Komarova, 2010; Scheiner et al., 2013; Trichilo et al., 2019). These models provide a robust foundation for our study as they have demonstrated good qualitative agreement with experimental observations from the literature. The Lemaire model successfully reproduces known behaviors of the bone remodeling system, including tight coupling between osteoblasts and osteoclasts, the catabolic effect of continuously elevated PTH and RANKL, the reverse, anabolic effect of increased OPG, and metabolic bone diseases such as estrogen deficiency and glucocorticoid excess. Subsequent developments by Pivonka and co-workers have incorporated bone volume fraction evolution and identified physiologically sensible parameter combinations, while Scheiner’s extension coupled bone cell dynamics with mechanical feedback, reproducing key features of mechanoregulation including postmenopausal osteoporosis progression and mechanical disuse responses. Trichilo et al. further validated the model framework by applying it to ovariectomized rats and comparing bone volume fraction predictions with experimental data. Given our aim to study the effect of the pulsatile glandular secretion pattern of PTH on the bone cell response, here we emphasize the mechanisms through which PTH is integrated and operates in the model by Lemaire et al. A detailed description of the underlying dynamics can be found in the original publication (Lemaire et al., 2004).

The BCPM framework models bone cell populations as averaged concentrations representing the aggregate behavior of bone multicellular units (BMUs) at various stages of the remodeling cycle. While individual BMUs and the respective cells undergo periodic activation, resorption, and formation phases (Kenkre and Bassett, 2018), the model captures the smeared effect of many simultaneously active remodeling sites (Lemaire et al., 2004). In steady state, this approach yields constant average cell concentrations that reflect homeostatic balance rather than the oscillatory dynamics of individual BMUs.

All parameters used in the model including respective description and reference can be found in Supplementary Table S2. Figure 1 shows a schematic illustration of the semi-coupled bone cell population model with the two-state receptor model including intracellular communication pathways and respective parameters.

In this section, we describe the results of the calibration strategy given in Section 2.3, followed by the cellular activity and bone cell dynamics of the final, semi-coupled model. Finally, we demonstrate the effect of selected pulse characteristics on bone cell dynamics.

This study presents a novel approach to simulating bone disease progression related to the human parathyroid gland. We analyzed alterations of parathyroid hormone release patterns within the framework of bone cell population models (BCPMs) – a central step toward more physiologically realistic modeling of disease progression and accurate pharmacokinetics-pharmacodynamic (PK-PD) models (Pivonka et al., 2024). While traditional BCPMs treat hormone concentrations as constants, we demonstrated that dynamic hormonal patterns contain valuable information about disease states that is lost when complex temporal dynamics are simplified to averaged reference constants. Indeed, it is not clear how changes in dynamic hormonal release patterns could be implemented in the mechanistic BCPM framework. To address this question, we proposed a semi-coupled approach of a two-state PTH-PTH1R receptor model with the original bone cell population model of Lemaire et al. (2004) incorporating the bone volume fraction balance equation (Pivonka et al., 2008). Despite increased complexity, the cell concentration dynamics remain qualitatively similar to the original model formulation, demonstrating the robustness across different model variants (Figures 8, 9). Extension to more complex BCPMs (Pivonka et al., 2008; Pivonka and Komarova, 2010; Scheiner et al., 2013; Trichilo et al., 2019) is straightforward.

The model successfully reproduces key bone cell responses including increased osteoclast activity during disease states and time-delayed osteoblast response to catabolic osteoclast action. We showed that catabolic responses, traditionally obtained by increasing the PTH activator function ( π PTH ) , can now be achieved by replacing π PTH with the calibrated cellular responsiveness α R or integrated activity α T from the two-state PTH-PTH1R model. While the latter quantities were previously given theoretical bone-related interpretations in the two-state receptor model (Martonová et al., 2023), our study represents their first implementation within an actual bone cell population modeling framework.

Our analysis of single pulse characteristics (Section 2.4) shows that maximum cell response occurs at equal duration of on- and off-phase. This demonstrates the inherent trade-offs between pulse height and width: very long PTH pulses result in extended but low activity due to receptor desensitization, while short pulses produce intense but brief responses. The highest catabolic response–quantified by maximal integrated activity α T and cellular responsiveness α R – occurs when on- and off-phases are of equal duration, balancing pulse height and width while allowing for receptor resensitization between pulses.

We compared our novel approaches with the original BCPM formulation to evaluate the importance of including pulsatile PTH dynamics. We implemented disease states in the original model using a “scaling parameter” that reflects how diseases would be modeled as constant PTH elevations rather than capturing the full pulsatile patterns. The parameter is the ratio of maximum PTH concentration in disease and healthy state (Equation 14), and we multiplied it with the baseline PTH concentration of the BCPM (Equation 15). This approach does not require a complex two-state receptor model and could be used directly linked with existing BCPMs. While the novel approaches yield comparable results overall, the original model’s inability to incorporate pulsatile characteristics becomes evident in the case of glucocorticoid-induced osteoporosis, where the pulsatile pattern is characterized by increased relative amplitude without significant changes in maximum concentration (Figure 7). The original model, which neglects these pulsatile dynamics, fails to capture the catabolic cell responses and subsequent bone loss (Figures 9, 10). This addresses fundamental challenges in mathematical modeling of endocrine systems, where dynamic hormone secretion patterns are often oversimplified as constant values. The importance of capturing hormonal pulsatility is particularly evident in PTH signaling, where pulsatile versus continuous administration produces opposing effects on bone metabolism (Veldhuis, 2008). Our semi-coupled model incorporates both pulsatile characteristics and cellular desensitization through the two-state receptor model, addressing key research gaps identified in endocrine system modeling (Brucker-Davis et al., 2001; Zavala et al., 2019; Zavala, 2022).

The new model qualitatively predicts expected catabolic responses for the majority of PTH-driven bone diseases (hyperparathyroidism, hypocalcemia, glucocorticoid-induced osteoporosis), specifically, both the increased osteoclast activity and the resulting bone volume loss that characterize these conditions. It does not predict the bone loss encountered in osteoporosis and postmenopausal osteoporosis, where the model maintains homeostasis. The reason for this might be that these diseases are not exclusively linked to alterations in PTH release patterns, but also involve more significant pathophysiological changes such as estrogen depletion in PMO which directly regulates RANKL production by osteoblasts and osteocytes and/or TGF- β activity in old-age OP (Hawkins et al., 2021; Lu and Tian, 2023; Ruiz-Lozano et al., 2024). This is supported by the respective PTH characteristics (Table 1) that do not deviate significantly from the healthy pattern, confirming our model behaves as expected: predicting bone loss when PTH deviations are large enough compared to the healthy pattern while remaining stable when PTH patterns remain relatively normal despite the presence of other pathological mechanisms driving bone loss.

Validating model predictions against clinical data presents inherent challenges, as bone loss measurements typically compare disease states to controls rather than tracking progression from onset. For bone volume predictions (Figure 10), our model shows conservative estimates across conditions: For glucocorticoid-induced osteoporosis, we predict 0.7% trabecular bone loss in the first 60 days compared to clinical observations of approximately 5% loss (Boling, 2004). In primary hyperparathyroidism, our predicted loss of 0.8% after 60 days could reasonably accumulate to the observed 4%–5% difference between PHPT and control subjects (Christiansen et al., 1992), considering more rapid initial bone loss. The model’s highest catabolic response occurs in hypocalcemia, consistent with experimental data showing substantial bone loss (19% in rat models after 6 months (Jakubas-Przewłocka and Przewłocki, 2005), though our predicted magnitude is lower. These discrepancies likely reflect fundamental differences in disease simulation approaches (pulsatile human PTH patterns versus induced disease states), experimental duration, species transferability, and the instantaneous disease onset in our model versus progressive development in biological systems. While these comparisons suggest the need for parameter optimization based on expected bone loss, the heterogeneity and limited availability of consistent clinical data currently constrains such validation efforts, a challenge that has motivated various cross-methodological data-driven calibration approaches in bone modeling (Araujo et al., 2014; Baratchart et al., 2022).

Similar challenges apply to validating bone cell concentration predictions (Figures 8, 9), as osteoblasts and osteoclasts are rarely tracked over time in clinical settings. Comparison with murine studies helps establish reasonable numerical expectations despite inherent species differences. For hyperparathyroidism induced in mice, Siddiqui et al. (2017) demonstrated a 2.3-fold increase in osteoclast number per bone surface, while our model predictions range from 2.25 to 3.9-fold increases in OCa concentration across different approaches, showing good agreement with experimental data. This alignment is further supported by the 3.3-fold increase in TRAP5b–an enzyme produced by osteoclasts–observed in the same study. Our predicted decrease in bone volume fraction (0.2%–0.37%) is conservative compared to the experimental 9% reduction in the cortex, despite reasonable osteoclast predictions and only moderate underestimation of osteoblast activity (1.7–2.4-fold vs 3.67-fold experimental increase). Notably, no significant trabecular BV/TV loss was observed experimentally, which is consistent with our conservative bone volume estimates. For glucocorticoid-induced osteoporosis, Hofbauer et al. (2009) reported a 1.55-fold increase in osteoclast number per bone surface, closely matching our novel model formulations (1.62–1.67-fold increase of OCa-concentration), while the original approach showed insufficient response (1.06-fold increase). For osteoblast dynamics, our model simulates an initial decrease in active osteoblast concentration, which aligns qualitatively with the reduced serum osteocalcin levels reported for GIO-induced mice (Hofbauer et al., 2009). The initial decrease after disease onset is followed by a compensatory increase in response to elevated osteoclast activity. While this cellular response pattern is phenomenologically correct, we acknowledge that GIO involves complex mechanisms beyond PTH pulsatility alterations that our model does not capture, as our primary objective was demonstrating the importance of pulsatile hormone dynamics rather than comprehensive GIO modeling. Despite this cellular-level agreement, our model again predicted minimal bone volume changes over 80 days compared to substantial experimental bone area loss over 4 weeks (Hofbauer et al., 2009). Hofbauer et al. reported that trabecular BMD remained unchanged in GIO-induced mice, with most pronounced decreases in cortical and subcortical compartments. Likely due to limited sample sizes, BMD effects did not reach statistical significance. The agreement of our cellular predictions (osteoclasts and osteoblasts) with experimental data, coupled with underestimation of corresponding bone loss, suggests that recalibrating bone formation and resorption rate constants with appropriate human data could improve structural predictions. While these rates are currently assumed constant following established BCPM practice (Lemaire et al., 2004; Pivonka et al., 2008; Scheiner et al., 2013), experimental studies have demonstrated variability in individual osteoclast resorptive activity (Kanehisa and Heersche, 1988), and mathematical modeling of injury repair has provided in vivo evidence for time-variable cellular activity rates (Lo et al., 2021). Careful consideration is needed since our model tracks average cellular concentrations and thus inherently represents averaged resorption and formation rates over cell populations. Nevertheless, the order-of-magnitude agreement in cellular responses across both conditions provide confidence in the model’s mechanistic foundation.

We acknowledge several key limitations of our current approach. First, PTH levels are prescribed externally rather than evolving from endogenous physiological feedback mechanisms in both the original and our novel BCPM. PTH is either maintained at constant levels (original model) or follows prescribed pulsatile patterns (novel approach), without incorporating the calcium-PTH feedback loop that naturally regulates PTH secretion in vivo (Chiavistelli et al., 2015). This approach is appropriate for our primary objective of studying how specific hormone patterns affect bone cell dynamics and remodeling activity under controlled conditions. However, an autonomous model–incorporating calcium homeostasis and PTH regulation through calcium-sensing receptors–would facilitate analysis of the underlying causes of dysregulated hormone patterns themselves. An autonomous model would eliminate the current assumption of instantaneous disease state onset and allow gradual disease progression through feedback dysregulation.

Second, we model disease states as immediate transitions from a healthy PTH pattern (Figures 8, 9). This type of approach was original suggested by Lemaire et al. (2004) and subsequently improved by Pivonka et al. (2008) to account for temporal changing disease patterns. Our current approach could be extended towards gradual changing PTH patterns over given time intervals. Additionally, changes in bone volume fraction are constant after steady-states of cell concentration are reached. This implies constant bone loss independent of disease duration, which is not physiological. Mechanostat model incorporation could address this limitation (Scheiner et al., 2013; Martínez-Reina and Pivonka, 2019).

Third, our square-wave pulses for the PTH secretion pattern based on the original formulations (Li and Goldbeter, 1989; Martonová et al., 2023) represent an idealized version of hormone release patterns. A more physiologically realistic approach would model PTH degradation as an exponential decrease after the onset of each pulse, rather than an immediate switch to zero concentration, while maintaining total secreted PTH.

Finally, while we focused on PTH1R signaling in osteoblasts, recent studies have found that PTH1R is also expressed on osteocytes, where PTH directly induces upregulation of RANKL gene production. The resulting increased RANKL/OPG ratio leads to higher osteoclast recruitment and activation. This is consistent with the conservative bone loss predicted by the model. Additionally, PTH binding to PTH1R expressed on osteocytes downregulates sclerostin production, a formation inhibitor, which thus enhances bone formation (Ben-awadh et al., 2014; Marino and Bellido, 2024). These osteocyte-mediated effects represent additional pathways through which PTH influences bone remodeling beyond the osteoblast responses captured in our current model.

These limitations suggest several research directions. The framework could be extended to more sophisticated bone cell population models distinguishing between modeling and remodeling processes (Trichilo et al., 2019). The activity function α ( t ) could potentially distinguish between catabolic and anabolic pathways directly, making explicit, separate pathways unnecessary. Future developments could explore gradual concentration changes between healthy and disease states or directly use the time-dependent cellular activity function α ( t ) instead of constant PTH effect quantification.

Beyond PTH dynamics, our approach offers a template for other biological contexts where both temporal patterns and receptor adaptation require consideration. The methodology applies to mechanical loading patterns during habitual movement or exercise, where cells respond to pulsatile mechanical stimuli and adapt to sustained loads. This framework could also be adapted to other signaling pathways such as RANKL-RANK-OPG binding. This demonstrates the broader applicability in modeling various biological regulatory systems beyond simple Hill functions and constant stimuli. The interaction of bone remodeling, calcium homeostasis, and PTH secretion represents an interesting approach for future model development that could bridge the gap between prescribed hormone patterns and the physiological mechanisms that generate them.