Abstract
Introduction:
Bone diseases significantly impact global health by compromising skeletal integrity and quality of life. In disease states linked to parathyroid hormone (PTH) glandular secretion, disrupted PTH patterns typically promote osteoclast proliferation, leading to increased bone resorption.
Methods:
While mathematical modeling has proven valuable in analyzing bone remodeling, current bone cell population models oversimplify PTH secretion by assuming constant levels, limiting their ability to represent disorders characterized by variations in PTH pulse characteristics. To address this, we present a novel semi-coupled approach integrating a two-state PTH receptor model with an established bone cell population model. Instead of conventional Hill-type functions, we implement a cellular activity function derived from the receptor model, incorporating pulsatile PTH patterns, cell dynamics, and intracellular communication pathways.
Results:
Our numerical simulations demonstrate the model’s capability to reproduce various catabolic bone diseases, providing realistic changes in bone volume fraction over a 1-year period. Notably, while direct implementation of PTH disease progression in the bone cell population model fails to capture diseases only characterized by altered pulse duration and baseline, such as glucocorticoid-induced osteoporosis, our semi-coupled approach successfully models these conditions.
Discussion:
This physiologically more realistic approach to endocrine disease modeling offers potential implications for optimizing therapeutic interventions and understanding disease progression mechanisms.
1 Introduction
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 (). 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 (). 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 (). 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 (; ). 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 (). 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 (). 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 (; ). 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 (). 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 (). Changes in the PTH secretion pattern have been associated with various diseases including primary and secondary osteoporosis, and hyperparathyroidism (; ; ; ; ), 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 (; ). 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 (). 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 , 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 (). The first approach by 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 (; ). The accumulation of signaling effects makes it difficult to isolate the contributions of individual biomolecular factors on specific cell types. The second approach by 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 (; ). 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.
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 (; ; ; ; ). 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 (); 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 () compared to cellular processes that take tens of hours to days (). 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 (; ; ). The activator/repressor functions are based on receptor occupancy as a function of PTH concentration and total number of receptors expressed on osteoblasts (). In the bone cell population models of and , , 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.
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 (). The work focused on clinically observed catabolic bone diseases related to perturbations of PTH glandular secretion. Following the approach proposed by , 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 (). 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.
2 Methods
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 , and further developed and refined by Pivonka and co-workers (; ; ; ). 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 ().
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 (), the model captures the smeared effect of many simultaneously active remodeling sites (). 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.
FIGURE 1
2.1 Mathematical model for bone cell dynamics
The bone cell population model by consists of a system of three ordinary differential equations. It describes the time-dependent behavior of osteoblastic precursor cells , active osteoblasts and active osteoclasts in molar concentrations asThe differentiation rates of uncommitted osteoblasts, osteoblast precursors and osteoclast precursors are denoted by and , respectively; and are the apoptosis rates of active osteoblasts and osteoclasts. In line with the source publication (), the model assumes that osteoblastic precursor cells do not undergo apoptosis but only differentiate into active osteoblasts, reflecting the biological understanding that once mesenchymal stem cells commit to the osteoblastic lineage, they progress through all differentiation stages ().
The effects induced by are assumed to depend on the concentration of active osteoclasts aswith proportionality constant and dissociation binding constant for and the respective receptor (). is stored in the bone matrix and released during resorption by active osteoclasts. It promotes -differentiation and -apoptosis, but acts as a repressor on -differentiation. The detailed derivation of (Equation 2) can be found in the original publication ().
The function describes the binding effect of the free ligand RANKL to the corresponding receptor activator nuclear factor B (RANK) aswhereas the molar concentration of RANK-RANKL complexes and RANK are denoted by and , respectively. The receptor-ligand binding triggers osteoclast precursor differentiation. Active osteoblasts secrete osteoprotegerin (OPG), which inhibits this process by binding to RANKL, thereby preventing RANKL-RANK interaction and subsequent osteoclast activation. Further details on (Equation 3) are given in the Supplementary Material of this article.
Parathyroid hormone (PTH) influences the RANK-RANKL-OPG pathway catabolically. It binds to its receptor on osteoblasts, increasing RANKL and decreasing OPG concentration. Consequently, more osteoclast precursors are differentiated into active osteoclasts, thus increasing bone resorption. Thus, the ratio of occupied RANK - and consequently – depends on the PTH effect, which is quantified by .
The Michaelis-Menten function describes the fraction of occupied PTH receptors and is derived from a one-state receptor model. The free ligand PTH binds to its receptor PTH1R forming a complex, whereas also the reverse reaction is possible, as shown in Figure 2. The kinetic parameters and quantify binding and unbinding tendencies.
FIGURE 2
The receptor-ligand binding reaches equilibrium long before the bone cells react, which leads to the assumption of a steady state. The function is defined asThe molar PTH concentration,depends on the basal synthesis rate , rate of external PTH injection and degradation rate . assumed that basal PTH concentration remains constant in both healthy and disease states. When PTH levels are elevated—whether due to disease or injected PTH—this is modeled as a sustained, constant increase from the basal level over a predefined time interval, rather than fluctuating daily. Bone diseases related to the parathyroid gland are characterized by alterations of the pulsatile characteristics, e.g., baseline secretion, pulse height, duration of each pulse and time between successive pulses. This work aims to develop the first bone cell population model that incorporates these characteristics.
2.2 Two-state receptor model for pulsatile PTH secretion
PTH1R exhibits multiple conformational states that can be generally classified into active (sensitized) and inactive (desensitized) states, each characterized by distinct signaling responses upon ligand binding (; ). The receptor can undergo conformational changes between these states independent of ligand binding, either through covalent modification or simple conformational rearrangement (; ). The response of osteoblasts to PTH-PTH1R binding depends on the respective conformation state of the receptor, with one state characterized by a shorter response and the other by a prolonged signaling response after receptor-ligand binding.
presented a model that accounts for this phenomenon of PTH1R, based on a general formulation of a two-state receptor model by , further detailed in . We describe the essential parts of the model below [a detailed description can be found in ].
Based on the ability of PTH1R to change conformation independent of ligand binding, the model assumes that both receptor conformational states can transform into each other regardless of their binding status (). The ligand PTH binds to both active (sensitized) and inactive (desensitized) receptors to form complexes in the corresponding states , with reversible unbinding reactions also possible. Both unbound receptors and ligand-receptor complexes can migrate between conformational states, allowing for dynamic transitions between all possible receptor states. All of the described reactions are governed by kinetic constants with (Equation 7), as shown in Figure 3.
FIGURE 3
The receptor and complex concentrations can be summarized in a concentration vector as . The respective time-dependent dynamics are given by a system of differential equations asThe coefficient matrix describes the kinetics of the receptor-ligand binding depending on the PTH ligand concentration asThe kinetic constants are given in Supplementary Table S1.
The ligand concentration L (Equation 8) depends on time , which makes it possible to include the pulsatile PTH pattern. The periodic ligand concentration is approximated as a piece-wise constant function, following the original approach () based on the two-state receptor model formulation (; ). This square-wave approximation has been validated against more realistic exponential decay profiles, demonstrating that square-wave stimulation provides a satisfactory approximation to periodic signals with exponential decay (). This piece-wise constant formulation is commonly employed in bone remodeling models for representing both endogenous hormone pulses () and administered drug injection patterns (; ). The periodic ligand concentration is defined asfor , where is the maximum number of periods during one simulation. The pulse shape is determined by various characteristics: tonic concentration ; pulse height/pulsatile concentration ; duration of on-phase ; and off-phase . It holds that .
The dimensionality of system Equation 6 can be reduced to three by normalizing the receptor and complex concentrations relative to total receptor concentration. When expressing the concentrations of active receptors , active complexes , inactive complexes , and inactive receptors as normalized values relative to the total receptor concentration, these quantities must satisfy the conservation relation . This conservation assumption is justified as receptor binding/unbinding dynamics reach steady-state much faster than cellular response timescales (; ), and the total receptor number is assumed to be much larger than fluctuations due to receptor production and degradation during the time periods relevant for cellular activity.
The model considers time-dependent ligand concentrations, which we define as pulsatile free PTH based on the glandular secretion pattern. The corresponding cellular activity is given by a linear combination of receptors and complexes in both states asThe activity constants with determine the weight of the respective concentrations for the cellular activity. The respective values were chosen in line with the original publication . The activity function (Equation 9) reflects the pulsatile ligand characteristics.
We derive two activity constants based on : the integrated activity and the cellular responsiveness . Both capture various characteristics of ligand input and activity output (see Figure 4). The integrated activity corresponds to the area under the curve of one activity pulse above baseline within one time period after cellular adaptation. The value of is computed by integrating over the interval for sufficiently large for the pulses to have reached steady-state. The composite trapezoidal rule is used for numerical integration.
FIGURE 4
The cellular responsiveness depends on but incorporates several characteristics of the ongoing stimuli: pulse shape; duration of pulse and off-phase; adaptation of cells to ongoing stimuli. It is defined asThe first factor of Equation 10 relates the integrated activity to the integrated activity of a step increase of the same magnitude as the original pulse . The second factor relates to the duration of one pulse and its off-phase (period ). Figure 4 represents an exemplary pulsatile PTH release pattern, the corresponding activity function , the derived integrated activity, , and remaining determinants of such as and .
Given that the receptor-ligand binding reaches equilibrium in a very short time (faster than cellular response in the bone cell population model), we propose the use of the activity constant and to quantify the effect of PTH on osteoblasts. The constants possess the following key features: reflection of the effect of receptor-ligand binding considering two conformation states; consideration of pulsatile, time-dependent behavior of ligand and cellular activity; constant quantification of cellular activity. Thus, we can replace in the bone cell population model with either or , to implicitly include pulsatile, time-dependent ligand concentration and the two conformation states of PTH1R.
2.3 Bone cell population model and two-state receptor: semi-coupling
We present a semi-coupled approach to integrate the two-state receptor model (Section 2.2) with the bone cell population model (Section 2.1). With this integration, the new bone cell population model includes information incorporated in the cellular activity: characteristics of pulsatile PTH; cellular adaptation to ongoing stimuli; effects of a two-state receptor. The approach involves replacing the Hill-type function with either the integrated activity or the cellular responsiveness . Receptor-ligand binding equilibrates considerably faster than bone cell dynamics, which justifies treating Hill-type functions as steady-state constants rather than dynamic functions, as explained in details in the source publication . This steady-state assumption was already implicit in the model; therefore, explicitly replacing these functions with constants is both mathematically and conceptually appropriate. Due to differences in magnitude between , , and , a scaling approach is necessary. We propose two distinct methods for this scaling.
The first method involves determining scaling parameters and based solely on healthy state values,withwhereas Equation 12 refers to the original Hill-type function for healthy state. The corresponding parameter values are given in Supplementary Table S2. The constants and are computed from Equation 10 and preceding formulation based on the healthy pulse characteristics (Table 1). These parameters are subsequently applied across all disease states. The kinetic parameter quantify binding and unbinding tendencies of the one-state receptor model (Figure 2).
TABLE 1
| State | [nM] | [nM] | [min] | [min] | [min] |
|---|---|---|---|---|---|
| Healthy | 6.4 | 4.2 | 10.6 | ||
| HPT | 7.6 | 3.5 | 11.1 | ||
| OP | 5.2 | 24.6 | 29.8 | ||
| PMO | 6.4 | 4.2 | 10.6 | ||
| HyperC | 9.41 | 6.18 | 15.59 | ||
| HypoC | 3.27 | 2.14 | 5.41 | ||
| GIO | 6.08 | 3.99 | 10.07 |
Parameters for healthy and disease states for PTH concentration.
The pulse shape is determined by the off-phase, , and the height of the pulse, . The sum of the duration of off-phase, , and on-phase, , is the period . The values are taken from , where disease states were computed using relative changes from healthy people based on experimental data. The considered disease states are Hyperparathyroidism (HPT), Osteoporosis (OP), Postmenopausal Osteoporosis (PMO), Hypercalcemia (HyperC), Hypocalcemia (HypoC), Glucocorticoid-induced Osteoporosis (GIO).
The second method is based on a comprehensive optimization approach across all states, determining parameters and through minimization,Here, index represents all seven physiological states detailed in Table 1. This optimization was implemented using Python’s SciPy library, using the default Broyden-Fletcher-Goldfarb-Shanno (BFGS) quasi-Newton method.
To establish comparable disease states between the pulsatile two-state receptor model (Section 2.2) and the constant-secretion reference model (Section 2.1), we introduce a scaling parameter for each disease given in Table 1,This parameter, illustrated in Figure 5, modifies the PTH concentration equation (Equation 5) for disease modeling in the reference model,whereas refers to the basal concentration of the original model formulation. Following these scaling approaches, (Equation 4) is replaced with either , , or to analyze the effect of the different constants and scaling approaches. The activity constants are computed independently (Section 2.2), reflecting the faster equilibration of PTH-PTH1R binding compared to bone cell response dynamics. Consequently, the RANK occupancy ratio in the bone cell population model becomes dependent on these scaled activities for both healthy and disease states.
FIGURE 5
To further validate and compare the model approaches, we analyze the temporal evolution of bone volume fraction, , describing the bone volume per total volume. Following , we define the change in bone volume fraction with time aswhere and represent the formation and resorption rates (bone volume formed and resorbed per unit cell concentration per unit time), respectively. This equation captures how bone volume fraction changes over time as a result of the competing processes of bone formation by active osteoblasts and bone resorption by active osteoclasts . Integration of this equation yields the temporal bone volume fraction profiles. In homeostasis, an equal fraction of bone volume is formed and resorbed. Thus, knowing the bone formation rate at a particular bone site (femoral neck, lumbar vertebra, radius) in homeostasis, the resorption rate can be computed asbased on the steady-state concentration of osteoblasts and osteoclasts .
2.4 Analysis of pulse characteristics
With the new model formulation, we can analyze the effect of different pulse characteristics on bone cell dynamics. The objective is to compare the magnitude of the activity constants of disease states to the healthy state with varying duration of pulse on- and off-phase, and , respectively. The period remains fixed to the physiological value according to Table 1 and it must hold that . To maintain physiological relevance, and must be fulfilled. For this analysis, the phases and may take every value possible that fulfills the above requirements. Thus, the on- and off-phases are determined according towhereas the same results can be obtained if and are switched in the above formulation as both cover the entire range of 0 to .
3 Results
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.
3.1 Calibration of Hill-type function and activity constants
We follow the approach given in Section 2.3 to include pulsatile PTH and PTH1R in two conformation states in the bone cell population model. We solve both Equation 11 and the minimization problem in Equation 13 to consider the different orders of magnitude of the Hill-type function , integrated activity and cellular responsiveness .
Both approaches to achieve comparable order of magnitude between , and result in similar values for the scaling parameters (Table 2). The aligned activity constants equal the formerly used function for the healthy state and show only a small deviation for postmenopausal osteoporosis and hypercalcemia (Figure 6). Hyperparathyroidism, hypocalcemia and glucocorticoid-induced osteoporosis show larger deviations.
TABLE 2
| Scaling parameter | Parameter value [-] | ASE [-] | RMSE [-] | Relative RMSE [-] |
|---|---|---|---|---|
| 2.07·10−2 | — | — | — | |
| 5.29·10−4 | — | — | — | |
| 3.03·10−2 | 1.15·10−3 | 1.28·10−2 | 9.22·10−2 | |
| 7.17·10−4 | 8.72·10−4 | 1.12·10−2 | 8.04·10−2 |
Results of scaling of and to achieve same order of magnitude as formerly used .
We identified the parameters for the healthy state and for all seven states (healthy and six disease states) using a minimization approach. The absolute squared error (ASE) refers to the result of the minimization problem (Equation 13). After identification of the scaling parameters using minimization, we calculate root mean square error (RMSE) and RMSE relative to the range of .
FIGURE 6
Despite relying on a single identified parameter across all states, both activity constants and show consistent trends with the Hill-type function (Figure 6). This consistency is also reflected in the small absolute error (ASE) for both minimization problems, which can be found in a range of to for both and (Table 2). The root mean square error (RMSE) relative to the range of Hill-type functions for the given states shows that the typical error is about 9% of the full span of for cellular responsiveness and 8% for integrated activity. The closest match is found for hypercalcemia, whereas and show the largest deviation from for glucocorticoid-induced osteoporosis.
3.2 Cellular activity and bone cell dynamics
The cellular activity function is computed for a PTH pattern in healthy and disease state according to Table 1. Figure 7 shows that for hyperparathyroidism (HPT), the tonic secretion increases more than double the maximum of healthy PTH concentration. Not only does HPT change the amount of PTH secretion, but also the duration of PTH secretion. Compared with a healthy state, the characteristic pulse lasts longer and the off-phase is shorter, ultimately resulting in a longer period (on- and off-phase). In glucocorticoid-induced osteoporosis (GIO), where PTH shows a lower basal secretion, its pulses reach almost the same maximum concentration as in healthy controls. This creates a larger relative amplitude of the pulses in GIO, despite similar absolute peak values.
FIGURE 7
The activity function reflects every change of the pulsatile PTH pattern. As shown in Figure 7, the cellular activity of HPT also has a longer period and on-phase of the pulse compared to the healthy state. The pulse height is elevated, whereas the activity pulse below as well as above baseline activity is longer than the healthy activity. For GIO, the activity pulse is longer compared to the healthy pattern, resulting in a higher fraction above baseline activity.
In both healthy and diseased conditions, cellular adaptation to the stimulus occurs during the first two periods, with all subsequent pulses looking identical. As the integrated activity and cellular responsiveness are computed based on the activity function , they also reflect any alterations of pulse characteristics. For example, both and increased noticeably for HPT, reflecting both prolonged period and activity pulse (see Table 3). One pulse remains longer above baseline and has a higher maximum value compared to the healthy state, leading to an increase in and . For GIO, both and are almost doubled, reflecting the higher maximum activity.
TABLE 3
| Activity constant | Healthy | HPT | OP | PMO | HyperC | HypoC | GIO |
|---|---|---|---|---|---|---|---|
| 36.6479 | 98.3844 | 35.1993 | 33.2851 | 6.7742 | 184.7505 | 62.0019 | |
| 0.9376 | 2.1014 | 0.4973 | 0.8550 | 0.1741 | 4.4375 | 1.6052 |
Values of integrated activity, , and cellular responsiveness, , for healthy state and all included disease states: hyperparathyroidism (HPT), osteoporosis (OP), postmenopausal osteoporosis (PMO), hypercalcemia (HyperC), hypocalcemia (HypoC), glucocorticoid-induced osteoporosis (GIO).
After scaling to achieve same order of magnitude, can be replaced with either , , or for healthy and disease state. The healthy state corresponds to a steady state of the model described by Equation 1. While transition from a healthy state to a pathological state in terms of PTH glandular secretion pattern might take some time (i.e., months to years), in numerical simulations the switch from healthy to disease is set instantaneously at a defined time point, as shown in Figure 8 for the case of HPT and in Figure 9 for GIO.
FIGURE 8
FIGURE 9
The results of the semi-coupled bone cell population model for HPT as disease state are shown in Figure 8. The cellular concentrations for , and have consistent dynamics across both the cellular responsiveness (, ) and integrated activity compared to the previously used Hill-type function . All five approaches demonstrate similar characteristic patterns, the main difference is the magnitude of the cellular concentrations in steady-state and disease case for calibration considering all states . The steady-state is equivalent for and in line with the identification of the scaling parameters based on only the healthy state.
Regarding the relative distance of to and , respectively, the concentration of is closer to during the steady-state (before and after disease state) using the original model formulation compared to the second calibration approach . The curve starts and levels off closer to concentration when using the activity constants. The homeostasis values of all three cell types are different for the model using either , or to quantify the PTH effect. After the immediate switch to a disease state, the increase in concentration varies across all cell types, models and scaling approaches. The original formulation shows the highest peak of all cell concentrations, whereas the cellular responsiveness, , results in the lowest increase.
The results of the semi-coupled bone cell population model for GIO as disease state are shown in Figure 9.
The bone cell concentrations show different dynamics for the original model formulation and the novel semi-coupled approach. The original model formulation results in cell concentrations in a range of 1⋅10−4 and reduced dynamics. The onset of the disease state is barely visible in the time interval between 20 and 80 days. In contrast to the formerly used , the change from homeostasis to GIO is clearly reproduced in the novel approach. The maximum concentration of and is found close to 2⋅10−3 when using either cellular responsiveness or integrated activity . The second scaling approach results in lower cell concentrations after onset of the disease state compared to the alternative approach.
This is reflected in the change of bone volume fraction during disease state (Figure 10). The initial bone volume fraction is chosen as 0.3 for trabecular bone (); the temporal evolution is described by Equations 16, 17. For the sake of conciseness and due to the comparable outcomes of both approaches, we present only the results of the calibration method accounting for all disease states and the formerly used as reference. For GIO, the novel model approach using both cellular responsiveness and integrated activity results in a loss of bone volume between 0.6% and 0.7% after 60 days of disease simulation. The loss of bone volume in the original model formulation is not significant for this disease state - in line with the reduced bone cell response (Figure 9). The highest catabolic effect during simulation of GIO is observed for .
FIGURE 10
The catabolic effect of hyperparathyroidism (HPT) and hypocalcemia (HypoC) are reflected in all model approaches and highest for the original formulation (Figure 10). Osteoporosis (OP), hypercalcemia (HyperC) and postmenopausal osteoporosis (PMO) show minor positive or no deviations from homeostasis across all approaches.
3.3 Pulse characteristics and homeostasis
Following the methodology presented in Section 2.4, we vary on- and off-phase of the pulse whereas all other characteristics remain fixed to physiologically healthy values. On- and off-phase are restricted to sum to the period (Equation 18).
Figure 11 shows the resulting values of integrated activity and cellular responsiveness w.r.t. the ratio of on- and off-phase. The curve is symmetrical around the point where equals , represented by the logarithmic ratio . The maximum of both and at this point corresponds to the highest catabolic response with healthy tonic and pulsatile PTH. Both activity constants are close to zero when is much larger than or the other way around.
FIGURE 11
For both and , the healthy value (cyan circle) appears near the maximum of both curves, occurring at a slightly positive ratio (see Figure 11). Hypocalcemia (purple star) exhibits the highest response in both metrics, which stand distinctly above other pathological conditions. Most disease states cluster around similar ratios, indicating a common temporal pattern in PTH and receptor dynamics across different pathological conditions. However, osteoporosis (orange triangle) presents a notable exception to this trend. Not only does it occur at a different ratio compared to other conditions, but it also shows the largest difference when comparing and responses.
Figure 12 shows the results of the bone cell population model for healthy state and HPT when replacing not only with the calibrated cellular responsiveness, but with the maximum cellular responsiveness resulting from the pulse characteristics and for healthy state (see Figure 11). This results in elevated baseline concentrations across all cell populations (, , and ), while showing a proportionally lower catabolic jump in hyperparathyroidism. We selected the cellular responsiveness instead of because its response curve is steeper around the maximum. This provides better discrimination of the effect of different pulse characteristics, as evidenced by the distinct positioning of, for example, osteoporosis (see Figure 11). The temporal dynamics remain consistent between both configurations.
FIGURE 12
4 Discussion
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 (). 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 incorporating the bone volume fraction balance equation (). 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 (; ; ; ) 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 , can now be achieved by replacing with the calibrated cellular responsiveness or integrated activity from the two-state PTH-PTH1R model. While the latter quantities were previously given theoretical bone-related interpretations in the two-state receptor model (), 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 and cellular responsiveness – 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 (). 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 (; ; ).
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 (; ; ). 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 (). 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 (), 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 (), 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 (; ).
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, 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, 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 (). 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. 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 (; ; ), experimental studies have demonstrated variability in individual osteoclast resorptive activity (), and mathematical modeling of injury repair has provided in vivo evidence for time-variable cellular activity rates (). 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 (). 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 and subsequently improved by 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 (; ).
Third, our square-wave pulses for the PTH secretion pattern based on the original formulations (; ) 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 (; ). 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 (). The activity function 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 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.
5 Nomenclature
5.1 Resource identification initiative
All simulations were performed using Python Programming Language (RRID:SCR_008394).
Statements
Data availability statement
The datasets presented in this study can be found in online repositories (). The names of the repository/repositories and accession number(s) can be found below: https://github.com/cmodiz/Bone-Models.git.
Author contributions
CM: Software, Conceptualization, Writing – review and editing, Writing – original draft, Formal Analysis, Visualization, Methodology. NMC: Conceptualization, Project administration, Supervision, Methodology, Writing – review and editing. SS: Conceptualization, Supervision, Writing – review and editing. JM-R: Supervision, Conceptualization, Writing – review and editing. JLC-G: Conceptualization, Supervision, Writing – review and editing. VS: Supervision, Project administration, Conceptualization, Writing – review and editing. SM: Writing – review and editing, Conceptualization, Supervision. PP: Conceptualization, Project administration, Methodology, Writing – review and editing, Writing – original draft, Supervision.
Funding
The author(s) declare that financial support was received for the research and/or publication of this article. We gratefully acknowledge funding received through the Australian Research Council (IC190100020; DP230101404; ARC-FT180100338).
Acknowledgments
We thank the main author of the source publication, Denisa Martonova, for sharing the original resources.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
The author(s) declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.
Generative AI statement
The author(s) declare that no Generative AI was used in the creation of this manuscript.
Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbioe.2025.1619276/full#supplementary-material
References
1
Aibar-AlmazánA.Voltes-MartínezA.Castellote-CaballeroY.Afanador-RestrepoD. F.Carcelén-FraileM. d. C.López-RuizE. (2022). Current status of the diagnosis and management of osteoporosis. Int. J. Mol. Sci.23, 9465. 10.3390/ijms23169465
2
AraujoA.CookL. M.LynchC. C.BasantaD. (2014). An integrated computational model of the bone microenvironment in bone-metastatic prostate cancer. Cancer Res.74, 2391–2401. 10.1158/0008-5472.can-13-2652
3
BaratchartE.LoC. H.LynchC. C.BasantaD. (2022). Integrated computational and in vivo models reveal key insights into macrophage behavior during bone healing. PLOS Comput. Biol.18, e1009839. 10.1371/journal.pcbi.1009839
4
Ben-awadhA. N.Delgado-CalleJ.TuX.KuhlenschmidtK.AllenM. R.PlotkinL. I.et al (2014). Parathyroid hormone receptor signaling induces bone resorption in the adult skeleton by directly regulating the rankl gene in osteocytes. Endocrinology155, 2797–2809. 10.1210/en.2014-1046
5
BilezikianJ. P.BandeiraL.KhanA.CusanoN. E. (2018). Hyperparathyroidism. Lancet391, 168–178. 10.1016/S0140-6736(17)31430-7
6
BiselloA.ChorevM.RosenblattM.MonticelliL.MierkeD. F.FerrariS. L. (2002). Selective ligand-induced stabilization of active and desensitized parathyroid hormone type 1 receptor conformations. J. Biol. Chem.277, 38524–38530. 10.1074/jbc.M202544200
7
BolingE. P. (2004). Secondary osteoporosis: underlying disease and the risk for glucocorticoid-induced osteoporosis. Clin. Ther.26, 1–14. 10.1016/s0149-2918(04)90001-x
8
BonadonnaS.BurattinA.NuzzoM.BugariG.RoseiE. A.ValleD.et al (2005). Chronic glucocorticoid treatment alters spontaneous pulsatile parathyroid hormone secretory dynamics in human subjects. Eur. J. Endocrinol.152, 199–205. 10.1530/eje.1.01841
9
Brucker-DavisF.ThayerK.ColbornT. (2001). Significant effects of mild endogenous hormonal changes in humans: considerations for low-dose testing. Environ. Health Perspect.109, 21–26. 10.1289/ehp.01109s121
10
CaprianiC.IraniD.BilezikianJ. P. (2012). Safety of osteoanabolic therapy: a decade of experience. J. Bone Min. Res.27, 2419–2428. 10.1002/jbmr.1800
11
ChelohaR. W.GellmanS. H.VilardagaJ.-P.GardellaT. J. (2015). PTH receptor-1 signalling—mechanistic insights and therapeutic prospects. Nat. Rev. Endocrinol.11, 712–724. 10.1038/nrendo.2015.139
12
ChiavistelliS.GiustinaA.MazziottiG. (2015). Parathyroid hormone pulsatility: physiological and clinical aspects. Bone Res.3, 14049. 10.1038/boneres.2014.49
13
ChristiansenP.SteinicheT.VesterbyA.MosekildeL.HessovI.MelsenF. (1992). Primary hyperparathyroidism: iliac crest trabecular bone volume, structure, remodeling, and balance evaluated by histomorphometric methods. Bone13, 41–49. 10.1016/8756-3282(92)90360-9
14
CookC. V.LightyA. M.SmithB. J.Ford VersyptA. N. (2024). A review of mathematical modeling of bone remodeling from a systems biology perspective. Front. Syst. Biol.4, 1368555. 10.3389/fsysb.2024.1368555
15
DattaN. S.Abou-SamraA. B. (2009). Pth and pthrp signaling in osteoblasts. Cell. Signal.21, 1245–1254. 10.1016/j.cellsig.2009.02.012
16
DingM.OdgaardA.HvidI. (1999). Accuracy of cancellous bone volume fraction measured by micro-ct scanning. J. Biomech.32, 323–326. 10.1016/S0021-9290(98)00176-6
17
GardellaT. J. (2020). The parathyroid hormone receptor type 1 (Humana, Cham), chap. 16. Contemporary endocrinology. 3rd edn.Cham, Switzerland: Springer Nature Switzerland AG, 323–347. 10.1007/978-3-319-69287-6_16
18
HarmsH. M.NeubauerO.KayserC.WüstermannP. R.HornR.BrosaU.et al (1994a). Pulse amplitude and frequency modulation of parathyroid hormone in early postmenopausal women before and on hormone replacement therapy. J. Clin. Endocrinol. Metab.78, 48–52. 10.1210/jcem.78.1.8288712
19
HarmsH. M.SchlinkeE.NeubauerO.KayserC.WüstermannP. R.HornR.et al (1994b). Pulse amplitude and frequency modulation of parathyroid hormone in primary hyperparathyroidism. J. Clin. Endocrinol. Metab.78, 53–57. 10.1210/jcem.78.1.8288713
20
HawkinsF.GarlaV.AlloG.MalesD.MolaL.CorpasE. (2021). “Chapter 5 - senile and postmenopausal osteoporosis: pathophysiology, diagnosis, and treatment,” in Endocrinology of aging. Editor CorpasE. (Amsterdam, Netherlands: Elsevier), 131–169. 10.1016/B978-0-12-819667-0.00005-6
21
Hernández-CastellanoL.HernandezL.BruckmaierR. (2020). Review: endocrine pathways to regulate calcium homeostasis around parturition and the prevention of hypocalcemia in periparturient dairy cows. Anim14, 330–338. 10.1017/s1751731119001605
22
HofbauerL. C.ZeitzU.SchoppetM.SkalickyM.SchülerC.StolinaM.et al (2009). Prevention of glucocorticoid-induced bone loss in mice by inhibition of rankl. Arthritis Rheum.60, 1427–1437. 10.1002/art.24445
23
Jakubas-PrzewłockaJ.PrzewłockiP. (2005). Assessment of changes due to the long-term effect of estrogen and calcium deficiency in the trabecular bone structure in rats. Clin. Exp.l Rheumatol.23, 385–388.
24
KanehisaJ.HeerscheJ. (1988). Osteoclastic bone resorption: in vitro analysis of the rate of resorption and migration of individual osteoclasts. Bone9, 73–79. 10.1016/8756-3282(88)90106-8
25
KenkreJ.BassettJ. (2018). The bone remodelling cycle. Ann. Clin. Biochem.55, 308–327. 10.1177/0004563218759371
26
KomarovaS. V.SmithR. J.DixonS.SimsS. M.WahlL. M. (2003). Mathematical model predicts a critical role for osteoclast autocrine regulation in the control of bone remodeling. Bone33, 206–215. 10.1016/S8756-3282(03)00157-1
27
LauffenburgerD. A.LindermanJ. J. (1993). Receptors: models for binding, trafficking, and signaling. New York, NY: Oxford University PressNew. 10.1093/oso/9780195064667.001.0001
28
LavaillM.TrichiloS.ScheinerS.ForwoodM. R.CooperD. M. L.PivonkaP. (2020). Study of the combined effects of pth treatment and mechanical loading in postmenopausal osteoporosis using a new mechanistic pk-pd model. Biomech. Model. Mechan.19, 1765–1780. 10.1007/s10237-020-01307-6
29
LemaireV.TobinF. L.GrellerL. D.ChoC. R.SuvaL. J. (2004). Modeling the interactions between osteoblast and osteoclast activities in bone remodeling. J. Theor. Biol.229, 293–309. 10.1016/j.jtbi.2004.03.023
30
LiY.GoldbeterA. (1989). Frequency specificity in intercellular communication. influence of patterns of periodic signaling on target cell responsiveness. Biophys. J.55, 125–145. 10.1016/S0006-3495(89)82785-7
31
LiM.ZhangX.LiS.GuoJ. (2024). Unraveling the interplay of extracellular domain conformational changes and parathyroid hormone type 1 receptor activation in class b1 g protein-coupled receptors: integrating enhanced sampling molecular dynamics simulations and markov state models. ACS Chem. Neurosci.15, 844–853. 10.1021/acschemneuro.3c00747
32
LoC. H.BaratchartE.BasantaD.LynchC. C. (2021). Computational modeling reveals a key role for polarized myeloid cells in controlling osteoclast activity during bone injury repair. Sci. Rep.11, 6055. 10.1038/s41598-021-84888-1
33
LocklinR. M.KhoslaS.TurnerR. T.RiggsB. L. (2003). Mediators of the biphasic responses of bone to intermittent and continuously administered parathyroid hormone. J. Cell. Biochem.89, 180–190. 10.1002/jcb.10490
34
LuL.TianL. (2023). Postmenopausal osteoporosis coexisting with sarcopenia: the role and mechanisms of estrogen. J. Endocrinol.259, e230116. 10.1530/JOE-23-0116
35
MarinoS.BellidoT. (2024). Pth receptor signalling, osteocytes and bone disease induced by diabetes mellitus. Nat. Rev. Endocrinol.20, 661–672. 10.1038/s41574-024-01014-7
36
Martínez-ReinaJ.PivonkaP. (2019). Effects of long-term treatment of denosumab on bone mineral density: insights from an in-silico model of bone mineralization. Bone125, 87–95. 10.1016/j.bone.2019.04.022
37
MartonováD.LavaillM.ForwoodM. R.RoblingA.CooperD. M. L.LeyendeckerS.et al (2023). Effects of pth glandular and external dosing patterns on bone cell activity using a two-state receptor model—implications for bone disease progression and treatment. PLOS One18, e0283544. 10.1371/journal.pone.0283544
38
[Dataset]ModizC. (2025). Bone models. Brisbane, Australia: Queensland University of Technology. 10.25912/RDF_1742263438470
39
OkazakiM.FerrandonS.VilardagaJ.-P.BouxseinM. L.PottsJ. T.GardellaT. J. (2008). Prolonged signaling at the parathyroid hormone receptor by peptide ligands targeted to a specific receptor conformation. Proc. Natl. Acad. Sci. U. S. A.105, 16525–16530. 10.1073/pnas.0808750105
40
PetersonM. C.RiggsM. M. (2010). A physiologically based mathematical model of integrated calcium homeostasis and bone remodeling. Bone46, 49–63. 10.1016/j.bone.2009.08.053
41
PivonkaP.KomarovaS. V. (2010). Mathematical modeling in bone biology: from intracellular signaling to tissue mechanics. Bone47, 181–189. 10.1016/j.bone.2010.04.601
42
PivonkaP.ZimakJ.SmithD. W.GardinerB. S.DunstanC. R.SimsN. A.et al (2008). Model structure and control of bone remodeling: a theoretical study. Bone43, 249–263. 10.1016/j.bone.2008.03.025
43
PivonkaP.Calvo-GallegoJ. L.SchmidtS.Martínez-ReinaJ. (2024). Advances in mechanobiological pharmacokinetic-pharmacodynamic models of osteoporosis treatment – pathways to optimise and exploit existing therapies. Bone186, 117140. 10.1016/j.bone.2024.117140
44
PotterL. K.GrellerL. D.ChoC. R.NuttallM. E.StroupG. B.SuvaL. J.et al (2005). Response to continuous and pulsatile pth dosing: a mathematical model for parathyroid hormone receptor kinetics. Bone37, 159–169. 10.1016/j.bone.2005.04.011
45
RobertsW. E.MozsaryP. G.KlinglerE. (1982). Nuclear size as a cell-kinetic marker for osteoblast differentiation. Am. J. Anat.165, 373–384. 10.1002/aja.1001650403
46
Ruiz-LozanoR.Calvo-GallegoJ. L.PivonkaP.McDonaldM. M.Martínez-ReinaJ. (2024). An in silico approach to elucidate the pathways leading to primary osteoporosis: age-related vs. postmenopausal. Biomech. Model. Mechan.23, 1393–1409. 10.1007/s10237-024-01846-2
47
RyserM. D.NigamN.KomarovaS. V. (2009). Mathematical modeling of spatio-temporal dynamics of a single bone multicellular unit. J. Bone Mineral Res.24, 860–870. 10.1359/jbmr.081229
48
RyserM. D.KomarovaS. V.NigamN. (2010). The cellular dynamics of bone remodeling: a mathematical model. SIAM J. Appl. Math.70, 1899–1921. 10.1137/090746094
49
SalariN.GhasemiH.MohammadiL.BehzadiM. h.RabieeniaE.ShohaimiS.et al (2021). The global prevalence of osteoporosis in the world: a comprehensive systematic review and meta-analysis. J. Orthop. Surg. Res.16, 609. 10.1186/s13018-021-02772-0
50
SchaeferF. (2000). Pulsatile parathyroid hormone secretion in health and disease. Chichester, England: John Wiley and Sons, Ltd., 225–243. 10.1002/0470870796.ch13
51
Schappacher-TilpG.CherifA.FuertingerD. H.BushinskyD.KotankoP. (2019). A mathematical model of parathyroid gland biology. Physiol. Rep.7, e14045. 10.14814/phy2.14045
52
ScheinerS.PivonkaP.HellmichC. (2013). Coupling systems biology with multiscale mechanics, for computer simulations of bone remodeling. Comput. Methods Appl. Mech. Eng.254, 181–196. 10.1016/j.cma.2012.10.015
53
SchmittC. P.HömmeM.SchaeferF. (2005). Structural organization and biological relevance of oscillatory parathyroid hormone secretion. Pediatr. Nephrol.20, 346–351. 10.1007/s00467-004-1767-7
54
SegelL. A.GoldbeterA.DevreotesP. N.KnoxB. E. (1986). A mechanism for exact sensory adaptation based on receptor modification. J. Theor. Biol.120, 151–179. 10.1016/S0022-5193(86)80171-0
55
SiddiquiJ. A.JohnsonJ.Le HenaffC.BitelC. L.TamasiJ. A.PartridgeN. C. (2017). Catabolic effects of human pth (1–34) on bone: requirement of monocyte chemoattractant protein-1 in murine model of hyperparathyroidism. Sci. Rep.7, 15300. 10.1038/s41598-017-15563-7
56
SunM.WuX.YuY.WangL.XieD.ZhangZ.et al (2020). Disorders of calcium and phosphorus metabolism and the proteomics/metabolomics-based research. Front. Cell Dev. Biol.8, 576110. 10.3389/fcell.2020.576110
57
TrichiloS.ScheinerS.ForwoodM.CooperD. M.PivonkaP. (2019). Computational model of the dual action of pth — application to a rat model of osteoporosis. J. Theor. Biol.473, 67–79. 10.1016/j.jtbi.2019.04.020
58
VeldhuisJ. D. (2008). Pulsatile hormone secretion: mechanisms, significance and evaluation. Dordrecht: Springer Netherlands, 229–248. 10.1007/978-1-4020-8352-5_10
59
WangB.YangY.Abou-SamraA. B.FriedmanP. A. (2009). Nherf1 regulates parathyroid hormone receptor desensitization: interference with β-arrestin binding. Mol. Pharmacol.75, 1189–1197. 10.1124/mol.108.054486
60
ZavalaE. (2022). Misaligned hormonal rhythmicity: mechanisms of origin and their clinical significance. J. Neuroendocrinol.34, e13144. 10.1111/jne.13144
61
ZavalaE.WedgwoodK. C.VoliotisM.TabakJ.SpigaF.LightmanS. L.et al (2019). Mathematical modelling of endocrine systems. Trends Endocrinol. Metab.30, 244–257. 10.1016/j.tem.2019.01.008
62
ZhuS.ChenW.MassonA.LiY.-P. (2024). Cell signaling and transcriptional regulation of osteoblast lineage commitment, differentiation, bone formation, and homeostasis. Cell Discov.10, 71. 10.1038/s41421-024-00689-6
Summary
Keywords
parathyroid hormone, parathyroid hormone/parathyroid hormone-related protein receptor, bone cell dynamics, disease modeling, pulsatile signal characteristics
Citation
Modiz C, Castoldi NM, Scheiner S, Martínez-Reina J, Calvo-Gallego JL, Sansalone V, Martelli S and Pivonka P (2025) Computational simulations of endocrine bone diseases related to pathological glandular PTH secretion using a multi-scale bone cell population model. Front. Bioeng. Biotechnol. 13:1619276. doi: 10.3389/fbioe.2025.1619276
Received
28 April 2025
Accepted
01 September 2025
Published
01 October 2025
Volume
13 - 2025
Edited by
Eiji Tanaka, Tokushima University, Japan
Reviewed by
Laurent Pujo-Menjouet, Université Claude Bernard Lyon 1, France
Etienne Baratchart, Lund University, Sweden
Updates
Copyright
© 2025 Modiz, Castoldi, Scheiner, Martínez-Reina, Calvo-Gallego, Sansalone, Martelli and Pivonka.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Corinna Modiz, corinna.modiz@hdr.qut.edu.au
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.