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 t m . The rescaled time is t − t m , where t 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.

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 η ( t ) and increased differentiation of osteoclasts through a term ω ( t ) . 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.

Each cell population and chemical concentration is scaled by reference values for ease of computation: P C and C are scaled by the number of preosteoclasts produced per day in the BMU; P B , B , S , and S c 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: f + X , x i = X X + x i , f − X , x i = x i X + x i , (1) respectively. For each type of interaction by species X , the threshold concentration where half of the interaction strength is achieved is the parameter x i , where i 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): E nat t = 1 , t ≤ t m 1 1 + t − t m / τ E , t > t m , (2) where t m is the time of onset of estrogen decline and τ E 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 k syn post-surgery and degraded at a first-order rate with rate constant κ E : E surg t = 1 , t ≤ t m 1 − k syn κ E exp − κ E t − t m + k syn κ E , t > t m , (3) where the degradation rate of estrogen is defined as κ E = ln ( 2 ) / t 1 / 2 and t 1 / 2 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 E surg ( t m + 30 d a y s ) = 15  pg/mL 156  pg/mL . The synthesis rate after surgery is calculated as k syn = κ E E surg ( t m + 30 d a y s ) . 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 ( P C ) and osteoclasts ( C ) are given by d P C d t = 1 − ω t f − E , e PC f + S c , s c PC P C (4) and d C d t = ω t f − E , e PC f + S c , s c PC P C − η C C , (5) respectively. Preosteoclasts are produced at a constant basal rate of one upon scaling and differentiate into osteoclasts at a rate of ω ( t ) . The presence of estrogen inhibits osteoclast differentiation, while sclerostin activates this differentiation with thresholds e PC and s c PC , respectively. This effect of sclerostin captures the role of RANKL implicitly. Osteoclast apoptosis occurs at a rate η C . 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, ω ( t ) , defined as ω t = ω PC , t ≤ t m , ω PC 1 + ω surg ⁡ exp − τ t − t m , t > t m . (6) Before estrogen decline t ≤ t m and for natural menopause according to the Jörg et al. (2022) model, differentiation occurs with a constant rate ω PC . At the onset of rapidly decreased estrogen due to surgical menopause (Equation 3; Figure 1c), the differentiation rate increases by ω surg , e.g., for a 10% increase, then ω surg = 0.1 . This increased rate of differentiation lasts for a period defined by the parameter τ , such that longer-lasting effects of surgery are defined by a smaller τ .

The preosteoblast ( P B ) and osteoblasts ( B ) populations are governed by d P B d t = 1 − ω PB f − S c , s c PB P B (7) and d B d t = ω PB f − S c , s c PB P B − η B + ω B B , (8) respectively. Preosteoblasts are produced at a constant basal rate of one upon scaling and differentiate into osteoblasts at a rate of ω PB inhibited by sclerostin, capturing the role of Wnt implicitly. Osteoblasts have an apoptosis rate of η B and are further differentiated into osteocytes at a rate of ω B .

The dynamic population of osteocytes ( S ) is described by d S d t = ω B B − η t S . (9) Osteocytes are derived from osteoblasts and are removed at a rate dependent on η ( t ) .

The second effect of surgical menopause is increased apoptosis of osteocytes via the new time-dependent rate, η ( t ) , defined as η t = η S , t ≤ t m , η S 1 + η surg ⁡ exp − τ t − t m , t > t m . (10) 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 η S . At surgical menopause onset, apoptosis increases by η surg , then returns to previous levels over a timescale of τ . 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 d S c d t = f − E , e Sc S − κ Sc S c . (11) where sclerostin is produced by osteocytes at a rate inhibited by estrogen with threshold e Sc and is degraded at rate κ Sc . Sclerostin affects bone formation through activation of osteoclastogenesis in Equations 4, 5 and inhibition of osteoblastogenesis in Equations 7, 8, respectively.

Bone density, B d , is determined by d B d d t = λ B B f − S c , s c Ω 1 + ν Ω f + R , r Ω − λ C C , (12) where the resorption factor R is assumed to be equal to the amount of osteoclasts present, i.e., R = C . 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 λ B . This production is modulated by sclerostin inhibition with a threshold of s c Ω and resorption factor R activation with a strength of ν Ω and a threshold of r Ω . Bone density decreases through bone resorption by osteoclasts at a rate λ C . To determine BMD, we scale bone density by bone mineral content, which is a constant B M C 0 , so we have B M D = B M C 0 B d .

To summarize, our model extensions introduce three new parameters fit to data: the factor for peak increase in osteocyte apoptosis due to surgery ( η surg ) , the factor for peak increase in osteoclast differentiation due to surgery ( ω surg ) , and the timescale during which these effects last ( τ ) . 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.

We are interested in understanding and parameterizing the impacts of natural and surgical menopause on the dynamics of bone mineral density, B d . Therefore, we solve the model ODEs and algebraic equations defined in Equations 1–12 using ode45 in MATLAB, with absolute and relative tolerances as 1 0 − 8 , 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 B d 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 B d . 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): arg min q ∑ j B M D data t j − B M D t j , q 2 . (13) Here, B M D data ( t j ) is the natural BMD data at measurement times t j , and B M D ( t j , q ) is the BMD predicted by the model at the same times using parameters q . For each parameter estimation procedure, lsqnonlin algorithm options are set to tolerances of 1 0 − 8 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 q NM = { e PC , e Sc } , and the estimated values of these parameters are listed in Table 3. With the estimated q NM , we then use surgical menopause BMD data to estimate surgical menopause parameters q SM = { η surg , τ , ω surg } 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 q SM = { 0,0,0 } , 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 q SM = { 5,1 / 365,5 } , signifying that the peaks of η ( t ) and ω ( t ) 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 q SM = { η surg , τ , ω surg } by 25% relative to the best-fit parameters. For instance, the upper bound is calculated using a 25% increase in η surg and ω surg and a 25% decrease in τ , 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 E = 1 or at the post-surgical menopause levels of E = k syn / κ E . We calculate the curve of steady BMD change rates in Figure 4 using the model with no new effects, namely, with η surg = ω surg = 0 . We perturb the sclerostin production rate by adding a control parameter α multiplying sclerostin production in Equation 11 so that its steady state is defined by S c = α f − ( E , e Sc ) S / κ Sc . The curve is defined by α ∈ [ 0.7 , 1.2 ] with E = k syn / κ E . The reference points for the long and short term model fits in Figure 4, are defined by α = 0.8 and α = 0.94 , 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).

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 e PC and e Sc , 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.

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., η surg = 0 and ω surg = 0 . 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 ( η surg ) , the increased differentiation rate of preosteoclasts ( ω surg ) , and the timescale over which the effects occur ( τ ) , 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, τ , 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; ω surg 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; τ = 0 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.

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).

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 ≈ 80% of their premenopause levels 25 years post-surgery. In the short-term model, the osteocyte population reaches ≈ 94% of its premenopause levels before slowly returning to 100 % . 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.