Predicting mechanoregulatory responses in bone during breast cancer metastasis: A Finite Element Analysis

doi:10.21203/rs.3.rs-3141613/v1

Download PDF

Research Article

Predicting mechanoregulatory responses in bone during breast cancer metastasis: A Finite Element Analysis

https://doi.org/10.21203/rs.3.rs-3141613/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Breast cancer metastasises to bone in 70–80% of patients with advanced disease. Bone cells contribute to tumour metastasis by activating bone resorption, which releases biochemical factors that stimulate tumour cell proliferation. The local mechanical environment of bone tissue is altered during early metastasis, prior to the formation of overt osteolytic metastasis. According to mechanoregulation theory, these changes might activate mechanobiological responses in bone cells and thereby contribute to osteolytic resorption. However, whether mechanobiological responses of bone cells drive osteolysis during metastasis is unknown. The objective of this study was to apply a computational mechanoregulation framework to predict how early changes in the bone mechanical environment contribute to osteolysis. Subject-specific finite element (FE) models were developed to predict the mechanical environment within bone tissue during early stage metastasis (3 weeks post-inoculation). We then applied a mechanoregulation algorithm to predict changes in bone tissue density as a function of the evolving mechanical environment due to tumour invasion. Substantial bone loss was predicted in the greater trochanter region, which coincides with experimental reports of regional bone loss in this animal model. Moreover, application of the mechanoregulation algorithm predicted that the mechanical environment evolved in a similar manner to that predicted through subject-specific finite element (FE) models. This is the first study to implement a computational mechanoregulation framework to predict the development of osteolysis. Our findings support the hypothesis that early changes in the physical environment of bone tissue during metastasis may elicit mechanobiological cues for bone cells and activate osteolytic destruction.

Bone

metastasis

finite element analysis

mechanobiology

computational modelling

Breast cancer commonly migrates from the primary tumour and metastasises to bone, and it has been proposed that breast cancer cells have a particular affinity to the physical properties of bone tissue, which provides an attractive environment for metastatic development (Paget 1889). Bone cells are instrumental in tumour metastasis, whereby osteoclasts activate bone tissue resorption, resulting in osteolysis. Osteolysis not only provides a space for tumour invasion but also releases biochemical factors from the bone matrix, which stimulates tumour cell proliferation.

In healthy bone, mechanosensitive osteocyte cells, which reside within the bone lacunae, release biochemical signals that signal to osteoblasts and osteoclasts to activate bone formation or bone resorption, depending on the mechanical demands on the tissue. Frost’s ‘Mechanostat’ theory proposed that bone cells remodel tissue in response to changes in the tissue strain. Specifically, it was proposed that there is a homeostatic strain range (‘lazy zone’), within which no net bone formation or resorption occurs, whereas above this threshold induces bone formation and bone resorption is activated below the lower strain threshold (Frost 1987). As bone tumours evolve, the local mechanical environment of bone tissue is altered (Arrington et al. 2006, Nazarian et al. 2008, Richert et al. 2015, Verbruggen et al. 2022). A recent pre-clinical animal study of metastatic mouse femurs has revealed the temporal nature of such changes, whereby bone mineral content and cortical thickness was significantly altered prior to the development of overt osteolytic lesions by 3 weeks of tumour development (Verbruggen et al. 2022), whereas osteolysis and significantly decreased bone mineral content were not reported until the later time point of 6 weeks post-inoculation (Verbruggen et al. 2022). In a following study, micro-CT-derived finite element (FE) models of the proximal femurs of these animal cohorts were developed and applied to study the mechanical environment within bone tissue during bone metastasis and osteolytic resorption (Verbruggen and McNamara 2023). Interestingly, in early metastasis there was a decrease in strain distribution within the bone tissue of the proximal femur, which coincided with the onset of early osteolysis, cortical thickening and mineralisation of proximal and distal femur bone (Verbruggen et al. 2022). It was proposed that this altered mechanical environment might stimulate resident osteocytes to activate osteoclast resorption, and thereby play a role in the cancer vicious cycle. However, how such an altered mechanical environment could drive the remodelling behaviour of resident bone cells remains unclear.

Mathematical and computational models have been applied to simulate and predict the role of biochemical processes and cell-cell signalling in tumour growth, and to a lesser extent have considered biophysical processes (Araujo et al. 2014, Cook et al. 2016, Kim and Othmer 2015, Rejniak and Anderson 2011, Tracqui 2009, Zhou and Liu 2014). The mechanoregulatory behaviour of bone tissue has been investigated through a combination of FE models and mechanoregulation theory. The basis of mechanoregulation theory is that the tissue is adapted to reach a homeostatic mechanical environment, and so mathematical relationships between the tissue strain and adaptation have been developed. A bone remodelling algorithm incorporating strain and microdamage stimuli was developed and applied to study trabecular remodelling in healthy bone and during osteoporosis (McNamara and Prendergast 2007, Mulvihill et al. 2008). A microstructural, load-dependent remodelling algorithm for trabecular bone adaptation, which was successfully validated by in-vivo micro-CT data, was capable of predicting formation and resorption of trabeculae in both healthy (Schulte et al. 2011) and ovariectomised (OVX) vertebrae (Schulte et al. 2013). Later, a mechanoregulation algorithm, driven by fluid velocity, accurately predicted cortical thickening in healthy cortical bone tissue in response to axial compression observed in experimental data (Pereira et al. 2015). A strain-dependent mechanoregulation algorithm was implemented to predict bone adaptation in a poroelastic FE models of a human femur (Villette et al. 2020), and was applied to investigate the role of femoral morphology for fracture behaviour. A more recent study sought to predict cortical and trabecular bone adaptation, regulated by strain energy density (SED) or maximum principal strain, in finite element models derived from longitudinal micro-CT scans (every 2 weeks) in mice (Cheong et al. 2020a). Interestingly, predictions of bone apposition were more accurate using maximum principal strain as a stimulus, whereas bone resorption was better predicted when driven by strain energy density (Cheong et al. 2020a). Another study applied a SED-driven remodelling algorithm to µCT-derived FE models of mouse tibiae, to investigate bone adaptation following ovariectomy and subsequent mechanical loading treatments (Cheong et al. 2020b). This approach was able to predict regions of bone apposition that were well correlated to experimental findings, whereas bone tissue resorption was not as precisely or spatially matched to experimental data (Cheong et al. 2020b). A recent study predicted bone resorption and formation, driven by microdamage and strain energy density respectively, in a FE model of human tibial bone fracture repair (Quinn et al. 2022). Machine learning has also been applied to analyse micro-CT imaging and changes in bone tissue stiffness within the trabeculae of osteoporotic patients to determine the key variables that influenced bone tissue resorption (yield stress, yield strain) (Sohail et al. 2019). These parameters were then implemented in an algorithm, whereby bone tissue adaption or resorption were governed by a strain stimulus which accurately predict changes in tissue morphology in a single 2D trabecular strut model (Sohail et al. 2019). Taken together, it has been established that computational modelling can predict mechanoregulatory responses and associated bone adaptation for healthy and osteoporotic bone models. However, mechanoregulation theory has not been applied to investigate the contribution of bone remodelling during bone metastasis. Applying such approaches may shed light on the potential role of mechanobiology in the development of osteolysis as bone metastasis progresses.

The objective of this study was to test the hypothesis that an altered physical environment of bone tissue could elicit mechanobiological cues for resident bone cells during metastasis, and thereby contribute to later osteolytic destruction. We implemented (a) subject-specific FE models to predict the mechanical environment within the bone tissue at the early stage of metastasis and (b) a mechanoregulation theory based on strain energy density to predict how material properties are altered as a function of the evolving mechanical environment within bone tissue during breast cancer-bone metastasis.

In this study an iterative approach is taken, whereby subject-specific FE models are implemented to predict the mechanical environment within the bone tissue at the early stage of metastasis. A mechanoregulation approach, which adapts the bone tissue properties on the basis of the mechanical environment (strain energy density), is applied to predict how material properties are altered as a function of the evolving mechanical environment within the bone tissue.

2.1 Obtaining bone tissue geometry and material properties

The current study builds upon µCT-derived FE models, developed in a previous computational study of strain distribution during breast cancer metastasis (Verbruggen and McNamara 2023). In brief, immune competent BALB/c mice, inoculated into a femur-adjacent mammary fat pad with 4T1 breast cancer cells, were euthanised, at either 3 or 6 weeks post-inoculation, and µCT scanning and reconstruction was conducted to develop solid models of the proximal femora.

In the current study, computational models from the 3-week murine cohort were developed by meshing with 4-noded linear elastic, homogeneous, tetrahedral elements, following mesh convergence analysis, using 3Matic and Abaqus software (Fig. 1A, 1B). A study compared 4- and 10-noded tetrahedral meshes in a human proximal femur FE model under compression and reported no significant differences in their output stress, strain or overall accuracy compared to experimental results (Ramos and Simoes 2006). Loading was applied to the femoral head surface equating to 120% of bodyweight, recorded for each animal during the experimental study (Verbruggen et al. 2022) and applied in the proximal-distal direction and 10.9% in the posterior-anterior direction (Charles et al. 2018) (see Fig. 1C). The distal surface was fixed in all directions (Fig. 1C). For the current study, each computational model was assumed an initially homogeneous, isotropic structure, with the mean bone mineral density assigned according to µCT imaging of the metastatic femoral bone (1.262 ± 0.02 g/cm³), converted to ash density (Verbruggen et al. 2022) and the corresponding Young’s modulus applied according to a power law relationship used in female mouse tibiae, according to the following equation (Lu et al. 2019).

$$E=10.5{\rho }^{2.29}$$

2.2 Mechanoregulation Theory

To introduce mechanoregulation theory, a USDFLD (user-defined field) subroutine was applied, whereby the density and Young’s modulus of each element within the model mesh was incrementally updated according to the rate of change of density predicted by a remodelling algorithm. Specifically, density adaption was predicted on the assumption that a mechanoregulatory response occurs, driven by a mechanical stimulus, strain energy density, U ^j (Huiskes et al. 1987, McNamara and Prendergast 2007), calculated according to:

$${U}^{j}=\frac{{{{E}^{j}(\epsilon }^{j})}^{2}}{2{\rho }^{j}}$$

where, for each location within a mesh, j, the tissue Young’s modulus, E ^j (MPa) and density, ρ ^j (g/cm³) are defined and the strain (ɛ ^j) is predicted by solving the finite element model. The stimulus, S_strain, is calculated by comparing to the reference strain energy density, U^ref, which is calculated at homeostatic equilibrium, according to the following equation:

$${\text{S}}_{strain}=\sum _{i=1}^{n}{f}_{i}\left(x\right)\left\{{U}_{i}\left(t\right)+{U}^{ref}\right\}$$

The homeostatic region was defined according to Frost’s ‘Mechanostat’ theory (1987), adapted for bone tissue mechanoregulation. Reference strain energy density, U^ref_, is calculated as a function of reference modulus E^ref, density ρ^ref and strain ɛ^ref, when bone remodelling is at equilibrium (see Eq. 4). This ‘lazy zone’ of equilibrium is defined by strain thresholds for net bone resorption and formation are ɛ^R and ɛ^F, respectively. The value of ɛ^ref is thus determined by the current value of ɛ ^j, according to conditions outlined below:

$${U}^{ref}=\frac{{E}^{ref}{\left({\epsilon }^{ref}\right)}^{2}}{2{\rho }^{ref}}$$

$$\left\{\begin{array}{c}\psi =-1, {{\epsilon }}^{ref} = {{\epsilon }}^{R} if 0< {{\epsilon }}^{\text{j}}< {{\epsilon }}^{R}\\ \\ \psi =0, {{\epsilon }}^{ref} = {{\epsilon }}^{j} if {{\epsilon }}^{R}< {{\epsilon }}^{\text{j}}< {{\epsilon }}^{F}\\ \\ \psi =1, {{\epsilon }}^{ref} = {{\epsilon }}^{F} if {{\epsilon }}^{\text{j}}> {{\epsilon }}^{F} \end{array}\right\}$$

where ψ is the mechanical signal ensuring density decreases under resorption conditions and increases upon bone formation. Finally, the rate of change of density in response to mechanical stimuli, dρi/dt, and a new resulting density, (ρ_i+1), are calculated using the following equations:

$$\frac{d{\rho }_{i}}{dt}= {C}_{1}\bullet {\psi S}_{strain}$$

$${\rho }_{i+1}= {\rho }_{i}+\frac{d{\rho }_{i}}{dt}$$

where a time constant, C₁, governs the rate of remodelling. This constant was determined separately for bone tissue resorption (C₁^R = 1600) and formation (C₁^F = 1.79) as described below. Maximal and minimal apparent density limits of 1.73 g/cm³ and 0.01 g/cm³ were assigned as standard, where density below 0.01 g/cm³ represents complete resorption (Scannell and Prendergast 2009, Van Rietbergen et al. 1993). These equations (5.1–5.5) are thus implemented sequentially as part of the remodelling algorithm to iteratively predict a new density and Young’s modulus for each element within a given model mesh, as illustrated in a flowchart (Fig. 1D).

2.3 Determining bone remodelling thresholds and constants

An upper strain threshold, beyond which murine bone tissue undergoes formation, is commonly reported to be approximately 1000µε (De Souza et al. 2005, Razi et al. 2015, Turner et al. 1994) and was applied to the remodelling algorithm for this study. In contrast, a wide variety of resorption threshold values have been theorised, ranging from 27.9–250µɛ (Cheong et al. 2020a, Frost 1987, Geraldes and Phillips 2014). Considering this range of thresholds, our parameter variation study (Supplementary Fig. 1B) was conducted and an interim value of 100µε was chosen for this model, below which bone resorption would occur.

The remodelling constants for bone tissue formation C1^F and bone tissue resorption C1^R were determined using data obtained from the experimental study from which these murine femurs were CT scanned (Verbruggen et al. 2022), and by performing parameter variation analysis (Fig. 2, Supplementary Fig. 2). Micro-CT analysis from this previous study provides bone tissue density values in the proximal femur cortical and trabecular bone regions of metastatic murine femurs at 3 weeks and 6 weeks post-inoculation (Verbruggen et al. 2022). These exact same regions were reconstructed for the current FE models, and were therefore highly beneficial when predicting the rate at which bone density increases or decreases over time in vivo. Using the remodelling algorithm described previously (Equations 2–7), a single trabecular strut from one metastatic femur FE model was subjected to high strain (1250µε) and C1^F adjusted until density increased from 1.123 g/cm³ to 1.178 g/cm³ to reflect healthy trabecular bone which had increased in mineralisation from 3 weeks and 6 weeks (Verbruggen et al. 2022). A parameter variation study was performed for different C₁^R values when simulating mechanoregulation in the FE models (Fig. 2), which revealed that a C₁^R value of 1.6x10³ predicts bone resorption within the greater trochanter regions, while higher constant values inaccurately predict the majority of the proximal femur model to be resorbed (Fig. 2). Importantly, overt osteolytic degradation in the greater trochanter region was reported previously (Verbruggen et al. 2022), in proximal femurs at 6 weeks post-inoculation. Therefore, C₁^R = 1.6x10³ was applied in the remodelling algorithm for all models in this study as this is more representative of in vivo outcomes.

2.4 Finite Element Analysis

This study sought to predict osteolysis in response to changes in the mechanical environment of metastatic bone, using finite element analysis (FEA) and mechanoregulation theory, and compare these to findings from previous experimental and computational studies. Qualitative analysis of the distribution of bone tissue density and strain energy density (SED) throughout the proximal femur was performed. In terms of quantitative analysis, maximum principal strain distribution was assessed, along with maximum principal strain and SED in these models. Density resorption was also analysed quantitatively, by calculating the percentage volume of each model that had decreased in density throughout the remodelling process. The decrease in density was considered notable if resorbed below 0.1 g/cm³, an interim value between mean model density (1.262 ± 0.02 g/cm³) and complete resorption (0.01 g/cm³), whereby elements with a value below this threshold have reduced in density by over 90%.

Changes in bone tissue density and mechanical strain within the micro-CT derived finite element models of metastatic (MET) proximal femurs at 3 weeks, and the predicted changes arising by 6 weeks after application of the mechanoregulation theory, were analysed. These results were compared to results of our previous FE study (Verbruggen and McNamara 2023), in which micro-CT derived FE models from the 3 week timepoint were identical except for their heterogeneous distribution of density, and also included a separately imaged 6-week cohort of metastatic proximal femur heterogeneous models for direct comparison. The subscript m describes results from homogeneous FE models with applied mechanoregulation theory (i.e. findings from the current study), while subscript s denotes the heterogeneous FE models imaged separately at 3 weeks and at 6 weeks post-inoculation (Verbruggen and McNamara 2023). In this way, MET_s models present the mechanical environment of bone tissue at two time points of the study and can serve as a direct comparison for the 3-week and 6-week time points predicted by the bone-remodelling algorithm, determined according to mechanoregulation theory (MET_m).

2.5 Statistical Analysis

Statistical analyses were performed using MiniTab (version 17) software. Each parameter was assessed for equal variance (F test) and student t-tests were implemented to determine whether averaged data was statistically significant between groups of equal variances. Welsh’s test was applied where sample groups had unequal variance. Results are displayed as mean ± standard deviation, with significance defined as a p value of < 0.05, and further significance also identified (p < 0.01, p < 0.001). Error bars in all bar charts represent standard deviation.

3.1 Predicted density changes in metastatic proximal femurs from 3 to 6 weeks

The bone remodelling algorithm was applied to each model and qualitatively assessed for changes in SED and bone mineral density distributions, beginning prior to the first increment (‘3 weeks’), and at every 7 increments of the algorithm, representing 1 week, to the final time point (‘6 weeks’), see Fig. 3 and Fig. 4. It should be noted that at this 6 week timepoint, the models had not reached homeostasis and remodelling was still active. The models predicted that SED was highest within the cortical bone tissue of the femoral neck and at the loaded surface of the femoral head in all models (Fig. 3). Interestingly, by 6 weeks application of the bone remodelling algorithm predicted elevated bone mineral density in these same cortical neck and femoral head regions (Fig. 4).

3.2 Predicted strain distributions at 3 and 6 weeks of metastasis

The homogeneous µCT-FE models of this study were compared to heterogeneous models previously reported, which had been µCT scanned at both 3- and 6-weeks post-inoculation and analysed under identical boundary conditions at these single time points (3wk METs, 6wk METs). Quantitative analysis of strain distributions within the models (Fig. 5A) show this method has succeeded in demonstrating increased strain distribution (i.e., positive skew) as bone tissue is resorbed. This increase was significant in all ranges up to 250µɛ, with the greatest difference in the 0–50µɛ range (3wk MET_m: 88.33 ± 3.96%, 6wk MET_m: 23.86 ± 3.96%, p < 0.001). Mean values of maximum principal strain and SED had increased significantly in both the MET_s and MET_m models between the two time points in question (Fig. 5B, C). A parameter variation analysis showed notable differences in strain distribution between models when the resorption threshold was altered (see Supplementary Fig. 1B). However, even the largest model volume variation (11.9%, 50–100µɛ range) would not have impacted the significant differences found between day 1 (9.89 ± 2.8%) and day 21 (34.88 ± 5.65%) within this strain range.

By the 6-week time point, application of the remodelling algorithm predicted that 26.51 ± 9.96% of the model volume decreased below a bone mineral density of 0.1 g/cm³. In comparison, our previous experimental study reported that bone volume in the trabecular regions of proximal femurs reported from prior micro-CT analysis had decreased by 21.15% between 3 weeks and 6 weeks post-inoculation of metastatic breast cancer cells (Verbruggen et al. 2022). In addition, the area of lowest density in each model consistently presented in the greater and lesser trochanter regions, as is illustrated in the anterior-posterior cross-sectional views of cortical and trabecular bone (Fig. 6A). This correlates with findings from our experimental study in which osteolysis was reported in the greater trochanter regions by 6 weeks post-inoculation (Verbruggen et al. 2022) (Fig. 6B, 7C). Of note, density was highest in the cortical bone tissues of the femoral neck and the inter-trochanter space (Fig. 6A).

This study is the first to implement a computational mechanoregulation framework to study and predict the development of osteolysis during breast cancer-bone metastasis. Subject-specific finite element (FE) models were developed and applied to predict an altered mechanical environment within bone tissue during early stage metastasis in an animal model (3 weeks post-inoculation). Our bone mechanoregulation algorithm was applied to represent 3–6 weeks of osteolysis development and successfully predicted osteolysis as a function of the altered mechanical environment. Specifically, we predicted a decrease in bone mineral density in the greater trochanter regions of the femur, the exact regions where osteolytic lesions were identified in an associated experimental study (Verbruggen et al. 2022). Moreover, application of the mechanoregulation algorithm predicted that the mechanical environment evolved in a similar manner to that predicted through subject-specific finite element (FE) models (Verbruggen and McNamara 2023). These results support the hypothesis that early changes in the physical environment of bone tissue during metastasis may elicit mechanobiological cues for bone cells and activate osteolytic destruction.

Some limitations in this study must be considered. Firstly, bone tissue was assumed to be linear elastic and homogeneous, which does not fully account for the heterogeneous and non-linear behaviour of bone tissue. However, linear elasticity has been assumed in previous µCT-FE models implementing mechanoregulation theory (Cheong et al. 2020a, Cheong et al. 2020b, Schulte et al. 2011, Schulte et al. 2013), necessitated by the need for computational efficiency. It was not possible to implement heterogeneous material properties at the outset in combination with the mechanoregulation framework presented here, because this resulted in large discrepancies between material properties of neighbouring elements and associated convergence challenges. Interestingly, a previous study predicted that homogenous models significantly overestimate apparent elastic moduli when compared to heterogeneous µCT-FE models (Renders et al. 2011), whereas a recent study reported that heterogeneous models did not predict failure loads and stiffness as closely to experimental results, when compared to homogeneous models (Oliviero et al. 2021). The approach predicted a lower strain distribution for the homogenous model at 3 weeks compared to our previous micro-CT derived FE models, which accounted for heterogeneous material properties (Verbruggen and McNamara 2023), but the material properties evolved to become heterogeneous over the simulation period according to the mechanoregulation theory, with the mechanical environment evolving in a similar manner to that predicted through subject-specific FE model. Further validation could be conducted by longitudinal Digital Volume Correlation (DVC) analysis (Palanca et al. 2022). Secondly, we assumed a constant remodelling rate (C₁) for cortical and trabecular bone, because manual delineation of cortical and trabecular bone (to implement separate C₁ constants) could introduce complex surface interactions not reflective of in vivo conditions, and, because there was a low volume of tissue above the bone formation threshold (< 0.001%), no sensitivity analysis was performed on the formation rate, C₁^F. Moreover, we did not account for biochemical signalling, which regulates interactions between tumour and bone cells. Previous studies have coupled biochemical signalling between osteoblasts and osteoclasts with mechanoregulation theory to study bone remodeling (Hambli 2014) or accounted for varying bone tissue resorption and formation rates according to experimental in vivo measurements (Cheong et al. 2020a). While we successfully demonstrated osteolysis over a 3-week period, the models did not reach homeostatic equilibrium in this timeframe, and therefore do not account for longer term mechanoregulatory responses. Variations in loading angle and magnitude during locomotion may arise (Charles et al. 2018) were not investigated. Future studies could investigate the influence of heterogeneity of material properties, biochemical signalling and associated remodelling activity, and loading on predicted bone tissue resorption. Nonetheless, the homogenous model presented here, in which bone remodelling is governed by mechanical stimuli alone, is in agreement with our experimental data and finite element models for 6-week-old metastatic animals (Verbruggen et al. 2022, Verbruggen and McNamara 2023).

A parameter variation study was conducted to establish a resorption (C₁^R = 1.6x10³), determined by parameter variation analysis (Supplementary Fig. 1), which could predict the degree of bone loss observed in experimental findings. Interestingly, application of the mechanoregulation theory predicted that principal strain and SED increased significantly from 3 to 6 weeks, which corresponds well with the predictions at the respective time-points from micro-CT derived heterogeneous finite element models (Verbruggen and McNamara 2023). Specifically, the identified resorption constant predicted that approximately 27% of the model volume would reduce to low bone density, which was within the same range as bone volume fraction reductions in the trabecular space (21.15%) reported experimentally (Verbruggen et al. 2022). At physiological loads, the predicted bone tissue resorption was more prominent within these models than bone tissue formation. This finding correlates well with similar studies of female C57BL/6 murine tibiae, wherein bone tissue formation only occurred at high loads (9–13N) (De Souza et al. 2005, Weatherholt et al. 2013). It is notable that implementing a lower resorption rate (C₁^R = 160), predicted an increase in strain distribution, maximum principal strain and SED (Supplementary Fig. 2). Further study on the resorption rates could shed light on the evolving mechanical environment and sensitivity to these remodelling parameters.

Physiological loads were applied to subject-specific proximal femur models and a bone mechanoregulation algorithm was implemented to predict element-specific density changes over 21 days. Qualitative analysis revealed that the lowest bone mineral density values were predicted to arise in the greater and lesser trochanter of the femurs. Interestingly, a second cohort of proximal femurs from the same animal model of bone metastasis were experimentally analysed (by micro-CT) 6 weeks after tumour-inoculation and these animals consistently presented with overt osteolytic lesions in the greater trochanter region (Verbruggen et al. 2022). Thus, the application of the mechanoregulation theory has accurately predicted the spatial nature of bone resorption within the proximal femur that corresponds to previous experimental results. Interestingly, a study that computationally investigated bone remodelling in a healthy human proximal femur over 5 years of physiological loads similarly predicted that density was lowest in the femoral trochanter (Hambli 2014). In the current study bone mineral density was predicted to be highest along the femoral neck region and femoral head surface, which corresponds to previously reported CT-imaged and DXA-imaged density maps of a healthy human proximal femur under stance and side-loading configurations (Dall’Ara et al. 2016). Moreover, strain distributions had increased in the femoral neck cortical bone in a human osteoporotic femur that had undergone bone loss and trabecular thinning (Verhulp et al. 2008). Interestingly, significant bone tissue resorption was reported experimentally by micro-CT analysis of the femoral head regions after 6 weeks of bone metastasis development, and indeed the femoral head was absent in 3 of 7 metastatic proximal femurs (Verbruggen et al. 2022). These may have been caused due to fracture along the femoral neck, and the results presented here predict that maximum principal strains are highest and continue to increase as metastasis progresses. Overall, application of the computational bone mechanoregulation framework presented here demonstrated the spatial nature of resorption by 6 weeks that correlated to the experimental findings.

In conclusion, this study is the first of its kind to predict, using mechanoregulation theory, bone remodelling in µCT-derived finite element models of metastatic bone tissue, prior to the development of overt osteolytic lesions. The bone remodelling algorithm predicted bone mineral density to decrease in regions that coincide with experimental studies of osteolysis. This study also reported changes in strain distribution from 3 weeks to 6 weeks, which were in keeping with predictions of micro-CT derived FE models of femurs at 3 weeks and 6 weeks post-inoculation of metastatic cells. Thus we propose that mechanobiology may play a role in the adaption of the bone tissue extracellular matrix to metastasis and contribute to the later development of osteolysis.

Acknowledgements: This publication has emanated from research conducted with financial support of the Irish Research Council Laureate Award Programme 2017/18 (MEMETic, IRCLA/2017/217), the Hardiman Postgraduate Research Scholarship, provided by the University of Galway (UG) in Ireland, and Science Foundation Ireland under Grant number [18/SPP/3522] as part of Precision Oncology Ireland, with co-funding from the National Breast Cancer Research Institute.

Author contributions: A.S.K.V.: Writing – review & editing, Writing – original draft, Visualisation, Methodology, Investigation, Formal analysis, Conceptualization. E.C.M.: Model design. R.M.D.: Resources, Funding acquisition, Model design. L.M.M.: Writing – review & editing, Writing – original draft, Supervision, Project administration, Investigation, Funding acquisition, Conceptualization.

Competing interests: The authors have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Araujo A, Cook LM, Lynch CC, Basanta D (2014) An integrated computational model of the bone microenvironment in bone-metastatic prostate cancer. Cancer Res 74(9):2391–2401
Arrington SA, Schoonmaker JE, Damron TA, Mann KA, Allen MJ (2006) Temporal changes in bone mass and mechanical properties in a murine model of tumor osteolysis. Bone 38(3):359–367
Charles JP, Cappellari O, Hutchinson JR (2018) A dynamic simulation of musculoskeletal function in the mouse hindlimb during trotting locomotion. Front Bioeng Biotechnol 6:61
Cheong VS, Marin C, Lacroix A, D. and, Dall’Ara E (2020a) A novel algorithm to predict bone changes in the mouse tibia properties under physiological conditions. Biomech Model Mechanobiol 19(3):985–1001
Cheong VS, Roberts BC, Kadirkamanathan V, Dall'Ara E (2020b) Bone remodelling in the mouse tibia is spatio-temporally modulated by oestrogen deficiency and external mechanical loading: A combined in vivo/in silico study. Acta Biomater 116:302–317
Cook LM, Araujo A, Pow-Sang JM, Budzevich MM, Basanta D, Lynch CC (2016) Predictive computational modeling to define effective treatment strategies for bone metastatic prostate cancer. Sci Rep, 6, pp. 29384
Dall’Ara E, Eastell R, Viceconti M, Pahr D, Yang L (2016) Experimental validation of DXA-based finite element models for prediction of femoral strength. J Mech Behav Biomed Mater 63:17–25
De Souza R, Matsuura M, Eckstein F, Rawlinson S, Lanyon L, Pitsillides A (2005) Non-invasive axial loading of mouse tibiae increases cortical bone formation and modifies trabecular organization: a new model to study cortical and cancellous compartments in a single loaded element. Bone 37(6):810–818
Frost HM (1987) Bone “mass” and the “mechanostat”: a proposal. Anat Rec 219(1):1–9
Geraldes DM, Phillips AT (2014) A comparative study of orthotropic and isotropic bone adaptation in the femur. Int J Numer Methods Biomed Eng 30(9):873–889
Hambli R (2014) Connecting mechanics and bone cell activities in the bone remodeling process: an integrated finite element modeling. Frontiers in Bioengineering and Biotechnology, 2, pp. 6
Huiskes R, Weinans H, Grootenboer H, Dalstra M, Fudala B, Slooff T (1987) Adaptive bone-remodeling theory applied to prosthetic-design analysis. J Biomech 20(11–12):1135–1150
Kim Y, Othmer HG (2015) Hybrid models of cell and tissue dynamics in tumor growth. Math Biosci Eng 12(6):1141–1156
Lu Y, Zuo D, Li J, He Y (2019) Stochastic analysis of a heterogeneous micro-finite element model of a mouse tibia. J Med Eng Phys 63:50–56
McNamara LM, Prendergast PJ (2007) Bone remodelling algorithms incorporating both strain and microdamage stimuli. J Biomech 40(6):1381–1391
Mulvihill BM, McNamara LM, Prendergast PJ (2008) Loss of trabeculae by mechano-biological means may explain rapid bone loss in osteoporosis. J Royal Soc Interface 5(27):1243–1253
Nazarian A, von Stechow D, Zurakowski D, Muller R, Snyder BD (2008) Bone volume fraction explains the variation in strength and stiffness of cancellous bone affected by metastatic cancer and osteoporosis. Calcif Tissue Int 83(6):368–379
Oliviero S, Roberts M, Owen R, Reilly G, Bellantuono I, Dall’Ara E (2021) Non-invasive prediction of the mouse tibia mechanical properties from microCT images: comparison between different finite element models. Biomech Model Mechanobiol 20(3):941–955
Paget S (1889) The distribution of secondary growths in cancer of the breast. Lancet, pp. 571–573
Palanca M, Oliviero S, Dall’Ara E (2022) MicroFE models of porcine vertebrae with induced bone focal lesions: Validation of predicted displacements with digital volume correlation. J Mech Behav Biomed Mater 125:104872
Pereira AF, Javaheri B, Pitsillides A, Shefelbine S (2015) Predicting cortical bone adaptation to axial loading in the mouse tibia. J Royal Soc Interface 12(110):20150590
Quinn C, Kopp A, Vaughan TJ (2022) A coupled computational framework for bone fracture healing and long-term remodelling: Investigating the role of internal fixation on bone fractures. Int J Numer Methods Biomed Eng 38(7):e3609
Ramos A, Simoes J (2006) Tetrahedral versus hexahedral finite elements in numerical modelling of the proximal femur. Med Eng Phys 28(9):916–924
Razi H, Birkhold AI, Weinkamer R, Duda GN, Willie BM, Checa S (2015) Aging leads to a dysregulation in mechanically driven bone formation and resorption. J Bone Miner Res 30(10):1864–1873
Rejniak KA, Anderson AR (2011) Hybrid models of tumor growth. Wiley Interdiscip Rev Syst Biol Med 3(1):115–125
Renders G, Mulder L, Van Ruijven L, Langenbach G, Van Eijden T (2011) Mineral heterogeneity affects predictions of intratrabecular stress and strain. J Biomech 44(3):402–407
Richert L, Keller L, Wagner Q, Bornert F, Gros C, Bahi S, Clauss F, Bacon W, Clézardin P, Benkirane-Jessel N, Fioretti F (2015) Nanoscale Stiffness Distribution in Bone Metastasis. World J Nano Sci Eng 05(04):219–228
Scannell PT, Prendergast PJ (2009) Cortical and interfacial bone changes around a non-cemented hip implant: simulations using a combined strain/damage remodelling algorithm. Med Eng Phys 31(4):477–488
Schulte FA, Lambers FM, Webster DJ, Kuhn G, Müller R (2011) In vivo validation of a computational bone adaptation model using open-loop control and time-lapsed micro-computed tomography. Bone 49(6):1166–1172
Schulte FA, Zwahlen A, Lambers FM, Kuhn G, Ruffoni D, Betts D, Webster DJ, Müller R (2013) Strain-adaptive in silico modeling of bone adaptation—a computer simulation validated by in vivo micro-computed tomography data. Bone 52(1):485–492
Sohail A, Younas M, Bhatti Y, Li Z, Tunç S, Abid M (2019) Analysis of trabecular bone mechanics using machine learning. Evolutionary Bioinf 15:1176934318825084
Tracqui P (2009) Biophysical models of tumour growth. Rep Prog Phys, 72(5)
Turner CH, Forwood M, Rho JY, Yoshikawa T (1994) Mechanical loading thresholds for lamellar and woven bone formation. J Bone Miner Res 9(1):87–97
Van Rietbergen B, Huiskes R, Weinans H, Sumner D, Turner T, Galante J (1993) The mechanism of bone remodeling and resorption around press-fitted THA stems. J Biomech 26(4–5):369–382
Verbruggen AS, McCarthy EC, Dwyer RM, McNamara LM (2022) Temporal and spatial changes in bone mineral content and mechanical properties during breast-cancer bone metastases. Bone reports, pp. 101597
Verbruggen AS, McNamara LM (2023) Mechanoregulation may drive osteolysis during bone metastasis: A finite element analysis of the mechanical environment within bone tissue during bone metastasis and osteolytic resorption. J Mech Behav Biomed Mater, pp. 105662
Verhulp E, van Rietbergen B, Huiskes R (2008) Load distribution in the healthy and osteoporotic human proximal femur during a fall to the side. Bone 42(1):30–35
Villette C, Zhang J, Phillips A (2020) Influence of femoral external shape on internal architecture and fracture risk. Biomech Model Mechanobiol 19(4):1251–1261
Weatherholt AM, Fuchs RK, Warden SJ (2013) Cortical and trabecular bone adaptation to incremental load magnitudes using the mouse tibial axial compression loading model. Bone 52(1):372–379
Zhou X, Liu J (2014) A computational model to predict bone metastasis in breast cancer by integrating the dysregulated pathways. BMC Cancer 14:618

No competing interests reported.

floatimage7.png
Supplementary Fig 1: Parameter sensitivity analyses. Distribution of maximum principal strain as a percentage of bone volume in one metastatic proximal femur model, predicted at 6 weeks, comparing between values of (A) bone tissue resorption constant, C₁^R_, and (B) resorption threshold strain value (ɛ_res) applied to one proximal femur model. The greatest difference in percentage volume for altered ɛ_res was between 75 and 200µɛ in the 50-100 range, at 11.09%. Generated using GraphPad
floatimage8.png
Supplementary Fig 2: Predicted changes in bone tissue strain at 3 and 6-week time points compared to bone remodelling algorithm predictions, where C₁^R is 160. Heterogeneous FE proximal femur models from a previous study at 3 and 6 weeks post-inoculation, loaded and analysed at a single time point, (3wk MET_s, 6wk MET_s) were compared to bone remodelling algorithm simulations of these same models with homogeneous material properties, at the 1^st and 21^st increments (3wk MET_m, 6wk MET_m). (A) Distribution of maximum principal strain as a percentage of bone volume, including inset of higher microstrain ranges (150-200µɛ to 300µɛ+). (B) Maximum principal microstrain and (C) SED (Pa) in both µCT-captured and mechanoregulation-driven FE models Generated using GraphPad

Download PDF

Version 1

posted

You are reading this latest preprint version

Predicting mechanoregulatory responses in bone during breast cancer metastasis: A Finite Element Analysis

Status:

Version 1

Abstract

Figures

1.0 Introduction

2.0 Methods

2.1 Obtaining bone tissue geometry and material properties

2.2 Mechanoregulation Theory

2.3 Determining bone remodelling thresholds and constants

2.4 Finite Element Analysis

2.5 Statistical Analysis

3.0 Results

3.1 Predicted density changes in metastatic proximal femurs from 3 to 6 weeks

3.2 Predicted strain distributions at 3 and 6 weeks of metastasis

4.0 Discussion

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1