In addition to Wnt and RANK-RANKL-OPG signaling, two cytokines play essential signaling roles in bone remodeling: macrophage colony-stimulating factor (MCSF) and TGF-β. MCSF stimulates osteoclastogenesis by binding to monocyte progenitor cells and preosteoclasts. This triggers intracellular cascades that induce NFATc1, the main transcription factor for osteoclastogenesis (Guder et al., 2020). The importance of MCSF stems from its role in stimulating the first stage of osteoclastogenesis (Figure 2), which RANKL does not stimulate, and the proliferation of osteoclast precursor cells. The role of TGF-β is less straightforward. Its regulatory effects are concentration dependent (Janssens et al., 2005; Wu et al., 2016). Moreover, it regulates both bone formation and bone resorption. Low concentrations of TGF-β stimulate osteoclast production while promoting preosteoblast migration and proliferation. High concentrations inhibit osteoclastogenesis. Also, high concentrations inhibit preosteoblast migration and late-stage osteoblast differentiation. Although these contradictory findings still puzzle researchers, the mechanism of changes in TGF-β concentration during remodeling is well understood. Inactive TGF-β is stored in the extracellular matrix of bone (Epsley et al., 2021). As osteoclasts remove bone, TGF-β is released and activated, increasing the concentration of TGF-β (Janssens et al., 2005; Matsumoto and Abe, 2011).

Numerous other cytokines also regulate bone remodeling. TGF-β is not the only cytokine released from the bone matrix during resorption; others include IGF-1, IGF-2, and BMP2 (Sims and Martin, 2020). Typically, cytokines are classified as osteoclastogenic or osteoblastogenic, though their roles may be concentration-dependent as described with TGF-β. Bone resorption is inhibited by anti-inflammatory cytokines such as interleukin (IL)-4, IL-10, IL-13, IL-18, and interferon (IFN)-γ (Walsh et al., 2006). Conversely, it is stimulated by pro-inflammatory cytokines such as IL-1, IL-1β, TNF-α, IL-6, IL-11, IL-15, and IL-17. These cytokines modulate RANKL and Wnt signaling to increase osteoclast activity (Walsh et al., 2006; Tilg et al., 2008). For example, TNF-α upregulates the expression of RANKL, DKK1, and sclerostin in osteocytes (Kitaura et al., 2020; Epsley et al., 2021). The influx of RANKL promotes osteoclastogenesis, while the influx of DKK1 and sclerostin inhibits osteoblastogenesis. Cytokines interact in a complex network with RANK-RANKL-OPG, Wnt, MCSF, and TGF-β to regulate bone remodeling.

For further detailed coverage of cytokines and growth factors that locally regulate bone cell function and thorough diagrams of signaling ligands and their associated receptors, intracellular kinases and transcription factors, and biological outcomes, readers are referred to Plotkin and Bivi (2014).

Immune cells contribute to bone homeostasis through cytokine expression and direct immune cell activity. Osteoclasts are derived from innate immune cells called monocytes (Saxena et al., 2021). Monocytes are more commonly differentiated into macrophages and dendritic cells. Studies have shown that these cells can transdifferentiate into preosteoclasts (Bonomo et al., 2016; Saxena et al., 2021; Srivastava and Sapra, 2022). Macrophages further modulate bone remodeling through the expression of inflammatory cytokines IL-1, IL-6, and TNF-α or bone formation factors IL-10, BMP-2, and TGF-β (Fischer and Haffner-Luntzer, 2022). In contrast, dendritic cells primarily stimulate osteoclastogenesis through RANK-RANKL activation of T cells, which upregulates T cell production of RANKL, IL-1, IL-6, IL-17, and TNF-α (Bonomo et al., 2016). However, not all T cells are osteoclastogenic.

Different populations of T cells affect osteoclasts and osteoblasts in different ways. Naïve CD4⁺ T cells can differentiate into osteoclastogenic subtypes, e.g., T helper (Th)17 and Th9 cells, or anti-osteoclastic subtypes, e.g., Th1, Th2, and T regulatory (Treg) cells, characterized by their cytokine expression profiles (Guder et al., 2020). For example, Th17 cells express high levels of IL-17, which upregulates RANK in preosteoclasts and RANKL in osteoblasts, increasing bone resorption (Fischer and Haffner-Luntzer, 2022; Srivastava and Sapra, 2022). Th17 cells also secrete IL-6, RANKL, and TNF-α to promote osteoclastogenesis and suppress osteoblast activity (Srivastava et al., 2018; Epsley et al., 2021). Cytokine profiles of Th1, Th2, and Treg cells contrast the profiles from Th17 cells. These cells secrete anti-osteoclastic cytokines such as IFN-γ, IL-4, TGF-β1, and IL-10 (Okamoto et al., 2017; Srivastava et al., 2018; Guder et al., 2020). However, following the pattern of Wnt, TGF-β, and other cytokines, T cell roles are not always clear. Activated Tregs secrete DKK1, which inhibits Wnt-mediated bone formation (Lehmann et al., 2021). This inhibitory effect contrasts studies showing Tregs increase Wnt-10b production by CD8⁺ T cells (Tyagi et al., 2018). Despite this, studies indicate a balance between Th17 and Treg cells is important for healthy bone remodeling, such that higher Th17 to Treg ratios contribute to rheumatoid arthritis and osteoporosis (Okamoto et al., 2017; Srivastava et al., 2018). Declining bone health is associated with many classic inflammatory diseases, such as periodontitis, rheumatoid arthritis, and aseptic prosthesis (Epsley et al., 2021). To obtain a more complete picture of bone remodeling, it is vital to consider these complex bone-immune interactions.

PTH is a systemic hormone that regulates calcium levels in the blood in part by triggering calcium release from the bone. Chief cells within the parathyroid gland produce PTH when serum calcium levels are low (Chaiya and Rattanakul, 2017). The increase in circulating PTH triggers bone remodeling, and subsequently, osteoclasts release calcium from the bone into the blood to maintain homeostasis (Peterson and Riggs, 2010; Coelho et al., 2016). This ability to stimulate remodeling has led to the development of synthetic PTH for osteoporosis treatment. However, PTH is another signaling factor with a dual role in osteoclastogenesis. Circulating PTH stimulates osteoclast activity by increasing the RANKL to OPG ratio but inhibits osteoclast formation by decreasing sclerostin and DKK1 (Silva and Bilezikian, 2015; Kenkre and Bassett, 2018; Wein and Kronenberg, 2018; Guder et al., 2020). Exogenous PTH alters bone remodeling differently depending on the administration schedule. Continuous administration decreases overall levels of bone density, whereas intermittent administration increases bone density levels (Coelho et al., 2016; Lemaire and Cox, 2019). This has led many researchers to develop mathematical models to understand the mechanisms of PTH regulation.

The dual role of PTH is currently understood to involve both stimulation of preosteoblast production and inhibition of preosteoblast differentiation (Chaiya and Rattanakul, 2017). As a result, preosteoblast cell populations increase while osteoblast populations remain unchanged. This causes an increase in the RANKL to OPG ratio since osteoblastic cells produce more OPG and less RANKL (Gori et al., 2000) as they mature. A short burst of PTH stimulates remodeling by increasing RANKL and suppressing OPG. High concentrations of PTH over a long period, as in the case of hyperparathyroidism, dysregulate bone remodeling due to the overproduction of osteoclasts. This leads to a larger resorption cavity that the limited number of osteoblasts cannot fill. This interrelationship between bone cells and PTH exemplifies the complexity of the bone remodeling process.

Estrogen and bone health have been closely linked for decades due to the correlation between postmenopausal estrogen decline and bone loss (Khosla et al., 2012). Although early research on the mechanism of estrogen regulation of bone remodeling was unclear, recent studies in osteoimmunology have improved our understanding. Estrogen deficiency increases bone turnover and unbalanced remodeling (Khosla et al., 2012). This occurs through estrogen-mediated inhibition of RANKL production and stimulation of OPG expression, which limits osteoclastogenesis (Eriksen, 2010; Florencio-Silva et al., 2015; Noirrit-Esclassan et al., 2021). Estrogen has also been shown to prevent apoptosis of osteoblasts and osteocytes (Florencio-Silva et al., 2015). This is consistent with studies showing that estrogen deficiency induces osteocyte apoptosis (Tomkinson et al., 1997; Tomkinson et al., 1998; Emerton et al., 2010; Khosla et al., 2012; Delgado-Calle and Bellido, 2022). Ovariectomized (OVX) murine experiments demonstrate that estrogen directly supports bone formation by upregulating Wnt-10b in bone marrow stromal cells (Perkins et al., 2023).

Further evidence of the estrogen-bone link is based on the presence of estrogen receptors (ER) on bone cells and targeted deletion studies. ERα is found on osteoclastic and osteoblastic cells, while ERβ is found on osteoblasts (Sharma et al., 2020). Targeted deletion of osteoblastic ERα in murine models led to low bone mass in both males and females (Almeida et al., 2017; Gao et al., 2021). The targeted deletion of ERα in osteoclasts and osteoclast progenitor cells increased osteoclast numbers in females but not in males (Almeida et al., 2017). Another study of ERβ deletion in mesenchymal stem cells found that bone mass increased only in female rodents (Almeida et al., 2017). These knockout studies indicate that estrogen signaling is vital to bone homeostasis in males and females, with sex-based differences in these signaling mechanisms.

Estrogen regulates bone remodeling through direct and immune-mediated mechanisms (Khosla et al., 2012). For instance, estrogen protects against T-cell-mediated bone loss by upregulating Wnt signaling. While mice with DKK1-expressing T cells experienced OVX-induced bone loss, knockout mice without DKK1-expressing T cells did not, and prior to OVX, these mice exhibited higher bone mass (Lehmann et al., 2021). The loss of bone in response to estrogen deficiency is recognized as a cytokine-driven process involving T cell populations such as Tregs and Th17 cells that results in the bone resorption activity of osteoclasts exceeding that of bone-forming osteoblasts (Pacifici, 2012). The anti-inflammatory effect of estrogen extends to macrophages and dendritic cells. Estrogen deficiency has been shown to induce the transdifferentiation of pro-inflammatory M1 macrophages into osteoclasts and increase the ratio of M1 to anti-inflammatory M2 macrophages (Saxena et al., 2021). Without estrogen, more dendritic cells were shown to express IL-7 and IL-15, which upregulates IL-17 and TNF-α production by T cells (Saxena et al., 2021). Furthermore, a cross-sectional clinical study of postmenopausal women showed that elevated inflammatory markers such as IL-6, IL-β, and TNF-α were negatively correlated with bone mass (Damani et al., 2023). These findings indicate estrogen regulates bone remodeling through immune-mediated effects and direct signaling within the bone microenvironment.

Pharmaceuticals can indirectly regulate remodeling while treating various diseases, or they can be designed to target mechanisms in the bone remodeling cycle intentionally. Glucocorticoids are anti-inflammatory agents that are used broadly but have negative effects on bone health by decreasing osteoblast and osteocyte populations and increasing osteoclast survival (Hardy et al., 2018). Despite this, they are often used to reduce chronic inflammation in rheumatoid arthritis that can otherwise lead to osteoporosis. Teriparatide, a PTH analog that has anabolic effects on the bone, has several proposed mechanisms for its action (Wein and Kronenberg, 2018). Since teriparatide requires expensive daily injections, it is used mainly for severe osteoporosis or those who need to use glucocorticoids long-term for other conditions (Hodsman et al., 2005).

Estrogen replacement therapy (ERT) is a pharmaceutical intervention directly designed to impact bone remodeling based on the association between estrogen decline and osteoporosis after menopause. Hormone replacement therapy (HRT) augments estrogen with progestogens. Selective estrogen receptor modulators (SERMS) act as estrogen receptor agonists in some tissues like bone and antagonists in other tissues, sometimes detrimentally (Ellis et al., 2015). These treatments have different and controversial risks associated with breast cancer, coronary heart disease, and stroke that impact their adoption based on individualized management of benefits and risks (Manson et al., 2013; Ellis et al., 2015; Hodis and Sarrel, 2018; Faubion et al., 2022; Onwude, 2022; Nudy et al., 2023).

Another group of medications to reduce bone loss are antiresorptives, which target signaling mechanisms of bone remodeling that contribute to osteoclast activity. Bisphosphonates are antiresorptives that are currently the most common treatments for bone loss. These drugs inhibit bone resorption by inducing osteoclast apoptosis and reducing osteoclast activity (Berkhout et al., 2015; Coelho et al., 2016; Aibar-Almazán et al., 2022). Bisphosphonates even alter bone remodeling after treatment is terminated because they bind to hydroxyapatite crystals on the surface of the bone matrix (Drake et al., 2008; Aibar-Almazán et al., 2022). They can be released from the surface in subsequent remodeling cycles (Coelho et al., 2016). Bisphosphonates are generally well tolerated but are most often discontinued due to gastrointestinal distress or concerns about side effects such as osteonecrosis of the jaw or spiral fractures of the femur midshaft (Aibar-Almazán et al., 2022). The monoclonal antibody denosumab is a newer antiresorptive. Denosumab inhibits resorption by blocking RANK-RANKL binding. It acts as an OPG mimic, binding to RANKL to prevent osteoclast activation. Although denosumab is more effective at preventing bone loss than bisphosphonates and used for metastatic cancers that target bone, there is a higher risk of osteonecrosis of the jaw (Aibar-Almazán et al., 2022; Fu et al., 2023) and higher frequency of second tumors in cancer patients on denosumab (Stopeck et al., 2010; Fizazi et al., 2011; Henry et al., 2011; Raje et al., 2018; Tovazzi et al., 2019).

Romosozumab is the newest pharmaceutical intervention for the bone remodeling cycle. Romosozumab is a monoclonal antibody that binds to sclerostin, allowing Wnt ligands to activate the canonical pathway, stimulate bone formation, and inhibit bone resorption (Figure 2). The disadvantage is that romosozumab is associated with more undesirable side effects than bisphosphonates including increased risk of adverse cardiovascular events (Asadipooya and Weinstock, 2019; Aibar-Almazán et al., 2022). Overall, targeted treatments of the bone remodeling cycle have poor compliance and high discontinuation rates due to a combination of high costs, unwanted side effects, and psychological factors (Aibar-Almazán et al., 2022). Viable new treatments need to eliminate or reduce these concerns.

Regardless of the approach, models of bone remodeling generally include the dynamics of osteoblast and osteoclast cells. While cell populations’ evolution and signaling interactions vary between models, osteoblasts always form bone, and osteoclasts always break down bone. The net effect of bone regulation by these cells can generally be represented in ODE form by (Frost, 1964b) d z d t = R f − R r = k f O B L − k r O C L (6) where z is bone volume fraction, often corresponding to the ratio of bone volume to tissue volume (BV/TV) bone histomorphometric quantity (Bouxsein et al., 2010; Dempster et al., 2013), TV is frequently the trabecular volume (Burr and Akkus, 2014), R _f is the formation rate, R _r is the resorption rate, k _f is the formation rate constant for the action of the osteoblasts, and k _r is the resorption rate constant for the osteoclast activity. The variables OBL and OCL in Eq. 6 usually represent changes from the steady state population, sometimes called active cell populations. Additionally, bone volume may be replaced with bone mass or other relevant bone properties. Bone volume, total osteoblast population, and total osteoclast population cannot have negative values.

A thorough understanding of osteoclast bone resorptive activity and population dynamics is crucial to predicting how much bone is resorbed during a remodeling cycle. The difference between the net formation and degradation terms determines the osteoclast population dynamics. The power law approach is used to mathematically represent these dynamics following Eq. 1 as (Komarova et al., 2003) d O C L d t = α OCL O C L g 11 O B L g 21 − β OCL O C L (7) where OCL represents the number of osteoclasts, OBL represents the number of osteoblasts, g ₁₁ represents autocrine (osteoclast to osteoclast) signaling action, g ₂₁ represents paracrine (osteoblast to osteoclast) signaling action, and α _OCL and β _OCL represent proliferation and degradation rate constants, respectively. The proliferation of osteoclasts (j = 1) depends on autocrine from osteoclasts (i = 1) and paracrine from osteoblasts (i = 2) signaling effects. The degradation rate of osteoclasts is assumed to be proportional to the current population.

Some power law models modify the signaling dynamics to account for a specific molecular factor by reformulating the population dynamics and recalculating general signaling exponents (Komarova, 2005; Graham et al., 2013). For instance, if a signaling factor F _OCL alters osteoclast proliferation, the osteoclast equation is modified to become d O C L d t = F OCL α OCL O C L g 11 O B L g 21 − β OCL O C L (8) with new values for g ₁₁ and g ₂₁ as compared to Eq. 7.

By the mass action kinetics approach following Eq. 5, the population of osteoclasts is given by (Pivonka et al., 2008) d O C L d t = α pOCL p O C L π act,pOCL X 1 − β OCL O C L π act,OCL X 2 (9) where OCL is the osteoclast population, pOCL is the preosteoclast population, π act,pOCL X 1 is the activation from signaling factor X ₁, π act,OCL X 2 is the activation from signaling factor X ₂, and α _pOCL and β _OCL represent differentiation and apoptosis rate constants, respectively. In Pivonka et al. (2008), differentiation and apoptosis terms are both activated and not inhibited. X ₁ is RANKL, and X ₂ is TGF-β. Different combinations of activating and repressing π terms are proposed in models from other publications. When following the mass action kinetics approach, researchers typically examine individual signaling factors during formulation. As a result, the overall structure of the osteoclast equation rarely undergoes drastic changes. Instead, new signaling factors are simply added through more π terms in the mass actions kinetics approach.

Uncommitted monocytes and preosteoclasts are rarely modeled as dynamic populations (thus, pOCL is a constant in Eq. 9). Osteoclasts are assumed to differentiate from a large pool of hematopoietic stem cells, so the uncommitted population is usually modeled as a fixed quantity. Although this assumption is reasonable for healthy bone remodeling, it loses validity when studying diseases where hematopoietic stem cell numbers are reduced (Weilbaecher, 2000; Matatall et al., 2016). Preosteoclasts are usually omitted for simplification under the assumption that remodeling is already occurring, i.e., the activation stage is assumed to occur instantaneously (Figure 1). However, this neglects the time needed to initiate this remodeling stage.

Mathematical models must include osteoblast cell dynamics to understand changes in bone formation rates. Osteoblast population balances are similar to those of osteoclasts given in Eqs 7, 9, and these balances are modeled by the power law approach following Eq. 1 as (Komarova et al., 2003) d O B L d t = α OBL O C L g 12 O B L g 22 − β OBL O B L (10)

and by the mass action kinetics approach following Eq. 5 as (Pivonka et al., 2008) d O B L d t = α pOBL p O B L π rep,pOBL Y 1 − β OBL O B L (11) where the parameters here correspond to osteoblast (cell type j = 2) dynamics and pOBL is the preosteoblast population. In Pivonka et al. (2008), differentiation is inhibited, and apoptosis is neither activated nor inhibited. Y ₁ is TGF-β. Unlike osteoclasts, osteoblasts are modeled with one or two consumption terms. The use of one consumption term encapsulates osteoblast apoptosis and its conversion to other bone cells, such as osteocytes and bone lining cells. When models include two consumption terms, one tracks osteoblast conversion to osteocytes. The other consumption term tracks osteoblast apoptosis and other osteoblast cell losses.

Another difference between osteoblast and osteoclast population balances is that preosteoblast dynamics are commonly modeled. A study analyzing a generalized model of bone remodeling highlights the importance of preosteoblast populations (Zumsande et al., 2011). For instance, preosteoblasts release key signaling molecules that initiate the resorption phase of bone remodeling. Additionally, preosteoblast cell dynamics must be modeled because the number of osteoblasts is dictated by the preosteoblast population after proliferation (Buenzli et al., 2012b). Since osteoblasts do not proliferate, a decrease in bone formation may result from fewer preosteoblasts. Preosteoblast population balances follow the same form as those of osteoblasts and osteoclasts. The formation term represents differentiation from uncommitted osteoblasts, while the consumption term represents conversion to osteoblasts.

As research continues to indicate that osteocytes are essential coordinators of bone remodeling, it is vital to include their dynamics and populations in mathematical models. Osteocyte population balances are less commonly found in mathematical models than osteoclast and osteoblast balances but generally follow similar principles. In power law models, osteocytes are modeled following Eqs 8, 10 as (Graham et al., 2013; Cook et al., 2022) d O C Y d t = F OCY α OCY O B L g 23 (12) where OCY is the osteocyte population, α _OCY is the rate of conversion from osteoblasts, and g ₂₃ is osteoblast signals that influence the production of osteocytes (cell type j = 3) via osteoblast embedding. The factor F _OCY represents osteocyte signaling that activates and terminates the bone remodeling cycle.

Differently, models following the mass action kinetics approach base their osteocyte population on the change in bone volume as (Martin et al., 2019; Martin et al., 2020; Calvo-Gallego et al., 2023) d O C Y d t = η d z d t (13) where η is the average concentration of osteocytes embedded in the bone matrix and z is the bone volume fraction. Note the lack of degradation terms for long-lived osteocytes in Eqs 12, 13. Some disease or injury conditions may explicitly induce loss of osteocytes, which can be incorporated by including a loss term or by reducing the osteocyte initial condition (Graham et al., 2013; Cook et al., 2022).

Most mathematical models of bone remodeling have overlooked reversal cells and bone lining cells. Their absence is reasonable, given that these cells’ importance and mechanistic behavior are not well understood. Moreover, the modeling work of Arias et al. (2018) suggests that reversal and bone lining cells are not required to capture the dynamics of bone remodeling.

These cells are only included in one model each. Trichilo et al. (2019) models a constant population of bone lining cells. They also include a dynamic parameter that describes bone lining cell differentiation into osteoblasts. This parameter varies with PTH dosage, introducing another mechanistic avenue for PTH to regulate osteoblastogenesis. The mononuclear cells modeled in Martin and Buckland-Wright (2004) remove collagen fibrils from the bone surface during resorption. The behavior of these mononuclear cells is in line with current understanding of the role that reversal cells play in bone remodeling, indicating that the so-called mononuclear cells in Martin and Buckland-Wright (2004) are reversal cells. While reversal and bone lining cells have historically been excluded, future mathematical models paired with experimental work could help provide mechanistic insights into their functions.

As mentioned earlier, Wnt plays an important role in modulating bone health. However, few models consider the details of this interaction. When mathematical models consider the Wnt/sclerostin interaction, Wnt levels are often excluded or assumed constant (Supplementary Tables S1–S3). Instead, these models focus on sclerostin levels (Figure 6) (Graham et al., 2013; Eudy et al., 2015; Lemaire and Cox, 2019; Martin et al., 2019). Two models with dynamic Wnt concentrations only allow Wnt to be altered through the presence of a tumor (Buenzli et al., 2012b; Farhat et al., 2017). Cook et al. (2022) includes a generalized dynamic concentration of Wnt-10b, where the amount of exogenous Wnt-10b (from dietary sources) influences BMU cell populations and bone volume.

While the RANK-RANKL-OPG and Wnt pathways are key regulators of the bone remodeling cycle, other cytokines modulate these signals and bone cell activity. Despite this, there is a distinct lack of variety in the cytokines considered in mathematical formulations of bone remodeling. Table 5 shows that only three cytokines aside from RANK-RANKL-OPG and TGF-β are explicted modeled: IL-6, MCSF, and IGF. One benefit of including other cytokines in bone models is the potential to explore their importance under various remodeling conditions, yet many models with cytokines lack this analysis. For example, IGF is included in Garzón-Alvarado (2012) to simulate osteosclerosis because tumors are known to increase IGF levels and thus increase osteoblast activity. However, this work does not directly analyze the effect of IGF on bone cell dynamics. Although the modeling and analysis of IL-6 and MCSF is limited, some studies analyze their role in remodeling using perturbation or sensitivity analysis.

Only two of the five models that include dynamic IL-6 levels analyze its effect on bone cell dynamics during the bone remodeling cycle (Table 5). IL-6 is included in Kroll (2000) in their study of the effects of PTH on bone dynamics, albeit in a simplified manner. Following Kroll (2000), Idrees et al. (2019) adapts Komarova et al. (2003) to include IL-6 in the simulation of intermittent versus continuous PTH treatment. Whereas Kroll (2000) scales theoretical parameter values so that osteoblast counts remain below 1000 cells, Idrees et al. (2019) performs a meta-analysis of various experiments to estimate parameter values statistically. This is an improvement over many other bone models that extract experimental values from multiple studies or rely on sparse and disparate clinical data sets.

The dynamics of IL-6 in PTH treatment models are simplistic compared to multiple myeloma models. Wang et al. (2011) accounts for IL-6 stimulation of multiple myeloma cells and IL-6 activation of RANKL. Wang et al. (2011) also performs a perturbation analysis to investigate the relative degree of RANKL activation by PTH versus IL-6 in homeostatic and diseased remodeling states and finds that PTH dominates over IL-6 under healthy bone remodeling conditions. However, IL-6 activates RANKL more than PTH under elevated IL-6 conditions, representing multiple myeloma disease. A limitation here is that the model lacks other cytokines and mechanisms that can alter RANKL activation, which may dominate over PTH and IL-6.

Wang et al. (2011) is extended by Ji et al. (2014) without any change in the representation of IL-6. However, the models differ in their parameter estimation and sensitivity analysis approaches. Wang et al. (2011) fits the model to achieve adult bone and cancer cell densities corresponding to literature values, while Ji et al. (2014) uses genetic algorithms. Although both studies include sensitivity analysis, Ji et al. (2014) performs sensitivity analysis at a fixed time point rather than over time as in Wang et al. (2011). Nonetheless, the results in Ji et al. (2014) support those in Wang et al. (2011). Upon varying 11 parameter values from 50% to 150%, Ji et al. (2014) shows that bone volume is most sensitive to TGF-β and the progression of multiple myeloma disease is most sensitive to IL-6. Given these results and the biological evidence that IL-6 contributes to bone pathophysiology, exploring its inclusion in future mathematical models of bone remodeling is pertinent.

The cytokine MCSF is a key activator of osteoclastogenesis, yet analysis of its effects via mathematical models of bone remodeling is more sparse than IL-6. Martin and Buckland-Wright (2004) is the first model that includes the role of MCSF. This model investigates the depth and duration of osteoclast erosion during resorption and models MCSF as a scalar variable but assumes it is always present. A sensitivity analysis of resorption depth indicates that changes in MCSF levels are equivalent to changes in maximum osteoclast activity, and both effects are minor compared to TGF-β. Lerebours et al. (2016) and Pivonka et al. (2013) are biomechanical models that include MCSF activation of uncommitted osteoclasts. However, they assume the MCSF concentration is constant, resulting in a continuous activator function term.

Proctor and Gartland (2016) investigates the kinetics of MCSF via network-based ODEs. As with Lerebours et al. (2016) and Pivonka et al. (2013), this work investigates the effect of mechanical loading and the effects of PTH and circadian rhythm. In contrast, the network analysis includes multiple parameters that capture the role of MCSF in remodeling (outlined in Table 5). Sensitivity analysis of the model shows that the secretion rate of MCSF by osteoblasts and preosteoclasts and the MCSF degradation rate result in a change of more than 5% in bone mass. When the rate of degradation doubles, bone mass increases by more than 60%. One limitation of this study is the assumption that MCSF secretion rates are considered equal across cell types. Still, these results warrant further mathematical investigation of MCSF in homeostatic and pathological bone remodeling. Ultimately, future mathematical models of bone biology should explore the complex and coordinated role of cytokines, growth factors, and hormones in bone remodeling.

Despite the established interest in osteoimmunology, few bone models include immune cells (Table 6). Most bone models that investigate the role of immune cells do so in the context of bone injury and repair (Trejo et al., 2019; Tourolle et al., 2021; Baratchart et al., 2022). Of particular interest is the model by Baratchart et al. (2022). Hypothesis testing of candidate models determines the interaction of monocytes, macrophages, injury factors, and inflammatory factors in bone cell dynamics. The model is supported by biological data, with the parameters calibrated with one set of experimental data and validated with another. These methods show how researchers can elucidate the mechanisms of complex bone-immune dynamics using mathematical models. However, bone injury and repair models describe the acute healing process of fractures, which has different signaling pathways than the continuous bone renewal or remodeling process for homeostasis and skeletal integrity over a lifetime. Therefore, bone healing models cannot be directly applied to the bone remodeling processes.

The few mathematical models of bone remodeling that incorporate immune cells are outlined in Table 6. Although the models by Akchurin et al. (2008) and Proctor and Gartland (2016) seemingly include immune cells (mononuclear cells and hematopoietic stem cells), these are simply different osteoclast progenitors. As mentioned earlier, these cells are often lumped into a general class of uncommitted osteoclasts or preosteoclasts.

Of the models listed in Table 6, only Islam et al. (2021) investigates the dynamic effect of multiple immune cells in bone remodeling. The work includes a three-compartment physiologically based PK model for differentiating naïve CD4⁺ T cells into Treg cells in the gut, blood, and bone. These Tregs then influence TGF-β production in the bone and induce Wnt-10b production. The physiologically based PK model is then linked to a bone remodeling model that includes the local effects of systemic changes in Wnt-10b (Cook et al., 2022). Since this is the only mathematical description of nonlocal immune effects on bone dynamics, significant opportunities remain for future research to explore multi-organ systemic interactions between the skeletal and immune systems.

Despite its documented importance in bone remodeling and estrogen-deficient osteoporosis, the incorporation of dynamic estrogen levels in mathematical models of bone remodeling is underwhelming (Table 7). Analysis of the 88 cell population-based bone modeling publications in Tables 1–4 reveals that roughly a third (30 of 88) mention or model estrogen. Of these 30 models, half (15) capture estrogen effects in their mathematical model, while the other half only mention estrogen briefly. Most models that mention estrogen cite evidence that estrogen deficiency is involved in osteoporosis or modulates bone remodeling (Kroll, 2000; Marathe et al., 2008; Pivonka et al., 2008; Marathe et al., 2011; Garzón-Alvarado, 2012; Ross et al., 2012; Buenzli, 2015; Chen-Charpentier and Diakite, 2016; Lee and Okos, 2016; Jerez et al., 2018; Bahia et al., 2020; Javed et al., 2020). The remaining articles that mention estrogen acknowledge that estrogen is not incorporated into their model or that integration of estrogen is an opportunity for future models (Peterson and Riggs, 2010; Buenzli et al., 2012b; Coelho et al., 2016; Martin et al., 2019). Altogether, this indicates that only about 15% of all bone remodeling models mathematically account for estrogen-induced biochemical changes in bone cell dynamics.

All models with estrogen effects are osteoporosis-specific models, and their respective mathematical representations are outlined in Table 7. Several models do not consider dynamic estrogen levels (Lemaire et al., 2004; Scheiner et al., 2013; Lemaire and Cox, 2019; Lemaire and Cox, 2019; Martin et al., 2019; Trichilo et al., 2019; Larcher and Scheiner, 2021). Instead, they model the effects of estrogen deficiency by altering RANKL, OPG, PTH, and TGF-β. Lemaire et al. (2004) and Lemaire and Cox (2019) manually lower the OPG and TGF-β production parameters for osteoporotic scenarios. Trichilo et al. (2019) and Martin et al. (2019) increase RANKL levels with a RANKL dosage term fitted to OVX rat data and clinical postmenopause data. Scheiner et al. (2014) accounts for more estrogen-deficiency effect by using π terms to capture disease-related increases in RANKL, decreases in mechanical loading sensitivity, and denosumab competition with RANK and OPG.

Two models with dynamic estrogen levels aim to determine the most effective therapeutic dose to prevent bone loss. In Rattanakul et al. (2003), periodic estrogen treatment is modeled with a linear increase in osteoclast removal. Chaiya and Rattanakul (2017) reformulates the model to be explicitly piecewise. Using the power law approach to illustrate, the osteoclast rate equation becomes d O C L d t = α O C L g 11 O B L g 21 − β O C L if  t ≠ n T α O C L g 11 O B L g 21 − β O C L − ρ O C L if  t = n T (14) where ρ is the parameter related to estrogen treatment, T is the prescribed dose time, and n is the treatment number. In addition to Eq. 14, Chaiya and Rattanakul (2017) adds a constant term to capture the osteoblast-stimulating effects of estrogen treatment on the osteoblast rate equation. The motivation of this work is to understand alternative treatment regimes for ERT because long-term continuous treatment, while effective in increasing bone volume, has been shown in some studies to increase the risk for breast cancer and heart disease (Levin et al., 2018).

Models of dynamic estrogen loss often use an exponential decay equation for estrogen concentration or an estrogen-dependent dynamic parameter (Schmidt et al., 2011; Post et al., 2013; Berkhout et al., 2015; Berkhout et al., 2016). Although Berkhout et al. (2015) notes that estrogen concentration could better explain disease dynamics, they opt for the decay equation instead due to the high uncertainty in their model. In addition to estrogen decay, Schmidt et al. (2011) and Post et al. (2013) modify the OPG production parameter following Lemaire et al. (2004). Javed et al. (2018) derives an alternate formulation for remodeling altogether to simplify the mass action kinetics models proposed by Lemaire et al. (2004) and Pivonka et al. (2008). Estrogen changes are represented by a hyperbolically scaled estrogen term that modulates the RANKL state variable. In contrast, Jorg et al. (2022) models estrogen as an age-dependent concentration with a characteristic time scale of menopause onset that is fit to clinical data. This model also considers how estrogen alters sclerostin levels.

Most models of estrogen dynamics only consider bone cell or RANK-RANKL-OPG interactions, despite accumulating evidence that estrogen modulates other mechanisms of bone remodeling such as Wnt signaling and immune-bone interactions (see Section 3.2.2). Future models need to incorporate these complex dynamics in the ongoing effort to improve mechanistic understanding of estrogen-deficient osteoporosis.

Many bone remodeling models explore the effects of one or more drugs on bone health (Riggs and Cremers, 2019; Ait Oumghar et al., 2020) (Supplementary Tables S1–S3). Researchers typically start by modeling a healthy or diseased remodeling cycle (or leveraging existing models) and then extend the process to include drug effects. For example, glucocorticoid therapies and their interactions with the bone remodeling cycle are modeled in Lemaire et al. (2004), Schmidt et al. (2011), and Lemaire and Cox (2019). These models alter one parameter related to a symptom of glucocorticoid treatment, specifically reduced osteoblast populations. This essentially involves reducing α _OBL in Eq. 10 or Eq. 11. However, reducing one parameter corresponding to an effect observed with glucocorticoid treatment is a simplistic approach that may miss important mechanistic impacts on the bone remodeling cycle.

Other models (e.g., Jorg et al. (2022)) study bisphosphonates, denosumab, or romosozumab. The antiresorptive drugs are modeled by combining the PKPD information of the drug of interest and an already-established mathematical model of the BMU. PK information consists of factors that explain how the drug disperses in the body. These are usually differential equations that track the amount of a drug in a target area. PD information describes how the drug interacts with the body. The effects can be shown directly through new parameters in the model or implicitly applied by changing an existing parameter.

Whereas the immune-bone connection gained traction in the 2000s, the link between gut and bone metabolism is more recent. So it follows that fewer mathematical models of bone remodeling consider gut-mediated impacts on bone health. Only one mathematical model of bone remodeling incorporates gut and immune cells (Islam et al., 2021). This model explores butyrate treatment of bone through T-cell-mediated changes in Wnt-10b. Although much is still unknown about the gut-bone connection, the Islam et al. (2021) model is initialized with data from mouse experiments that complement the mathematical model. Sensitivity analysis and in silico hypothesis generation link the calculated parameters to experimental conditions that can be modified to explore new treatments. This highlights the benefit of experimentally supported mathematical models of bone remodeling. The multi-compartment modeling approach of Islam et al. (2021) and Peterson and Riggs (2010) provide examples of how mathematical models of bone remodeling may explore relationships of systemic multi-organ effects.

Similar to pharmaceutical modeling, most cancer models start by modeling normal bone homeostasis and supplement it with an equation for tumor dynamics (Marathe et al., 2008; Ayati et al., 2010; Araujo et al., 2014; Ji et al., 2014; Coelho et al., 2016; Farhat et al., 2017). Many cancer models also add or adjust parameters such as RANKL, TGF-β, and PTH, which are known to be modified by tumors. The populations of these cancer tumor cells (T) are usually modeled in one of two ways. The first modeling method is based on growth curves, as in d T d t = γ T d e n s i t y λ − η T (15) where γ and η are growth and decay parameters, Tdensity is a relationship between the current and maximum cancer cell population, and λ is an additional relationship term to capture the effects of other cancer interactions considered (Ayati et al., 2010; Buenzli et al., 2012b; Coelho et al., 2016; Zhang and Mager, 2019; Miranda et al., 2020). For example, in the Coelho et al. (2016) model, λ in Eq. 15 corresponds to the concentration of osteoclasts.

The second common way to model cancer populations follows the mass action kinetics approach. Here, the populations are controlled by different signaling factors represented by π terms (Wang et al., 2011; Ji et al., 2014). In a publication that uses both modeling methods, the growth curve is better for early cancer, and the mass action kinetics structure is better for established cancer (Farhat et al., 2017).

The mechanisms of tumor growth and metastasis have important spatial considerations. This is the primary motivation behind existing bone remodeling PDEs and ABMs of cancer (Ayati et al., 2010; Ryser et al., 2012; Araujo et al., 2014), which track the movement of cancer cells in space.