Mathematical modeling of bone remodeling after surgical menopause

TL;DR

This paper introduces a mathematical model to understand how bone loss changes after surgical menopause, focusing on the role of osteocytes in bone remodeling.

Contribution

The novelty is a new model capturing increased bone loss and long-term rebound after surgical menopause, incorporating osteocyte dynamics.

Findings

The model shows increased bone loss in the first 15 years after surgical menopause.

A rebound in bone mineral density is predicted in the long-term post-surgery.

Osteocyte regulation is identified as a key factor in bone remodeling after estrogen deficiency.

Abstract

Osteoporosis is a skeletal pathology characterized by decreased bone mass and structural deterioration resulting from an imbalance in bone metabolic processes. Estrogen deficiency in postmenopausal women leads to an increased risk of osteoporosis, while women who have undergone complete oophorectomies display an even higher risk due to the sudden decrease in estrogen. Some evidence indicates that bone loss slows in the period beyond 15 years after surgery; however, there is substantial uncertainty in clinical data. To explore the effects of surgically induced menopausal transition, here we propose a mathematical model for the bone cell dynamical responses to sudden estrogen deficiency, which extends an existing model for osteoporosis due to aging and natural menopause. Using data on key effects observed in female mice and humans after bilateral oophorectomy, this new model considers the role of osteocytes embedded within the mineralized bone matrix in regulating bone remodeling, which results in net bone loss after surgical menopause. The model parameter values in natural and surgical menopause were estimated from aggregated human clinical data from existing longitudinal studies. The new model effectively captures the previously unmodeled increase in bone loss during the first 15 years post-surgical menopause and the rebound in bone mineral density in the long-term. With this model, effects of treatments on targeting osteocyte dynamics could be explored in the future.

Figures

(a,b) Bone mineral density (BMD) measured in the lumbar spine of women (except Looker et al. (1998) are from the hip). Data are normalized to premenopausal levels from the time of menopause onset, and error bars represent the standard deviation. Normalization process for each dataset is detailed in Section 2.1. The data are from Pansini et al. (1995); Ohta et al. (2002); Hadjidakis et al. (2003); Looker et al. (1998); Hibler et al. (2016); Chittacharoen et al. (1999); Yasui et al. (2007). Dashed lines show linear fits to the data with slopes listed in the legends. (a) Natural menopause and (b) surgical menopause. (c) Comparison of estrogen decline in natural and surgical menopause, shown over 30 years after menopause onset. The inset shows the decline over the 3 days following menopause onset.

Schematic of mathematical model species in the bone remodeling process. The transformations of cells through differentiation, production, and degradation are illustrated by black dashed arrows. Signaling-related inhibition interactions are shown by red flat-head arrows, and activation interactions are shown as green solid arrows. Parameters that govern interactions are shown near the corresponding arrows. The schematic shows estrogen as an inhibitor of osteoclast differentiation and sclerostin production. Surgical menopause-induced changes are depicted by red dashed arrows with scissors icons. Created with BioRender.com.

(a) Reparameterization of the Jörg et al. (2022) model using more sources of natural menopause (NM) data. The dashed red curve shows the original model, the solid red curve shows the new parameterization, and the black solid curve shows the model with sudden estrogen loss alone, without any new effects. (b) Extension of the model to surgical menopause. New effects are parameterized using long-term (blue curve) and short-term (green curve) data. The shaded regions highlight model sensitivity to parameters. (c) Cell populations from the surgical menopause model using the three parameter sets shown in panel (b). NM: Natural menopause. SM: surgical menopause.

Steady BMD change as a function of sclerostin production levels at steady state (relative to the premenopause level). Steady-state cell concentrations are calculated assuming that no new effects are included in the model. Sclerostin production is then adjusted by a percentage of its premenopause level to calculate the resulting steady-state change in BMD.

Tables

TABLE 1. Summary of bone mineral density (BMD) data curated for parameterization of model. Data are taken from cross-sectional studies that use dual-energy X-ray absorptiometry to measure BMD.

Study source	BMD source	Number of women (type of menopause)
Pansini et al. (1995)	L2-L4 vertebrae	160 (natural), 67 (surgical)
Ohta et al. (2002)	L2-L4 vertebrae	20 (natural), 20 (surgical)
Hibler et al. (2016)	L1-L4 vertebrae	53 (surgical)
Hadjidakis et al. (2003)	L2-L4 vertebrae	177 (natural), 210 (surgical)
Chittacharoen et al. (1999)	L1-L5 vertebrae	309 (natural), 102 (surgical)
Chittacharoen et al. (1999)	L2-L4 vertebrae	141 (surgical)

TABLE 2. Initial conditions taken at steady state at 30 years before menopause onset. All variables are in dimensionless form.

Variable	Meaning	Initial value
PB	Preosteoblasts	2.16×102
PC	Preosteoclasts	4.10×103
C	Osteoclasts	4.20×101
B	Osteoblasts	1.07×102
S	Osteocytes	6.12×102
Sc	Sclerostin	1.12×104
Bd	Bone density	1

TABLE 3. Natural menopause (NM) model parameters with parameters taken from Jörg et al. (2022) or estimated using the procedures outlined in Section 2.3. ND: no dimensions.

Parameter	Meaning	Value	Dimension	Source
ωPC	Differentiation rate of PC→C	0.93	day−1	Jörg et al. (2022)
ωPB	Differentiation rate of PB→B	0.32	day−1	Jörg et al. (2022)
ωB	Differentiation rate of B→S	6.4×10−4	day−1	Jörg et al. (2022)
ηB	Apoptosis rate of B	8.678×10−3	day−1	Jörg et al. (2022)
ηC	Apoptosis rate of C	1.096×10−4	day−1	Jörg et al. (2022)
ηS	Apoptosis rate of S	1.1×10−4	day−1	Jörg et al. (2022)
λB	First-order bone formation rate	1.29×10−6	day−1	Jörg et al. (2022)
λC	First-order bone resorption rate	3.82×10−6	day−1	Jörg et al. (2022)
κSc	Degradation rate of Sc	0.05	day−1	Jörg et al. (2022)
τE	Characteristic time of estrogen decline	2.6	year	Jörg et al. (2022)
ePC	Threshold for estrogen inhibition of PC→C	0.2556	ND	Estimated, NM
eSc	Threshold for estrogen inhibition of Sc production	10.59	ND	Estimated, NM
scPC	Threshold for sclerostin activation of PC→C	8.6×106	ND	Jörg et al. (2022)
scPB	Threshold for sclerostin inhibition of PB→B	1.63×102	ND	Jörg et al. (2022)
scΩ	Threshold for sclerostin inhibition of Bd	3.04×103	ND	Jörg et al. (2022)
rΩ	Threshold for resorption activation of Bd	1.02×103	ND	Jörg et al. (2022)
BMC0	Steady-state BMC	0.8	ND	Jörg et al. (2022)
νΩ	Maximum relative effect of resorption	1.08×103	ND	Jörg et al. (2022)

TABLE 4. Surgical menopause (SM) model parameters from sources, calculated as described in Section 2.2, and estimated using the procedures outlined in Section 2.3 for short-term (15 years or less post-surgery) and long-term (up to 30 years post-surgery) data. ND: no dimensions.

Parameter	Meaning	Value	Dimension	Source
κE	Estrogen degradation rate post-surgery	6.1996	day−1	Ginsburg et al. (1998)
ksyn	Estrogen synthesis rate post-surgery	0.6	day−1	Calculated
Esurg(tm+30)	Estrogen concentration 30 days post-surgery	0.096	ND	Bellanti et al. (2013)
ηsurg	Apoptosis of S post-surgery (short-term)	5	ND	Estimated, SM short-term
τ	Timescale of post-surgery dynamics (short-term)	9.7×10−3	day-1	Estimated, SM short-term
ωsurg	Enhanced C differentiation post-surgery (short-term)	1.86	ND	Estimated, SM short-term
ηsurg	Apoptosis of S post-surgery (long-term)	0.4174	ND	Estimated, SM long-term
τ	Timescale of post-surgery dynamics (long-term)	0	day−1	Estimated, SM long-term
ωsurg	Enhanced C differentiation post-surgery (long-term)	0.2155	ND	Estimated, SM long-term

Full text

Introduction

1

Bone tissue is continuously resorbed and formed through the bone remodeling process, which maintains healthy tissue and repairs micro-fractures in the skeleton. At homeostasis, the biomechanical, biochemical, and cellular mechanisms involved in remodeling of the adult skeleton maintain bone mass. However, any alteration in this complex bone turnover cycle can result in changes in bone (i.e., pathologic or anabolic) (Allen and Burr, 2014). The cellular mechanisms that determine bone health occur in functional locations called basic multicellular units (BMUs) within which cellular interactions contribute to bone tissue remodeling through a continuous cycle of activation, resorption, and formation (Robling et al., 2006). The three main types of cells that contribute to this cycle are osteoclasts, osteoblasts, and osteocytes. Osteocytes, widely recognized for their strain sensitivity, are the most abundant of these cells and signal to recruit other cells to the BMU to initiate bone resorption and formation (Creecy et al., 2021). One such signal is sclerostin (Suen and Qin, 2016; Plotkin and Bellido, 2016; Delgado-Calle et al., 2017), which negatively regulates bone formation by inhibiting Wnt (Li et al., 2005; Krause et al., 2010; Atkins et al., 2011; Bellido, 2014; Kim et al., 2020), thus reducing osteoblastogenesis. Sclerostin has also been shown to contribute to bone resorption by upregulating secretion of receptor activator of nuclear factor $[eqn]$ B (RANK) ligand (RANKL) from osteocytes (Wijenayaka et al., 2011; Suen and Qin, 2016), which stimulates osteoclastogenesis (Nakashima et al., 2011; Tomkinson et al., 1998; Fujiwara et al., 2016; Plotkin and Bellido, 2016; Karlafti et al., 2019). Osteoclasts degrade the bone protein matrix and solubilize the mineral hydroxyapatite, and osteoblasts initiate bone matrix formation by forming osteoid tissue, which is later mineralized into bone (Kenkre and Bassett, 2018). Up to 20% of the osteoblasts within osteoid tissue differentiate further into osteocytes (Parfitt, 1976; Martin et al., 1998; Plotkin and Bellido, 2016; Delgado-Calle and Bellido, 2022). After the resorbed bone tissue is replaced through bone formation, osteocytes signal to BMU cells to slow bone formation (Kenkre and Bassett, 2018).

Perturbations of the bone remodeling process can cause an imbalance between catabolic and anabolic activity, leading to bone pathologies characterized by substantial bone loss. In particular, osteoporosis is a low-bone-density disease caused by such an imbalance and leads to increased fracture risk, reducing quality of life for those affected and imposing a significant financial burden on the global economy (Harvey et al., 2010). Osteoporosis is prevalent in postmenopausal women, with estimates suggesting that between 30% and 40% of women over age 50 are affected by low bone mass or osteoporosis (Rachner et al., 2011; Sarafrazi et al., 2021). Estrogen-deficient bone loss has been linked to several metabolic processes in bone remodeling. The presence of estrogen has been shown in human studies to prevent apoptosis of osteocytes and osteoblasts, and estrogen reduces levels of sclerostin (Mödder et al., 2010). Estrogen also reduces the impact of osteoclasts by preventing osteoclastogenesis and increasing osteoclast apoptosis (Florencio-Silva et al., 2015). When estrogen levels temporarily drop during perimenopause or are permanently low after menopause, osteoclast differentiation increases and causes more bone removal (Hsu et al., 2024); meanwhile, despite an initial increase in osteoblasts with the decline in estrogen, the bone-forming osteoblast activity is unable to match the pace of the increase in resorption, resulting in less bone mass (Jilka et al., 1998; Almeida et al., 2007; Khosla, 2013; Seeman, 2013; Karlamangla et al., 2021).

For women with intact ovaries, highly variable and declining estrogen production throughout perimenopause and menopause leads to increased bone loss. During late perimenopause and the early postmenopausal period, human cohort studies have shown that bone is lost at a rate of 2%–2.4% per year in the spine and 1.2%–1.7% per year in the hip (Finkelstein et al., 2008; Shieh et al., 2016). During this time, women lose approximately 25% of their trabecular bone (honeycomb-shaped bone structures) and 15% of their cortical bone (dense bone tissue on the outside of bone structures) (Finkelstein et al., 2008). This rapid bone degradation lasts about 5–10 years. After this period, bone is lost at a much slower rate of 0.5% per year (The North American Menopause Society, 2021).

Another cause of estrogen loss is the surgical removal of the ovaries. A bilateral oophorectomy is usually performed to reduce the risk of cancer or treat non-malignant ovarian diseases, such as endometriosis or benign cysts (Cohen et al., 2012; Challberg et al., 2011; Aitken et al., 1973). This sudden onset of estrogen deficiency increases a patient’s risk for osteoporosis (Rodriguez and Shoupe, 2015). After oophorectomy, both the age of the patient at surgery and the usage of hormone replacement therapies influence the likelihood of developing osteoporosis. Oophorectomy before age 45 has been reported to increase the risk of osteoporosis; however, those who underwent an oophorectomy after age 45 had bone density similar to that of women with intact ovaries (Aitken et al., 1973). A review of available early postmenopause data found that the earlier menopause occurred, whether natural or surgical, the lower the resulting bone density became (Gallagher, 2007). Abnormal bone scans were identified in 71% of women who underwent a preventative bilateral oophorectomy (Cohen et al., 2012). This study did not find a difference between those who underwent surgery before or after the age of natural menopause; however, the authors pointed out that there was a large difference in follow-up ages between the two groups, indicating that the age of the woman when the bone scans were taken is also an important factor that should have been accounted for (Cohen et al., 2012). Fakkert et al. (2017a) provided a systematic review of bone mineral density (BMD) following surgical menopause. They cautioned about data bias in the reported data. They concluded that while surgical menopause substantially decreases BMD, this decline becomes indistinguishable from that observed after natural menopause, once the age of natural menopause is reached (Fakkert et al., 2017a). Oophorectomy-induced bone loss may be prevented with hormone replacement therapy, but many women have an aversion to taking estrogen due to perceived risk (Challberg et al., 2011). Overall, the etiology of bone loss leading to osteoporosis due to surgical menopause needs further exploration.

The impacts of surgical menopause on mechanisms involved in bone remodeling have been explored using ovariectomized animal models, where animals undergo either an ovariectomy (OVX) procedure or a control procedure (i.e., sham operation) that mimics surgery but keeps the ovaries intact, and via in vivo studies from human biopsies or animal cells in culture. Several studies have shown that osteocyte apoptosis increased after estrogen withdrawal using in vitro experiments (Brennan et al., 2014); sheep, rat, and mouse in vivo experiments (Brennan et al., 2011; Tomkinson et al., 1998; Huber et al., 2007; Emerton et al., 2010; Florencio-Silva et al., 2018); and ex vivo experiments in human bone (Tomkinson et al., 1997). In in vivo animal studies, this increased apoptosis persisted for at least 31 months post-OVX (Brennan et al., 2011). However, the study by Florencio-Silva et al. (2018) found differing results: the number of osteocytes was significantly lower only immediately after surgery, and osteocyte counts increased by 1 month post-OVX. Osteocyte death has been proposed as a key marker of poor bone quality (Milovanovic and Busse, 2020), and the OVX experiments mentioned above suggest a marked change in osteocyte number after surgical menopause in contrast to natural menopause. Other in vitro animal studies showed that estrogen-deficient osteocytes release higher levels of RANKL (Sipos et al., 2009; Hofbauer et al., 2004), which stimulates the differentiation rate of osteoclast precursors (McNamara, 2021; Naqvi et al., 2020; Choi et al., 2008). It has not yet been quantified how sudden estrogen loss interacts with all cell types in combination to affect overall BMD and fracture risk.

Although estrogen plays a central role in bone health, few mathematical models in systems biology focus on the mechanisms of estrogen’s impact on bone (Cook et al., 2024). There is extensive literature on mathematical models of bone remodeling, as reviewed by Pivonka and Komarova (2010) and Cook et al. (2024). Most models of biochemical and cellular species dynamics use either power-law approximations (Komarova et al., 2003; Graham et al., 2013; Cook et al., 2022) or mass-action kinetics (Lemaire et al., 2004; Pivonka et al., 2008). Other models focus on how mechanical loading and morphology can affect bone remodeling (Lerebours et al., 2016; van Oers et al., 2008; Scheiner et al., 2013; 2014; Larcher and Scheiner, 2021). A recent review from our team (Cook et al., 2024) identified published mathematical models that incorporate explicit effects of estrogen in bone remodeling and postmenopausal treatments (Rattanakul et al., 2003; Schmidt et al., 2011; Post et al., 2013; Berkhout et al., 2015; 2016; Chaiya and Rattanakul, 2017; Javed et al., 2018; Jörg et al., 2022) and others with implicit estrogen effects (Lemaire et al., 2004; Scheiner et al., 2013; 2014; Lemaire and Cox, 2019; Trichilo et al., 2019; Martin et al., 2019; Larcher and Scheiner, 2021). In particular, these models investigate postmenopausal osteoporosis by adjusting parameters affected by estrogen decline or senescence. Many of these parameters describe mechanisms in the RANKL pathway (Pivonka et al., 2010; Lemaire et al., 2004; Scheiner et al., 2014; Lemaire and Cox, 2019; Martínez-Reina et al., 2021). A recent paper by Ruiz-Lozano et al. (2024), published after our review, distinguishes the effects of aging from those of estrogen decline by examining men and women. For aging, the model includes the effects of increased production of sclerostin and its impact on Wnt signaling and anabolic and catabolic effects of transforming growth factor (TGF)- $[eqn]$ , and estrogen loss is modeled implicitly via time-dependent RANKL and osteoprotegerin (OPG) changes. Our review (Cook et al., 2024) provides further details on the cells and signaling molecules included in various models, and an overview of the network of complex signaling interactions involving sclerostin, Wnt, and RANK-RANKL-OPG, along with their cellular sources. Beyond estrogen effects, models that incorporate osteocyte effects (Graham et al., 2013; Martínez-Reina et al., 2021; Cook et al., 2022; Jörg et al., 2022; Ruiz-Lozano et al., 2024) show promise in capturing age- or menopause-related changes in osteocyte signaling.

While mathematical modeling has been useful for studying postmenopausal osteoporosis, none of these models explicitly considers the effects of surgical menopause on bone cell populations. The recent model by Jörg et al. (2022) includes estrogen’s effects on osteoclasts and sclerostin, includes resorption signals (lumped effects for TGF- $[eqn]$ , bone morphogenetic protein, and the RANK-RANKL-OPG pathway, among others), and incorporates osteocyte dynamics. This model is based on realistic human time frames and includes pharmacological treatments. In particular, the mathematical model is parameterized using BMD data from patients who experienced natural menopause (Looker et al., 1998) and estimates parameters using datasets that incorporate different treatment protocols. However, the natural menopause data used to parameterize this model contains only two postmenopausal data points. Furthermore, as in the other models discussed in Cook et al. (2024), the model in Jörg et al. (2022) does not consider an abrupt decline in estrogen, which occurs in the surgical menopause scenario, nor does this model incorporate important mechanisms involved in surgical menopause, such as the impact of osteocyte death on BMD. A few prior models allow for osteocyte death (Graham et al., 2013; Cook et al., 2022), but they do not account for the effects of estrogen or aging. A limitation of most previous models is that they do not account for osteocyte apoptosis (Martínez-Reina et al., 2021; Ruiz-Lozano et al., 2024) or for osteocyte changes at all. While osteocytes are long-lived, osteocyte apoptosis is upregulated in estrogen deficiency (Tomkinson et al., 1997; 1998; Milovanovic and Busse, 2020). So, incorporating this mechanism is important for activating bone remodeling (Khosla et al., 2012), which is enhanced in surgical menopause (Peris et al., 1999; Fakkert et al., 2017b; Matsuno et al., 2025).

In this paper, we (1) aggregate BMD data sources from natural menopause patients and then reparameterize a subset of parameters in the Jörg et al. (2022) model to ensure that the mechanisms of the model reflect the broader BMD trends after natural menopause, (2) include new estrogen dynamics to describe the sudden and dramatic loss of estrogen due to surgical menopause, and (3) extend the mathematical Jörg et al. (2022) model for the dynamical responses of BMU bone cells to the case of estrogen deficiency during the surgical menopausal transition using information about the critical impacts observed in female mice and humans after removal of the ovaries. The new model considers the role of embedded osteocytes in regulating osteoclast differentiation and inducing enhanced bone resorption after surgical menopause. With this new model, we perform parameter exploration to determine which mechanisms are most important for capturing trends in surgical menopause data. This model could be used to explore medical interventions to correct the imbalances in bone remodeling after surgical menopause in a population at higher risk for early onset of osteoporosis.

Methods

2

Curated bone mineral density (BMD) data

2.1

While the impact of gradual estrogen decline on the bone remodeling process was investigated by Jörg et al. (2022), the dataset used to parameterize their model was sparse after menopause (Looker et al., 1998). In the Looker et al. (1998) dataset, only two BMD measurements from the proximal femur were recorded after menopause; furthermore, these data were aggregated by age group and did not specify menopause onset, which can lead to underestimates in BMD loss. To improve the accuracy of the Jörg et al. (2022) model in the natural menopause scenario and to parameterize our new surgical menopause model, we gather a larger set of published postmenopausal data. Due to the available datasets, we aggregate lumbar spine measurements from women without hormonal replacement for both natural and surgical menopause. In Table 1, we show data aggregated from cross-sectional studies that measured BMD using dual-energy X-ray absorptiometry in women undergoing both surgical and natural menopause.

To compare BMD measurements across datasets, each dataset’s BMD values are normalized by their respective values at menopause onset, and time is rescaled to the age at menopause onset, defined as $[eqn]$ . The rescaled time is $[eqn]$ , where $[eqn]$ is age in years. We plotted the average normalized BMD measurements and standard deviations for natural menopause (Figure 1a) and surgical menopause (Figure 1b). Motivated by linear bone loss estimates from cohort studies, we provide illustrative linear fits to data from the first 15 years postmenopause, but note that these fits were not used in our model. The slopes show a 1.54% decrease in BMD per year in the first 15 years for natural menopause (Figure 1a) and a 2.03% reduction in BMD per year over the same period for surgical menopause (Figure 1b). We note that the rate BMD loss postmenopause varies between different locations in the body with the lumbar spine having a reportedly greater loss than the femur per year (initially 1.67% and 3.12% loss per year, respectively) (Zhai et al., 2008); this is reflected in the difference between our aggregated datasets and Looker et al. (1998) (Figure 1a). Interestingly, the BMD measurements from Hadjidakis et al. (2003) suggest a rebound in BMD long after surgical menopause. We lack lumbar spine data from natural menopause due to the women’s advanced ages. The longest observational study measured femoral bone BMD up to 25 years postmenopause (Moilanen et al., 2020). The authors found that at this location in the body, bone loss during natural menopause occurs at a steady rate, with a total loss of 10% relative to baseline. Because we lack lumbar spine data for 20 years after natural menopause, it is unclear whether BMD slows or rebounds only in the lumbar spine after surgical menopause, or whether it also occurs after natural menopause.

(a,b) Bone mineral density (BMD) measured in the lumbar spine of women (except Looker et al. (1998) are from the hip). Data are normalized to premenopausal levels from the time of menopause onset, and error bars represent the standard deviation. Normalization process for each dataset is detailed in Section 2.1. The data are from Pansini et al. (1995); Ohta et al. (2002); Hadjidakis et al. (2003); Looker et al. (1998); Hibler et al. (2016); Chittacharoen et al. (1999); Yasui et al. (2007). Dashed lines show linear fits to the data with slopes listed in the legends. (a) Natural menopause and (b) surgical menopause. (c) Comparison of estrogen decline in natural and surgical menopause, shown over 30 years after menopause onset. The inset shows the decline over the 3 days following menopause onset.

Modeling bone remodeling after natural and surgical menopause

2.2

To understand how natural and surgical menopause differentially impact the bone remodeling system, we build on the mathematical model in Jörg et al. (2022), which assumes well-mixed chemical and cellular species within a BMU. Under this assumption, the dynamics of bone cell populations, chemical signals, and hormones are described by ordinary differential equations (ODEs). These species ultimately affect bone formation by altering the rates of bone production and degradation. First, we describe the model in Jörg et al. (2022), which tracks the cell densities of preosteoblasts, osteoblasts, preosteoclasts, osteoclasts, and osteocytes, as well as the sclerostin concentration, total bone density, and the bone mineral content (BMC). The precursor cells are continuously replenished and undergo apoptosis, which is influenced by various chemical signals. Osteoblasts may differentiate into osteocytes, which produce a chemical signal, sclerostin. Sclerostin upregulates osteoclast differentiation and downregulates osteoblast and osteocyte differentiation, thereby affecting bone density. The roles of the signalling molecules Wnt and RANKL are not explicitly included in the model; instead, their effects are captured by the impacts of sclerostin, and other signaling models were lumped into the resorption signals. Several estrogen effects were included in the Jörg et al. (2022) model: the inhibition of sclerostin production (Mödder et al., 2010) and the downregulation of osteoclastic bone resorption via suppressing osteoclast differentiation (Kameda et al., 1997). Therefore, a decrease in estrogen, through natural or surgical menopause, increases osteoclast and sclerostin levels and reduces osteoblast levels, which collectively lead to bone loss.

We now describe our model that extends the work of Jörg et al. (2022) to account for surgical menopause. To explore the impact of surgical menopause on bone remodeling, we incorporate two new mechanisms into the model: increased apoptosis of osteocytes through a term $[eqn]$ and increased differentiation of osteoclasts through a term $[eqn]$ . We include these effects through apoptosis and differentiation rates, which depend on the time since surgical menopause. This time dependence captures inflammation and metabolic responses outside the BMU. In Figure 2, we present a schematic of the model, illustrating how different cells and chemical signals interact and influence bone formation. The mechanisms affected by estrogen are shown in Figure 2, with new surgical menopause effects highlighted by red dashed arrows with scissor icons, representing surgical removal of the ovaries.

Schematic of mathematical model species in the bone remodeling process. The transformations of cells through differentiation, production, and degradation are illustrated by black dashed arrows. Signaling-related inhibition interactions are shown by red flat-head arrows, and activation interactions are shown as green solid arrows. Parameters that govern interactions are shown near the corresponding arrows. The schematic shows estrogen as an inhibitor of osteoclast differentiation and sclerostin production. Surgical menopause-induced changes are depicted by red dashed arrows with scissors icons. Created with BioRender.com.

Each cell population and chemical concentration is scaled by reference values for ease of computation: $[eqn]$ and $[eqn]$ are scaled by the number of preosteoclasts produced per day in the BMU; $[eqn]$ , $[eqn]$ , $[eqn]$ , and $[eqn]$ are scaled relative to the number of preosteoblasts produced per day; and estrogen is scaled by the initial concentration of estrogen at the onset of menopause.

Activation and repression signaling interactions are modeled with saturating Hill-type functions:

$$[eqn]$$

respectively. For each type of interaction by species $[eqn]$ , the threshold concentration where half of the interaction strength is achieved is the parameter $[eqn]$ , where $[eqn]$ denotes the target of the signaling interaction.

For simplicity, we assume estrogen is described by an algebraic equation, and its form depends on the type of menopause investigated. The normalized estrogen concentration over time during natural menopause is (Jörg et al., 2022):

$$[eqn]$$

where $[eqn]$ is the time of onset of estrogen decline and $[eqn]$ is the characteristic time of estrogen decline. We capture the sudden and rapid decrease in relative estrogen concentration due to oophorectomy surgery using the following equation, which is derived assuming estrogen is still synthesized at a reduced constant (zero-order) rate of $[eqn]$ post-surgery and degraded at a first-order rate with rate constant $[eqn]$ :

$$[eqn]$$

where the degradation rate of estrogen is defined as $[eqn]$ and $[eqn]$ is the half-life of estrogen, which is 161 min in postmenopausal women (Ginsburg et al., 1998). The initial concentration of estrogen before surgery is 156 pg/mL, and the estrogen level stabilizes and reaches a steady state by about 30 days post-surgery to an average value of 15 pg/mL (Bellanti et al., 2013). Thus, we set the post-surgery normalized estrogen concentration as $[eqn]$ . The synthesis rate after surgery is calculated as $[eqn]$ . Figure 1c shows the difference in the estrogen decrease in the case of natural menopause compared to surgical menopause, described by Equations 2, 3, where the inset illustrates the timescale of rapid estrogen decline in surgical menopause.

The changes in cell number over time of the preosteoclasts $[eqn]$ and osteoclasts $[eqn]$ are given by

$$[eqn]$$

and

$$[eqn]$$

respectively. Preosteoclasts are produced at a constant basal rate of one upon scaling and differentiate into osteoclasts at a rate of $[eqn]$ . The presence of estrogen inhibits osteoclast differentiation, while sclerostin activates this differentiation with thresholds $[eqn]$ and $[eqn]$ , respectively. This effect of sclerostin captures the role of RANKL implicitly. Osteoclast apoptosis occurs at a rate $[eqn]$ . Note that we do not include the apoptosis of any precursor cells (preosteoclasts or preosteoblasts) or any effect of estrogen on osteoclast apoptosis, as these were estimated to have a negligible impact in Jörg et al. (2022).

The first effect of surgical menopause is included in a new time-dependent differentiation rate for preosteoclasts to osteoclasts, $[eqn]$ , defined as

$$[eqn]$$

Before estrogen decline $[eqn]$ and for natural menopause according to the Jörg et al. (2022) model, differentiation occurs with a constant rate $[eqn]$ . At the onset of rapidly decreased estrogen due to surgical menopause (Equation 3; Figure 1c), the differentiation rate increases by $[eqn]$ , e.g., for a 10% increase, then $[eqn]$ . This increased rate of differentiation lasts for a period defined by the parameter $[eqn]$ , such that longer-lasting effects of surgery are defined by a smaller $[eqn]$ .

The preosteoblast $[eqn]$ and osteoblasts $[eqn]$ populations are governed by

$$[eqn]$$

and

$$[eqn]$$

respectively. Preosteoblasts are produced at a constant basal rate of one upon scaling and differentiate into osteoblasts at a rate of $[eqn]$ inhibited by sclerostin, capturing the role of Wnt implicitly. Osteoblasts have an apoptosis rate of $[eqn]$ and are further differentiated into osteocytes at a rate of $[eqn]$ .

The dynamic population of osteocytes $[eqn]$ is described by

$$[eqn]$$

Osteocytes are derived from osteoblasts and are removed at a rate dependent on $[eqn]$ .

The second effect of surgical menopause is increased apoptosis of osteocytes via the new time-dependent rate, $[eqn]$ , defined as

$$[eqn]$$

Similarly to Equation 6, we assume that before menopause and for natural menopause according to the Jörg et al. (2022) model, differentiation occurs at a constant rate of $[eqn]$ . At surgical menopause onset, apoptosis increases by $[eqn]$ , then returns to previous levels over a timescale of $[eqn]$ . This timescale is assumed to be equal to the timescale of increased osteoclast differentiation in Equation 6.

The production of the signaling molecule sclerostin is governed by

$$[eqn]$$

where sclerostin is produced by osteocytes at a rate inhibited by estrogen with threshold $[eqn]$ and is degraded at rate $[eqn]$ . Sclerostin affects bone formation through activation of osteoclastogenesis in Equations 4, 5 and inhibition of osteoblastogenesis in Equations 7, 8, respectively.

Bone density, $[eqn]$ , is determined by

$$[eqn]$$

where the resorption factor $[eqn]$ is assumed to be equal to the amount of osteoclasts present, i.e., $[eqn]$ . Osteoclasts inhibit bone formation, and osteoblasts contribute to bone formation. From Equation 12, the rate of bone density change increases proportionally to osteoblast number at a rate of $[eqn]$ . This production is modulated by sclerostin inhibition with a threshold of $[eqn]$ and resorption factor $[eqn]$ activation with a strength of $[eqn]$ and a threshold of $[eqn]$ . Bone density decreases through bone resorption by osteoclasts at a rate $[eqn]$ . To determine BMD, we scale bone density by bone mineral content, which is a constant $[eqn]$ , so we have $[eqn]$ .

To summarize, our model extensions introduce three new parameters fit to data: the factor for peak increase in osteocyte apoptosis due to surgery $[eqn]$ , the factor for peak increase in osteoclast differentiation due to surgery $[eqn]$ , and the timescale during which these effects last $[eqn]$ . The model parameters and species definitions for natural menopause are listed in Table 3, the initial conditions for the variables are listed in Table 2, and the model parameters and species definitions for surgical menopause are listed in Table 4.

Model solution and parameter estimation

2.3

We are interested in understanding and parameterizing the impacts of natural and surgical menopause on the dynamics of bone mineral density, $[eqn]$ . Therefore, we solve the model ODEs and algebraic equations defined in Equations 1–12 using ode45 in MATLAB, with absolute and relative tolerances as $[eqn]$ , from 30 years before menopause onset to 30 years after menopause onset. Thus, the initial condition is 30 years before menopause onset, whether natural or surgical. For each simulation, we initialize the dynamic species based on their steady-state values at the initial condition. The Jörg et al. (2022) model includes premenopausal BMD decline. To determine the steady state solutions to Equations 4, 5, 7–9, 11 while $[eqn]$ is fixed at a value of 1, we solve the system of equations with time derivatives set to 0 using the fsolve function in MATLAB. The other species do not depend on $[eqn]$ . This initialization is called within the parameter estimation routine for the natural menopause case, as the fitted parameters affect the initial conditions. These updated values are reported in Table 2.

To estimate parameters in both the natural and surgical menopause mechanisms, we use the lsqnonlin function in MATLAB, which solves the nonlinear least-square objective function in Equation 13 using the Levenberg-Marquardt algorithm (Levenberg, 1944; Marquardt, 1963):

$$[eqn]$$

Here, $[eqn]$ is the natural BMD data at measurement times $[eqn]$ , and $[eqn]$ is the BMD predicted by the model at the same times using parameters $[eqn]$ . For each parameter estimation procedure, lsqnonlin algorithm options are set to tolerances of $[eqn]$ and a maximum of 10,000 function evaluations and iterations. We use MATLAB version R2024b in Windows on a PC with 11th Generation Intel Core i7-11700 8 Core processor. We have provided model code and files in a repository at https://github.com/ashleefv/ SurgicalMenopauseBone (Nelson et al., 2025b).

With the natural menopause BMD data described in Section 2.1, we first aim to re-estimate model parameters related to the impact of estrogen on osteoclastogenesis and the production of sclerostin by osteocytes for natural menopause. This corresponds to parameters $[eqn]$ , and the estimated values of these parameters are listed in Table 3. With the estimated $[eqn]$ , we then use surgical menopause BMD data to estimate surgical menopause parameters $[eqn]$ in the 15 years (short-term) and in the 30 years (long-term) after onset of menopause due to surgery. For the surgical menopause cases, upper and lower bounds are used to constrain the parameter space. The lower bounds are $[eqn]$ , which signify no effect, permanent effect, and no effect, respectively. The upper bounds are more subjective and were selected to yield only reasonable responses; we would need cell population data to further interrogate these with different upper bounds. The upper bounds used are $[eqn]$ , signifying that the peaks of $[eqn]$ and $[eqn]$ at menopause onset are 5 + 1 = 6 times higher than baseline values in natural menopause and the surgery induced effects last a minimum of 1 day. The estimated values of these parameters are listed in Table 4.

We calculate sensitivity by varying the surgical menopause model parameters $[eqn]$ by 25% relative to the best-fit parameters. For instance, the upper bound is calculated using a 25% increase in $[eqn]$ and $[eqn]$ and a 25% decrease in $[eqn]$ , which corresponds to longer-lasting effects of surgery.

We determine the BMD sensitivity to sclerostin levels presented in Section 3.3 by calculating the steady state levels of cells and chemical concentrations, defined by Equations 4–11. We then calculate the constant rates of BMD loss or production using Equation 12 with the steady-state values of osteocytes, sclerostin, and osteoblasts. Estrogen is assumed to be at a fixed value either at premenopause levels of $[eqn]$ or at the post-surgical menopause levels of $[eqn]$ . We calculate the curve of steady BMD change rates in Figure 4 using the model with no new effects, namely, with $[eqn]$ . We perturb the sclerostin production rate by adding a control parameter $[eqn]$ multiplying sclerostin production in Equation 11 so that its steady state is defined by $[eqn]$ . The curve is defined by $[eqn]$ with $[eqn]$ . The reference points for the long and short term model fits in Figure 4, are defined by $[eqn]$ and $[eqn]$ , respectively, reflecting that the osteocyte levels in each case obtain minimum values of 80% and 94% of the premenopause levels in our results (Figure 3).

(a) Reparameterization of the Jörg et al. (2022) model using more sources of natural menopause (NM) data. The dashed red curve shows the original model, the solid red curve shows the new parameterization, and the black solid curve shows the model with sudden estrogen loss alone, without any new effects. (b) Extension of the model to surgical menopause. New effects are parameterized using long-term (blue curve) and short-term (green curve) data. The shaded regions highlight model sensitivity to parameters. (c) Cell populations from the surgical menopause model using the three parameter sets shown in panel (b). NM: Natural menopause. SM: surgical menopause.

Results

3

Using the mathematical framework introduced in Section 2.2, we study how mechanisms mediated by estrogen loss affect bone remodeling in both natural and surgical menopause. We focus on how BMD is affected decades after menopause and compare the BMD dynamics to both natural and surgical menopause data described in Section 2.1. We use data from natural and surgical menopause patients (Figure 1) and the mathematical models outlined in Section 2.2 with the methods described in Section 2.3 to parameterize the models and compare their results with the relevant data. Using new parameters and mechanisms, we identify key pathways that drive BMD decline and rebound and propose new treatment directions based on our results.

Parameterized model of natural menopause captures BMD behavior in larger lumbar spine dataset

3.1

To assess the model output, we compare the model-predicted BMD with the experimental data described in Section 2.1. In Figure 3a for the natural menopause case, we present the average BMD for natural menopause patients in red markers for the larger dataset and distinguish the femur BMD data (Looker et al., 1998) that were used to parameterize the Jörg et al. (2022) model mathematical model with open red markers. The Jörg et al. (2022) model results (red dashed curve) fit the Looker et al. (1998) data well but do not fit the additional natural menopause data measured in the lumbar spine, which show a faster BMD decline.

To improve the model’s accuracy in predicting lumbar spine BMD dynamics after menopause, we estimate key parameters associated with natural menopause using aggregated natural menopause data. In particular, we estimate functional thresholds that depend on estrogen in both preosteoclast to osteoclast differentiation and sclerostin production from osteocytes, which are $[eqn]$ and $[eqn]$ , respectively. Our reparameterized model of natural menopause (shown as a red solid curve) captures this rapid BMD decrease, and this behavior is within the experimental error (Figure 3a). We terminate natural menopause simulations at approximately 20 years post-menopause onset, matching the time period of our dataset.

With new mechanisms, model of surgical menopause reproduces BMD trends in both short- and long-term data

3.2

Figure 3b shows the average relative BMD of surgical menopause patients 15 years post-surgery (black markers) and long-term data beyond 15 years (blue markers), indicating a rebound in BMD. Note that these are the same data in Figure 1b, where the error bars are shown.

To illustrate that new mechanisms are necessary to capture the surgical menopause behavior, we first show the model dynamics (Figures 3a,b black curve) using newly estimated bone parameters for natural menopause from Section 3.1 and the estrogen decline in surgical menopause from Equation 3 without including any new effects on cell dynamics, i.e., $[eqn]$ and $[eqn]$ . The effect of sudden estrogen loss alone enhances BMD loss compared to natural estrogen decline in Figure 3a. Indeed, the model with no new cell dynamics closely matches the surgical menopause data for the first 1–3 years post-surgery but overpredicts the extent of BMD loss in later years (Figures 3a,b black curve). This result suggests that surgical menopause involves both an increased rate of bone loss in the short term and a subsequent slowing of bone loss in the long term, which is not captured by simply including the onset of a sudden estrogen decline in the natural menopause model of Jörg et al. (2022).

We estimate the new parameters modulated by surgical menopause in our model (Table 4): the percentage increase osteocyte apoptosis rate $[eqn]$ , the increased differentiation rate of preosteoclasts $[eqn]$ , and the timescale over which the effects occur $[eqn]$ , using data in Figure 1b. We obtain root mean squared errors in BMD prediction of 4.59% and 4.63% compared to the surgical menopause data up to 15 years and 30 years, respectively. Using the short and long time scale datasets resulted in different long-term BMD dynamics (Figure 3b). The short-term surgical menopause model better captures the data for the 15-year period compared to the case without new effects (Figures 3a,b black curve). The long-term surgical menopause model also fits these data well and has a BMD rebound not observed in the short-term case. The sensitivity of the model predictions to the parameter values is shown by the corresponding shaded area around each curve. Across these sensitivity regions, the model’s behavior over the 2 years post-surgery is insensitive to the new effects, whereas parameter variations substantially alter its long-term predictions.

To understand other differences in these parameterized surgical menopause models, we show the dynamics of the osteoclast, osteoblast, and osteocyte cell populations as well as the levels of sclerostin after menopause (Figure 3c). The no-new-effects case and the short-term case show that all species (except sclerostin) return to approximately their premenopausal levels within 10 years. This is expected because the timescale of surgical effects, $[eqn]$ , is small for the short-term model. The short-term effect is observed in the osteoclast (osteocyte) population, with a sharp increase (decrease) at the onset of menopause, followed by a rebound. The sudden increase in osteoclast population yields steeper dips in Figure 3b, i.e., relative BMD in the short-term surgical menopause model is lower than other models immediately after menopause onset.

In the long-term surgical menopause case, we see a smaller increase in osteoclastogenesis after menopause onset compared to the short-term case; $[eqn]$ is an order of magnitude smaller in the long-term case than the short-term case. Osteocyte density also decreases more slowly in the long-term surgical menopause model, but continues after the onset of menopause; $[eqn]$ for long-term surgical menopause effects, making the impact of surgery permanent. Interestingly, the percentage of osteoblasts decreases after surgery, then rebounds to a slightly higher level for an extended period.

Osteocyte and sclerostin dynamics are key drivers of slowing BMD loss

3.3

In our model, the signalling molecule sclerostin is produced solely by osteocytes, with its production rate increasing as estrogen levels fall. Sclerostin levels play a crucial role in regulating bone formation and resorption. Figure 4 shows the sensitivity of BMD production rate as a function of the sclerostin production rate (method detailed in Section 2.3). Higher sclerostin production rates lead to increased bone loss, whereas sclerostin production rates below 96% of the premenopause level lead to bone formation. In the absence of additional modeling effects, the production rate of sclerostin postmenopause is 10% higher than premenopause (Figure 3c). This leads to a bone loss of −0.75% per year postmenopause compared to −0.22% per year premenopause (marked on Figure 4).

Steady BMD change as a function of sclerostin production levels at steady state (relative to the premenopause level). Steady-state cell concentrations are calculated assuming that no new effects are included in the model. Sclerostin production is then adjusted by a percentage of its premenopause level to calculate the resulting steady-state change in BMD.

Our model extension, which simulates a surgically induced temporary increase in osteocyte apoptosis, reduces sclerostin production, slows bone loss, and captures the reduced bone loss or rebound observed in clinical data. The long-term surgical menopause model yields a continually declining osteocyte population with osteocyte levels reaching $[eqn]$ 80% of their premenopause levels 25 years post-surgery. In the short-term model, the osteocyte population reaches $[eqn]$ 94% of its premenopause levels before slowly returning to $[eqn]$ . If an 80% or 94% lower level of osteocytes persisted at steady state, this would yield BMD growth of 1.2% or 0.15% per year, respectively, as marked on Figure 4. Although these osteocyte levels are not obtained at steady state because the bone system is still responding to surgery, this explains how these osteocyte values alter sclerostin production and, subsequently, BMD, resulting in a reduced rate of bone loss or even a rebound in BMD.

Discussion

4

We present a mathematical model of bone remodeling to quantify the effects of estrogen loss in surgical menopause. Since experimental data suggest that surgical and natural menopause affect different mechanisms of bone remodeling, we extended an existing mathematical model of bone remodeling to incorporate increased osteoclast differentiation and osteocyte apoptosis. The objective of this framework is to understand and capture trends seen in newly aggregated physiological data: (1) surgical menopause leads to an increased loss of BMD in the short term, and (2) this loss slows or even rebounds by 10 or more years post-surgery.

The BMD predictions after reparameterization to the larger natural menopause BMD dataset better capture the overall trends of the natural menopause data compared to the previous mathematical model (Jörg et al., 2022). These parameters, $[eqn]$ and $[eqn]$ , influence the strength of estrogen signaling on osteoclastogenesis and sclerostin release by osteocytes. Compared to the previous mathematical model, our parameter fitting results in smaller threshold values for the estrogen signaling pathway related to osteoclastogenesis, resulting in more osteoclast differentiation for a higher concentration of estrogen. These parameter changes capture the overall larger decrease in BMD shown in the newly aggregated data.

Our new surgical menopause model that incorporates mechanisms impacted by the sudden loss of estrogen and inflammation resulting from surgery, including increases in osteocyte apoptosis and osteoclastogenesis follows the sharp decrease in BMD in the first 15 years post-surgery and a rebound in BMD 15–30 years post-surgery, consistent with the clinical data. To understand which mechanisms underlie this varied behavior, we fit our surgical menopause model to two datasets: BMD data from up to 15 years and up to 30 years post-surgery. The short-term data fit indicates a higher osteocyte apoptosis rate and a lower osteoclast differentiation rate compared to the long-term data fit. Osteocyte levels differ substantially between the long-term and short-term model calibrations, with the long-term parameters yielding slower but permanent rates of osteocyte apoptosis due to surgical menopause. These results show that the model is capable of fitting the available data, but long-term predictions should be made after the model is trained on long-term data.

Our mathematical study does have some limitations. Here, we do not consider the mechanical mechanisms underlying the initiation and regulation of bone remodeling, which are mediated by osteocytes’ response to mechanical strain (Bonewald and Johnson, 2008; Santos et al., 2009). However, sclerostin levels have been shown to decrease with mechanical loading (Robling et al., 2008), so another potential mechanism for decreasing sclerostin could be to incorporate strain into our mathematical framework. We model bone remodeling within an individual BMU and do not account for strain-induced deformation of the bone matrix. This work also does not incorporate hormonal interventions or bone remodeling treatments. Most osteoporosis treatments are classified as anti-resorptive agents that inhibit osteoclast-mediated resorption to prevent further bone loss, but they can cause side effects. For example, bisphosphonates are widely prescribed antiresorptive treatments to promote osteoclast apoptosis (Rodan and Fleisch, 1996; Rogers et al., 1999; Russell et al., 1999; Berkhout et al., 2015), but they can lose efficacy over time and lead to osteonecrosis and atypical femoral fractures (Shane et al., 2014; Khosla et al., 2007; Whitaker et al., 2012). Another anti-resorptive treatment is the anti-RANKL monoclonal antibody denosumab, which inhibits osteoclastogenesis (Khosla and Hofbauer, 2017). Other antibody treatments targeting sclerostin produced from osteocytes show promising results of increased bone formation and decreased bone resorption in rodent models (Li et al., 2009; Zhang et al., 2016; Tian et al., 2011) and in postmenopausal women with bone loss (McClung et al., 2014; Recker et al., 2015; Cosman et al., 2016; Zhang et al., 2016; Chavassieux et al., 2019; McClung, 2017). However, bone formation only lasts a few months (Ominsky et al., 2017). Therefore, treatment with sclerostin inhibitors is only recommended for 6–12 months (McClung, 2017). Our parameter estimates indicate promising avenues for treatment in patients undergoing surgical menopause, suggesting that targeting osteocyte or osteoblast dynamics may support long-term BMD preservation. Implementing treatment methods on this system requires further investigation.