CT-based finite element simulating spatial bone damage accumulation predicts metastatic human vertebrae strength and stiffness.

Текст сохранён. Часть иллюстраций пока доступна только в первоисточнике; загрузка будет проверена повторно.

1 Introduction

In 2022, more than 1.9 million new cancer cases were estimated in the US ( Rebecca et al., 2022 ). Spinal bone metastases affect more than 70% of patients with advanced cancer ( Cronin et al., 2018 ). Nearly half of the spinal metastasis patients undergo vertebral radiation therapy (RT) ( Begg et al., 2011 ), with up to 40% of these patients suffering clinically significant vertebral fractures (VF) post-RT ( van der Velden et al., 2017 ; Faruqi et al., 2018 ). Once VF occurs, patients may suffer severe impairment of quality of life ( Pond et al., 2014 ), higher health costs ( McDougall et al., 2016 ), and neurological deficits in up to 50% of patients with VF ( Oster et al., 2013 ), shortening patient survival ( Saad et al., 2007 ; Oster et al., 2013 ) and 3-year life expectancy ( Oefelein et al., 2002 ; Pond et al., 2014 ). Clinical guidelines for predicting PV risk remain subjective with low accuracy VF ( Yao et al., 2017 ). Establishing the risk of pathologic vertebral fractures before catastrophic pain or neurologic deficits occur remains an unmet, critical clinical need in managing patients with spinal metastatic disease.

The invasion of vertebral bone with bone metastases destroys vertebral anatomy ( Taneichi et al., 1997 ). It disrupts the bone’s cellular homeostasis ( Wu et al., 2018 ), causing remarkable degradation of the vertebral bone composition ( Burke et al., 2016 ; Burke et al. 2017 ; Burke et al. 2017 ) and microarchitecture ( Burke et al., 2016 ; Burke et al., 2017 ; Bailey et al., 2020 ; Bailey et al., 2022 ). Using a pre-clinical model of osteolytic and mixed vertebral bone tissue, ( Atkins et al., 2019 ), showed this diminished bone quality to be associated with a higher accumulation of diffuse and linear microdamage and microfractures at regions of high stress when exposed to mechanical loading. These changes underly the well-documented worsening of the pathologic bone tissue’s mechanical properties ( Nazarian et al., 2008 ; Burke et al., 2018 ; Bailey et al., 2020 ; Stadelmann et al., 2020 ).

Derived from CT data, computational finite element (FE) models allowed detailed insights into the effect of a range of osteolytic defects, simulated as idealized void-based geometries, on the failure strength of cadaver vertebra ( Whyne et al., 2001 ; Tschirhart et al., 2004 ; Tschirhart et al. 2006 ; Tschirhart et al. 2007 ; Alkalay and Harrigan, 2016 ). In multiple myeloma patients with VF ( Campbell et al., 2017 ) and in cancer patients with osteolytic vertebral lesions ( Costa et al., 2019 ), subject-specific finite element (FE) highlighted the bone metastases’ deleterious effect on the patient’s vertebral mechanical competence, defined as lower vertebral stiffness and work to yield, as compared to cancer patients without observed vertebral fracture. Recent work employing digital image correlation ( Palanca et al., 2023 ) demonstrated osteolytic lesions to cause significant alteration in the spatial strain distribution within the vertebrae and endplate, suggesting this effect as the mechanism for the observed vertebral failure. In vertebrae containing osteosclerotic lesions, Palanca et al. (2023) reported strain peaks within the vertebral volume to be largely focused around the osteosclerotic lesion and were associated with vertebral failure. However, the effect of osteosclerotic lesions on the mechanism of VF remains poorly understood.

A recent review of the CT-based finite element modeling for pathologic vertebrae highlighted the limitations of idealized and simplified lesion geometries and material models employed when developing patient-specific prediction of vertebral strength ( Molinari and Falcinelli, 2022 ). Such development is crucial for developing clinical tools tailored to individualized treatments and diagnosis. Incorporating material characterization of osteolytic and osteosclerotic human bone, Stadelmann et al. (2020) compared a homogenized FE (hFE) model with that of a micro FE model for predicting the strength and stiffness of cadaveric vertebrae containing a wide range of osteolytic, osteosclerotic, and mixed bone metastases obtained from cancer patients. Although the study demonstrated the hFE approach utility for predicting the experimental strength, the authors did observe bone metastasis type, i.e., osteolytic vs. osteosclerotic vs. mixed (having both osteolytic vs. osteosclerotic within the vertebral volume), to markedly affect the FE models prediction accuracy for strength and stiffness. In total, these findings suggest that osteolytic and osteosclerotic lesions may differentially affect the process of damage evolution within the vertebral bone ( Garcia et al., 2009 ) that affects vertebral failure ( Jackman et al., 2016 ).

Bone is a multiscale, hierarchal biomaterial with complex mechanical responses at each scale ( Schwiedrzik et al., 2013 ). A common approach for studying the effect of damage evolution on bone strength is to use FE simulations based on so-called continuum damage models, in which a single scalar variable D evolves irreversibly from 0 (no damage) to at most 1 (completely damaged, i.e., the material’s stiffness is fully diminished) according to a damage evolution law. Zysset and Curnier (1996) developed a general material model integrating isotropic damage and time-independent plastic flow, later expanded by Dall’Ara et al. (2013) with the central assumption that plastic flow and damage accumulation are intrinsically related. Incorporated into a three-dimensional continuum model of anisotropic trabecular bone, the model demonstrated strong correlations between predictions and experimental results for none-pathologic bone ( Garcia et al., 2009 ; Wolfram et al., 2011 ) and vertebral strength ( Chevalier et al., 2008 ; Dall’Ara et al., 2010 ). To describe a multi-axial yield and failure criterion for trabecular bone based on BV/TV, a scalar function of stress demarcates the boundary between the intact and damaged states of the bone. This function will consider the trabecular bone’s anisotropic and heterogeneous nature in the form of a halfspace generalization of the Hill criterion.

Alternative damage-based approaches employing elastic-plastic material failure ( Molinari and Falcinelli, 2022 ) and phase field model (PFM) ( Preve et al., 2024 ) have been proposed for investigating the failure of transpedicular screws in metastatic vertebrae. The former relies on a maximum stress criterion with the failed element assigned a negligible modulus to simulate failure. The PFM model relies on a critical energy release rate obtained using a power-law equality based on the bone density properties, assumed as locally soft for osteolytic bone and locally stiff for osteoblastic lesions with both lesion’s geometries modeled as spherical regions. In brittle, porous materials under compressive, the PFM model demonstrated the failure of load to be influenced by the size of the voids, with small voids promoting damage nucleation and enhancing the bridging of macro-pores through micro-crack formation ( Cavuoto et al., 2024 ). In contrast, macro-pores affect the overall material response and drive the propagation of large fractures. These findings present a possible mechanism for the failure of osteolytic bone that does not demonstrate critical lytic foci but rather a “moth-eaten” presentation common to many such vertebrae ( Bailey et al., 2022 ).

The study aims were to 1) extend the constitutive bone model proposed by Johnson et al. (2010) , implemented in the ΣMIT FE framework ( Radovitzky et al., 2011 ) to include a continuum bone damage computational framework. 2) Evaluate the new model’s accuracy for predicting the strength and stiffness of human vertebrae containing osteolytic osteosclerotic and mixed bone metastases 3) evaluate whether a single set of calibrated bone material parameters can be established to accurately predicts vertebral strength and stiffness independent of bone metastatic type. For this study, we simulated the measured strength and stiffness of ten cadaveric pathologic vertebrae containing osteolytic, osteosclerotic, and mixed bone metastases that were mechanically tested in a previous study ( Stadelmann et al., 2020 ). We employed a machine learning approach to create a global material calibration scheme that would apply to all the bodies and tested the performance of specimen-specific vs. “global” material calibration on the model’s vertebral strength and stiffness prediction. We hypothesized that 1) the damage-based FE framework simulates the observed values for strength and stiffness and 2) a single set of model material constants provides a close simulation of the measured strength and stiffness across all bone metastasis types.

2 Materials and methods

2.1 Experimental measurements

For this study, we selected to simulate ten cadaveric vertebrae containing osteolytic, osteosclerotic or mixed bone metastases from a sample of 45 pathologic vertebrae that we mechanically tested in a previous study ( Stadelmann et al., 2020 ). The vertebrae selected present a wide range of bone volume fraction values with either a high or a low measured strength or stiffness or a high difference in these mechanical properties. Hence, they represent the most challenging cases for modeling the measured experimental mechanical behavior.

2.1.1 Specimens

Ten vertebral samples were obtained as part of our previous study ( Stadelmann et al., 2020 ) from nine cadaver spines (three female, six male, age 49–71 years, mean age 54) from donors with solid cancer (three breast, three lung, two prostate, and one kidney), Table 1 . As part of the previous study ( Stadelmann et al., 2020 ), each spine was imaged in a clinical CT scanner, and upon radiographic review, vertebral segments exhibiting osteolytic, osteosclerotic, or mixed bone metastases were extracted. Following the protocol proposed by Dall’Ara et al. (2010) , the vertebral end plates were sectioned to obtain plane-parallel segments, and the sectioned vertebrae prepared for micro CT scanning ( Stadelmann et al., 2020 ).

2.1.2 MicroCT imaging

Stadelmann et al. (2020) : in brief, each vertebral section was imaged at 24.5 μm isotropic voxel size (μCT100, Scanco Medical, Switzerland) using the following parameters (tube voltage: 70 kV, tube current: 200 μA, integration time: 500 ms). A Gauss filter (Sigma: 0.8, Support: 1) ( Chevalier et al., 2009 ) was applied to reduce high-frequency noise within the images, and the vertebra was segmented using a custom script (IPL, Scanco Medical, Switzerland). We computed the vertebral bone mineral density (BMC) for each image within the segmentation volume and, per image, applied an adaptive threshold algorithm to compute an optimal bone segmentation threshold ( Ding et al., 1999 ; Burghardt et al., 2007 ). The mean threshold value (429 ± 56 mgHA/cm 3 ) was applied to segment the bone tissue within the vertebral volume, and the resulting μCT data was used to derive the overall (i.e., cortical + trabecular regions) bone volume fraction (BV/TV).

2.1.3 Mechanical experiment

Stadelmann et al. (2020) : in brief, our study followed the standardized mechanical test protocol described by Dall’Ara et al. (2010) . First, the specimen’s center of mass was computed from the μ CT images, digitally shifted anteriorly by a distance equal to 10% of the anteroposterior width of the vertebral bottom surface, and the image-derived location of force application and outer contour of the bottom surface printed on a sheet of paper to allow precise location of the specimen in the test system. This methodology aims to induce an anterior wedge-shaped fracture, a common fracture pattern in osteoporotic and cancer patients ( Dall’Ara et al., 2010 ).

Each vertebral level was equilibrated in 0.9% NaCl saline solution for 1 h at room temperature, then placed and positioned based on the printed information on a fixed steel platen secured to the hydraulic testing system (858 Mini Bionix II, MTS, Eden Prairie, United States). With the specimen positioned, the cranial plate, having a ball joint mechanism to allow unconstrained deformation of the sample, was lowered until a tare load of 25 N was recorded to confirm contact. Both steel plates had their contact surfaces sandblasted to prevent the sample from sliding during the compression test. The vertebra was tested under monotonic uniaxial compressive displacement at a rate of 5 mm/min ( Chevalier et al., 2009 ) until either a failure was registered or the built-in load cell (model: 662.20D-04) maximum force (15 kN) was reached. We defined vertebral strength (Fexp) as the maximum measured compressive force. Experimental stiffness was derived as the coefficient of the linear regression model fitted to 20%–80% of the load-displacement curve elastic region prior to the vertebral yield strength.

2.2 Establishing the damage-based FE framework

2.2.1 Establishing elastic-plastic-damage model for vertebral bone

The model to simulate the pathologic vertebral bone elastic and plastic response was developed based on the constitutive law proposed by Johnson et al. (2010) . Specifically, the modification incorporated a) bone volume fraction, B V T V , to capture the heterogeneity of bone volume density throughout the vertebrae and b) a three-dimensional elastic-plastic damage constitutive law for bone tissue ( Garcia et al., 2009 ). The model consists of viscoelastic and viscoplastic components in series ( Figure 5 ), in which the deformation gradient is decomposed into the viscoelastic, F V E , and viscoplastic, F V P , components as Eq. 1 : F = F V E F V P (1)

The viscoelastic component was modeled using a Maxwell–Wiechert model with viscosity parameters of η 1 and η 2 to capture two separate low and high-strain rate viscoelastic mechanisms within the microstructure of the bone. The rate sensitivity at lower and higher strain rates is likely due to the collagen fibers’ nature and the presence of interstitial fluid within the microstructure, respectively ( Johnson et al., 2010 ). The total Cauchy stress, T , of the model, can be calculated by simply summing the individual Cauchy stress as, Eq. 2 : T = T E 0 + T E 1 + T E 2 (2) where the Cauchy stress of the fully elastic branch, T E 0 , is calculated as, Eq. 3 : T E 0 = F V E S E 0 F V E T det F V E (3) and the Cauchy stress of the viscoelastic branches, T E 1 and T E 2 , can be obtained by employing the central difference method to the deviatoric viscoelastic strain rate, ε ˙ V E ′ , expressed as, Eq. 4 ε ˙ V E ′ = 3 2 η i B V T V k T V E ′ i + 3 2 μ 0 i B V T V k T ˙ V E i For i = 1 , 2 (4) where the bone is assumed to be an isotropic material with the shear, μ 0 i , and bulk, k 0 i , moduli, Eqs 5 , 6 : μ 0 i = E 0 i 2 1 + ν (5) k 0 i = E 0 i 3 1 − 2 ν (6)

In this way, the second Piola–Kirchoff stress of the fully elastic branch, S E O , in Eq. 3 can be calculated as, Eq. 7 S E O = 2 μ 0 0 B V T V k ε V E ′ + K 0 0 t r ε V E I (7) where I is the identity matrix. Also, the deviatoric viscoelastic strain, ε V E ′ , can be obtained as, Eq. 8 ε V E ′ = ε V E − 1 3 t r ε V E I (8)

Similarly, after calculating the total Cauchy stress from Eq. 2 , the deviatoric Cauchy stress, T ′ , can be calculated as, Eq. 9 : T ′ = T − 1 3 t r T I (9)

From there, the magnitude, τ , and the direction, N , of the deviatoric Cauchy stress can be obtained as, Eqs 10 , 11 : τ = 1 2 t r T ′ T T ′ (10) N = 1 τ 2 T ′ (11)

The model viscoplastic behavior is calculated based on the co-directionality assumption for the plastic flow and the deviatoric total Cauchy stress, which results in the following rate of plastic stretching, Eq. 12 D V P = τ 3 S 0 m N (12) with m as the inverse of the slope of the log-log relationship between the strain rate and plastic yield stress and the natural logarithm of S 0 as the y -intercept of this relationship. From there, the rate of change of the plastic deformation gradient is calculated as follows, Eq. 13 : F ˙ V P = F V E − 1 D V P F (13)

2.2.2 Damage model

The damage model is derived based on the work by Dall’Ara et al. (2013) based on the elastic-plastic-damage model proposed by Garcia et al. (2009) , in which a Halfspacewise Hill criterion is defined as, Eq. 14 : Y S , D = S : F ± S − r D (14)

The damage hardening law, which adjusts the yield surface, is defined as, Eq. 15 : r D ≔ 1 − 1 − α e − k D (15)

With α: the ratio between the yield and ultimate stress. The fourth-order tensors, F ± , defines the different yield properties for compression and tension and is expressed as, Eq. 16 : F ± = 1 T F 2 ∑ i = 1 3 1 σ 0 ± 2 M i ⊗ M i − ∑ i , j = 1 ; i # j 3 χ 0 ± σ 0 ± 2 M i ⊗ M j + ∑ i , j = 1 ; i # j 3 χ 0 ± 2 τ 0 2 M i ⊗ ¯ _ M j (16) where σ 0 ± is the uniaxial tensile and compressive strengths, τ 0 is the shear strength, χ 0 ± is the stress interaction coefficients, and TF is the tissue function used to scale the model for dense bone (BV/TV > 0.5), Eq. 17 : TF BV / TV , E max , E 0 , k = BV / TV k if BV / TV ≤ 0.5 BV / TV k + E max − E 0 E 0 B V / T V − 0.5 0.5 2 k if BV / TV > 0.5 (17)

M is the second-order fabric anisotropy tensors defined as: M i ≔ m i m i ⊗ m i (18)

where i = 1,2,3. The eigenvectors m i correspond to the normal direction of the orthotropic symmetric planes, and the eigenvalues m i express the extent of anisotropy since we model the bone as isotropic, m i = 1.

2.3 Homogenized finite element (hFE) model

Based on Pahr and Reisinger (2020) , the 24.5 μm μCT data was resampled to a resolution of 0.3185 mm isotropic per voxel (Medtool ver. 1.4™, Dr. Pahr Ingenieurs e.U). We selected this image resolution to represent the in-plane clinical CT data resolution (0.3185 mm) obtained under our ongoing patient study with a slice thickness used in clinical CT images (0.625 mm). For each image voxel, bone mineral density (BMD) was assigned based on an empirical conversion law relating CT HU to a BMD value followed by a threshold operation (390 mg/cm 3 ). We used Medtool 1.4™, a computational geometry and 3D mesh generation software library ( https://www.cgal.org ), to mesh the masked volumes with a tetrahedral mesh (element size = 1.0 mm). Mesh quality was improved using the HealMesh software library ( Mauch et al., 2006 ) geometrical and topological mesh optimization algorithms. Subsequently, each mesh element was assigned a local BV/TV by interpolating the segmented and masked image (Medtool 1.4™).

2.3.1 Material model parameter calibration

The boundary value problem was initially solved for a selected hFE model with a set of initial reference material model parameters. An iterative procedure was employed where the most relevant material parameters were incrementally adjusted from their reference values in the direction in which they would reduce the RMSE of the simulated load-displacement curve relative to the experimentally measured ( Burke et al., 2018 ). The process was stopped when RMSE was equal to or lower than 5%. The initial trial values of the viscoelastic, viscoplastic cortical bone and damage model parameters were taken directly from Siegel et al. (2012) . To reduce the degrees of freedom in the calibration process, we chose to independently calibrate only material model parameters deemed relevant for applying quasi-static, monotonic compressive loading. The following set of material model parameters were adjusted proportionally. 1) The Young’s moduli for the viscoelastic branches, E 0 1 and E 0 2 , were adjusted by the same factor as that of the equilibrium response, E 0 0 . The viscosity parameters for the viscoelastic branches, η 0 1, and η 0 2 , were also adjusted by the same factor as E 0 0 with an additional factor of 1/100 to reduce rate-dependency effects further for our quasi-static study. Similarly, the uniaxial tensile strength σ 0 + and shear strength τ 0 were reduced by the same factor as the uniaxial compressive strength σ 0 − . Non-dimensional parameters such as the parameter governing the relationship between strain rate and yield stress m , the stress interaction coefficients χ 0 ∓ , the two parameters defining the viscous damage behavior χ 0 ∓ and ζ, and the ratio between the yield and ultimate stress α, were kept constant. The resulting material model parameters for each vertebra are reported in Table 2 .

2.3.2 Finite element (FE) model simulations

A C++ driver to model the quasi-static loading of vertebrae under compression was implemented within ∑MIT ( http://summit.mit.edu ), a computational solid mechanics framework developed by Professor Raúl Radovitzky’s group ( Radovitzky et al., 2011 ), providing parallel computing large-scale simulations. For the simulation, we imposed a monotonically increasing uniaxial compressive displacement along the craniocaudal direction on the upper surface of the vertebral mesh model and an encastre boundary condition on the inferior surface of the vertebral mesh model. At each iteration step, the imposed displacement and compressive force resultant, calculated as the sum of the residuals at each node of the cranial boundary, were saved, resulting in a simulated compressive load-displacement curve. The apparent vertebral stiffness (K) was calculated as the slope of the curve linear region. Failure strength (F max ) was computed as the maximal compressive load predicted. In addition, the spatial distribution and magnitude of local stresses, strains, and simulated damage magnitude were computed. The simulation was set to stop when the local damage D of any quadrature point reached 1 (completely damaged). The model damage simulations at failure corresponding to the nine vertebrae for which load-displacement curves were illustrated in Figure 3 are presented in Supplementary Figure SA1 .

2.3.3 Mesh resolution study

The boundary value problem was solved initially for the hFE 1 mm mesh model, a mesh size selected based on our previous studies ( Stadelmann et al., 2020 ). To study the effect of mesh size on the results, we produced additional FE mesh models at 0.5 mm and 2 mm for each vertebral sample. The hFE model was employed for each vertebral sample to simulate vertebral strength and stiffness at each mesh size. Figure 3 shows the load-displacement curves of 10 sample cadaver vertebral bodies with 2, 1, and 0.5 mm mesh sizes. Quantitative values of vertebral strength and stiffness at each mesh size are detailed in Supplementary Table SA1 .

2.4 Vertebral specific material model parameter calibration

Similar to the previous section, we initially solved the boundary value problem for the hFE model 1 mm mesh model with the initial reference viscoelastic, viscoplastic cortical bone, and damage model based on material model parameters for bone ( Johnson et al., 2010 ; Dall’Ara et al., 2013 ). We used an iterative procedure where the most relevant material parameters were incrementally adjusted from their reference values in the direction in which they would reduce the RMSE of the simulated load-displacement curve relative to the experimentally measured one ( Stadelmann et al., 2020 ) to less than 5%. For complete details, see Supplementary Appendix . The resulting material model parameters are reported in Table 2 .

2.5 Optimization method for computing “global” material parameters

The vertebra’s geometry and internal microstructure affect its ability to withstand loads, specifically its stiffness and strength values. Each of the study’s ten vertebrae varies in geometry and internal microstructure parameters due to bone metastasis and individual variation in geometry and material properties. We aimed to investigate whether a single set of material parameters could be applied to simulate vertebral strength and stiffness across the three bone metastasis types. Such a model is highly desirable as it will obviate the need to achieve patient-specific calibration. To achieve this goal, we define an objective function as: O b j = λ × R M S E + 1 − λ × S D (19) where RMSE (Root Mean Squared Error) represents the average error across all samples, and SD denotes the standard deviation as a measure of the variance or spread of errors across all the samples. λ is a hyperparameter that controls the trade-off between the two objectives, here set as = 0.5. The normalized error between predicted and experimental mechanical responses for each vertebra is defined below and obtained using the trapezoidal rule, which approximates the area under the curve by dividing it into trapezoids and summing their areas, Eq. 19 : e i = ∫ 0 δ f a i l u r e i δ F ^ i d δ (20) where δ F ^ i = F h F E i − F exp i / F exp i is the normalized error, and δ f a i l u r e i is the displacement at the failure point for the i th vertebra. Normalizing the error (dividing the error for each vertebra by the area under the experimental force-displacement curve of that vertebra) ensures consistency in comparing errors across different experimental data sets. This process is particularly important as we process heterogeneous data sources with different stiffness and strength value ranges. The material model parameters are selected within the material parameters obtained for each vertebra. This strategy ensures that the initial parameter values are within a feasible range based on the available data while allowing for exploration of the parameter space to find a common set of parameters applicable to all vertebrae.

We performed the FE analysis for all ten vertebrae using the selected material model parameters, calculated the error between the FE predictions and experimental measurements using Eq. 20 , and evaluated the selected set of material model parameters performance by calculating the objective function Eq. 19 . We used the pair of material parameters and the calculated objective function to train a surrogate model, with the surrogate model implemented using a three-layer neural network framework provided by PyTorch having two hidden layers containing 20 neurons and 10 neurons, respectively. Using the trained surrogate model computed, we applied the gradient descent method to optimize the material parameters to minimize the objective function. Using the trained surrogate model for an initial guess, we computed the gradient of the objective function concerning the material parameters and adjusted the parameters iteratively in the direction that minimizes the objective function. This process was performed using PyTorch’s automatic differentiation capabilities, efficiently computing gradients with respect to the material parameters instead of the trained surrogate model’s internal parameters (weights and biases). Table 2 presents the obtained values.

2.6 Statistical analysis

Statistical analysis was performed in JMP Pro (14.3, SAS Institute, Inc.). Descriptive statistics were used to summarize the specimens’ demographic, bone, and mechanical properties. The strength and stiffness data were not normally distributed based on the normality test results (Shapiro-Wilk test). Based on the test for unequal variability between the mesh size data, Welch’s ANOVA was used to test for the effect of 1) bone metastases on the experimental strength and stiffness of the osteolytic and osteosclerotic vertebrae, 2) mesh size on the differences in strength and stiffness and 3) computed errors with respect to the experimental data. Post hoc comparisons were performed using the non-parametric Steel-Dwass test.

The sampling of multiple vertebrae per spine can introduce clustering (non-independence) of the data. We fitted linear mixed-effects models (LMMs) under different assumptions about the correlation structure among segments from the same spines to test this effect. According to the Akaike information criterion (AIC), the independence structure best fits the data ( O’Brien and Fitzmaurice, 2004 ). We applied linear regression to test the association between the simulated and measured strength and stiffness based on this finding. Statistical significance was set at the 5% level.

3 Results

3.1 Experimental study

Vertebral demographics, primary cancer, bone material, and vertebral mechanical properties are listed in Table 1 . Corresponding experimental load-displacement curves are presented in Figure 2 . Osteolytic vertebrae strength, a mean (standard deviation) of 2.65 (1.09) kN, and stiffness 10.75 (1.51) kN/mm were lower than osteosclerotic vertebrae [7.44 (2.43) kN, p = 0.0069 and 14.45 (2.96) kN/mm, p = 0.0405, respectively]. The vertebrae with mixed bone metastases exhibited high strength, 10.80 kN, and stiffness, 36.84 kN/mm.

3.2 Bone metastasis affects vertebral structure

Presents sagittal and axial μ CT image data for an osteolytic, mixed, and osteosclerotic vertebral sample used in the study. The osteolytic vertebrae presented a thinning of the cortex and rarefaction of the trabeculae ( Figure 1 : M1 and M2) with focal bone loss (cavities) through the entire Vertebra ( Figure 1 ). Osteosclerotic vertebrae demonstrate a marked increase in BV/TV, as can be observed from the near solid tissue ( Figure 1 : P2 and P2). Per its classification, Vertebra with a mixed lesion demonstrates regions of high bone density with regions containing lytic lesions ( Figure 1 : V1 and V2). Descriptive analysis found osteosclerotic vertebrae with higher bone mineral density (BMD), mean (standard deviation) value of 213.19 (60.44) mgHA/cm 3 and bone volume fraction (BV/TV), 12.80 (7.13)%, than osteolytic vertebrae, (108.095 (23.75) mgHA/cm 3 and 9.06 (0.70)%, respectively). ANOVA analysis found these differences not statically significant at the 5% level. The mixed lesion vertebrae showed high BMD (324.68) mgHA/cm 3 and BV/TV (18.61)%, a value significantly higher than osteolytic vertebrae BMD ( p = 0.039) Figure 2 .

3.3 Mesh resolution study

Figure 3 presents experimental and predicted load-displacement curves at 2, 1, and 0.5 mm mesh models for MD15A-L1. Corresponding mid-sagittal von Mises stresses at vertebral strength are presented for each mesh size simulation. Supplementary Tables SA1, SA2 detail the specimen-specific strength, stiffness, and corresponding error (%) values at 2, 1, and 0.5 mm mesh models.

3.3.1 Strength simulations

FE-simulated strength, a mean (standard deviation) of 6.82 (3.16) kN, showed excellent agreement with the measured strength, 6.80 (3.16) kN, Supplementary Table SA1 , simulation error compared to the experimental strength, Root Mean Squared Error (RMSE) = 0.04. Mesh refinement, 0.5 mm, underestimated strength, 5.86 (2.81) kN/mm Supplementary Table SA1 , with higher simulation error (RMSE = 0.42). A coarser mesh size, 2 mm, overestimated strength, 7.33 (3.48), Supplementary Table SA1 , (RMSE = 0.32). ANOVA analysis found these differences significant ( p < 0.0001). Post-test comparisons showed 0.5 mm and 2 mm RMSE values significantly higher than the 1 mm model ( p = 0.0005 and p = 0.0017, respectively), with 0.5 mm model RMSE higher than the 2 mm model ( p = 0.0160).

3.3.2 Stiffness simulations

The predicted stiffness at 1 mm mesh model, 18.50 (11.47) kN/mm, closely agreed with the experimental values, 18.60 (11.23) kN/mm, Supplementary Table SA2 , with RMSE = 0.23. Either refinement, to 0.5 mm or coarsening, to 2 mm, of the mesh models resulted in higher predicted stiffness [19.48 (11.23), RMSE = 0.75 and 19.86 (11.79), RMSE = 1.21, respectively] kN/mm, Supplementary Table SA2 . ANOVA analysis found these differences significant ( p < 0.0192). Post-test analysis showed the error values for the 2 mm mesh higher than the 1 mm model ( p = 0.0101).

3.4 hFE predicts metastatic vertebrae strength and stiffness

Figure 4 presents the mean load-displacement curve for experimental and simulated specimen-specific ( Figure 4A ) and global ( Figure 4B ) experiments. Detailed values for specimen-specific and global material calibration-based prediction of strength and stiffness values and corresponding computed error concerning experimental values are detailed in Supplementary Tables SA3, SA4 , respectively.

3.4.1 Specimen-specific material calibration

Regression analysis showed specimen-specific material calibration simulations explained 99% of the variance of measured strength ( p < 0.0001, Figure 4A .I) and stiffness ( p < 0.0001, Figure 4A .II). The rheological model ( Figure 5 ) showed the osteosclerotic vertebrae long-term modulus, E 0 0 , [3.5 (0.92)]GPa, viscoelastic low-, E 0 1 : [0.96 (0.25)]GPa and high-, E 0 2 : [5.19 (1.35)]GPa modulus, and viscoelastic strain rate low-, η 0 1 : [2.89 (0.76)MPa∙s, and high, η 0 2 : 0.49 (0.13)kPa∙s] parameters lower than osteolytic vertebrae elastic [ E 0 0 : 4.2 (0.28), E 0 1 : 1.15 (0.07), and E 0 2 : 6.09 (0.41)]GPa and viscoelastic [ η 0 1 : 3.43 (0.24)MPa∙s and η 0 2 :0.59 (0.04)kPa∙s] parameters, respectively, Supplementary Table SA5 . Mixed-lesion vertebra model parameters were comparable to osteolytic vertebrae.

3.4.2 Global material calibration

Regression analysis demonstrated that the simulation using global material calibration explained 83% of the variance in measured strength ( p = 0.0002) and 92% in stiffness ( p < 0.0001), Figure 4B . Compared to the specimen-specific calibration, analysis of variance (ANOVA) found the “global” model reduction in strength prediction accuracy to be significant ( p = 0.0060). We found no such statistically significant difference in the stiffness prediction at the 5% level.

3.5 FE simulation suggests bone metastasis type affects bone damage evolution differentially

In the osteolytic vertebra, the FE simulation suggested the damage initiating at yield strength at the posterior cortex ( Figure 6A .I), evolving within the vertebral body through a large region of low bone modulus to failure ( Figure 6A .II). In the osteosclerotic vertebra, the FE simulation suggested damage initiated at yield strength at the posterior cortex ( Figure 6B .I), the damage evolving largely confined to regions of low bone modulus principally bypassing regions of high bone modulus ( Figure 6B .II). Upon reaching simulated failure, the simulated damage was confined to the pre-damaged region in the osteolytic vertebra ( Figure 6A .III,IV) while, in the osteosclerotic vertebrae, the simulated damage evolved within bone regions with low and heterogeneous distribution of bone modulus ( Figure 6B .III,IV).

4 Discussion

This study developed a constitutive, damage-based finite element model to simulate the strength and stiffness of cadaveric pathologic vertebrae containing osteolytic, osteosclerotic, and mixed bone lesions. Derived from CT data comparable to standard clinical CT image resolution and implemented within the ΣMIT computational solid mechanics framework ( Radovitzky et al., 2011 ), the model accurately predicted the vertebrae’s strength and stiffness under applied compressive loading independent of the metastatic bone lesion type. We further demonstrated the feasibility of a machine learning approach to establish a unified set of material calibration parameters, providing a strong prediction of pathologic vertebral strength and stiffness. This finding suggests that such a unified set, once validated in a larger data set, might reduce the need to incorporate density phantoms during patient radiation treatment planning, which would require significant investment in the recertification of CT protocols. The FE model damage simulation suggested that osteolytic and osteosclerotic lesions differentially affect damage evolution within the bone network.

4.1 Selection of vertebral specimens for hFE model validation

We aimed to simulate the measured compressive strength and stiffness of human vertebrae with osteolytic, osteosclerotic, and mixed bone metastases from a group of 45 cadaveric vertebrae that were obtained from donors with breast, lung, renal, and prostate cancers and mechanically tested in our previous study ( Stadelmann et al., 2020 ). The vertebrae selected demonstrated either markedly higher strength, stiffness, and BMD values, for example, VA15-S-T11 [(+84.5, +49.8, and 73.1)%, respectively), markedly lower, [GA1-T11: (−69.7, −60.7 and −50.4)% respectively] or asynchronous values for mechanical and material properties, [PA15-T8 (−13.8, −59.4, +30.8)%, respectively], than the median strength and stiffness values measured for the 45 vertebrae tested ( Stadelmann et al., 2020 ). Although little data exists on the strength of human vertebrae with osteosclerotic or mixed bone metastases, the current study osteolytic vertebral strength, a mean (SD) of 2.6 (1.1)kN ( Stadelmann et al., 2020 ), agrees with strength values reported for osteolytic cadaver vertebrae by Costa et al. (2019) . Therefore, the pathologic vertebrae selected for this study form a challenging sample for validating the FE framework proposed.

4.2 Establishing the hFE modeling approach

We postulated that integrating bone tissue viscoelastic and viscoelastic response may improve the computational damage model’s ability to simulate the bone metastases-mediated changes in vertebral strength and stiffness. For this study, we modified the elastic element in each of the viscoelastic model branches to incorporate 1) the damage evolution model ( Dall’Ara et al., 2013 ) to capture the spatial bone structure material softening ( Section 4.2 ) and 2) the nonlinear relationships between trabecular bone material properties ( Section 4.1 ) introduced via bone fraction value (BV/TV) ( Zysset, 2003 ) within the model’s elasticity-density relationship for trabecular bone tissue.

Based on vertebral-specific material calibration, the resulting computational framework successfully explained 99% of the variation in the vertebrae measured strength and stiffness invariant to the type of bone metastases. This prediction improves upon our previous hFE model ( Stadelmann et al., 2020 ) employing “state of the art” FE modeling of human vertebrae mechanical behavior ( Pahr et al., 2014 ), which explained 71% of the variance in the current specimens’ measured strength ( p = 0.0023) and 52% of measured stiffness ( p = 0.0180). This confirms that the FE model, as implemented, captures the factors that determine strength and stiffness.

4.2.1 Mesh size affects predicted vertebral strength and stiffness

In agreement with Jones and Wilcox (2007) , mesh size affected our model strength and stiffness simulation accuracy. Regression analysis showed the 1 mm calibrated mesh size to yield an excellent agreement between the predicted strength, stiffness, and measured values. Higher mesh size, 2 mm, affected higher averaging of the bone microstructure and thickening of the vertebral cortex region, likely to yield higher bone density per element. This increase resulted in the model’s overestimated simulation strength, a mean of 4.1%, and degraded accuracy ( Supplementary Table SA1 ). Refining mesh size to 0.5 mm yields higher structural fidelity for bone microstructure and vertebral cortex anatomical detail. Although the finer mesh model is expected to produce improved prediction, the higher mesh model underestimated measured strength, a mean of −17.8%, with significantly lower simulation accuracy, Supplementary Table SA1 . We assume that the finer mesh contains stiffer and softer elements due to the broader range of local bone volume densities leading to “localized” damage accumulation within the finer mesh regions of low density. This “localized” damage affects higher numerical instability within the mesh, which “softens” the simulated stress-strain response, a finding in agreement with Werner et al. (2019) reporting FE mesh refinement to soften the post-yield behavior of trabecular bone under large deformation. Given these findings, we elected to perform our simulations at 1 mm mesh element size.

4.2.2 Optimization approach for establishing FE model parameters

Capturing the spatial variation in the bone tissue density-elasticity relationships requires mapping CT Hounsfield unit (HU) data to bone density values using a calibrated bone density phantom. However, introducing the density phantom in the imaging field will require re-certifying the radiotherapy planning and imaging protocols for each CT imaging and clinical linear accelerator system, a highly costly and complicated process. Phantomless calibration is based on skeletal muscle or adipose tissues’ known linear attenuation properties to create a standard calibration curve for converting trabecular HU to vBMD the HU values ( Gudmundsdottir et al., 1993 ). Although this technique strongly correlates with phantom-based vBMD analysis in non-cancer patients, it has not been validated in cancer patients, and its utility has yet to be shown to be generalized on multiple CT scanners. Asynchronous calibrations in which the bone-density phantom-based calibration is performed daily as part of the CT system quality assurance separate from the patient scan may alleviate some of these challenges. We, therefore, investigated whether establishing a single set of “optimized” material parameters based on initial asynchronous calibrations using a bone density phantom could simulate vertebral strength and stiffness across the wide range of bone metastases tested in this study.

For this purpose, we applied a “Surrogate optimization” approach ( Espinosa et al., 2023 ) to find a single set of parameters to minimize the summation of errors across all samples while ensuring that the error for each sample is within a similar range. This resulting common objective creates an approximation model that mimics the behavior of the true objective function when the objective function is computationally expensive, time-consuming, or difficult to evaluate directly, allowing for faster parameter space exploration. For this purpose, we used neural networks to construct the surrogate model with the optimization algorithm leveraging gradient-based optimization techniques to find a single set of material model parameters that minimizes error for all ten vertebrae while ensuring that the error is not getting too large. The optimized “global” model explained 83% of the measured strength variance ( p = 0.0002) and 92% of the measured stiffness variance ( p < 0.0001).

A possible explanation for the model’s stronger prediction for stiffness vs. strength might be associated with the power coefficient k (Eq. 4 ) assumed to be the same for the model’s morphology-elasticity and morphology-yield relationships. Although stiffness, yield, and ultimate strength of trabecular bone were reported to be strongly correlated ( Goulet et al., 1994 ), it may be necessary to define k differently for the two relationships. Overall, this work supports the promise of the global parameter approach to simulate vertebral body mechanics at an actionable level of confidence.

4.3 Bone metastases phenotype affects the pattern of vertebral bone damage

The effect of bone metastases type on bone damage accumulation underlying the failure and post-failure behavior of pathologic vertebrae is poorly understood. Atkins et al. (2019) reported that bone metastasis type (osteolytic vs. mixed) affects the pattern of bone tissue microdamage, with pathologic bone exhibiting a significantly higher magnitude of diffuse microdamage compared to the healthy bone for a comparable compressive load. Micro-FE simulations ( Choudhari et al., 2016 ; Atkins et al., 2019 ) found that regions of higher stress and strain within the pathologic bone correlated with the extent and pattern of diffused micro-damage, the damaged bone tissue showing significantly higher stress magnitude than healthy bone tissue for comparable compressive loading ( Atkins et al., 2019 ). Our FE simulation demonstrates that bone metastases type produces a distinct difference in damage evolution within the vertebral body and that, for both bone metastases types, the damage evolved localized to regions of low bone modulus ( Figure 6 ). Our results are in qualitative agreement with these studies.

The FE framework rheological model revealed osteolytic bone with higher elastic and viscoelastic modulus values than osteosclerotic bone ( Figure 5 ), which agrees with our previous experimental study ( Stadelmann et al., 2020 ). Burke et al. (2016 , 2017) reported osteolytic bone metastases alter the bone’s collagen fibril organization and degrade tissue mineralization, likely to increase collagen fibrillar sliding at the collagen/mineral interface under applied stress, affecting the bone plasticity resulting in higher bone microdamage ( Vashishth et al., 2000 ). Our simulations suggest that in osteolytic bone, damage accumulation is primarily associated with inelastic strain accumulation and viscous property, forming new internal surfaces and voids in the tissue ( Arthur Moore and Gibson, 2002 ), affecting the magnitude of diffused damage in the bone.

In osteosclerotic bone, the FE framework rheological model predicted the bone tissue with lower elastic and viscoelastic modulus despite the predicted higher compressive and shear strength at failure. Although little is known about the mechanism of osteosclerotic bone failure, higher bone mineralization, bone connectivity, and thicker bone trabecula yield higher bone volume to tissue volume ( Bailey et al., 2022 ), resulting in a plate-like, dense bone network, which may explain the rheological model prediction of higher strength at the tissue level. The simulation suggested that osteosclerotic vertebrae’s stiffness and overall failure strength were dominated by damage evolution localized to regions of low-bone modulus within the bone network, largely bypassing regions of high modulus associated with high BV/TV. This finding supports our experimental study ( Alkalay et al., 2021b ), which shows that osteosclerotic and mixed (comprising osteosclerotic and osteolytic regions) vertebrae exhibit similar measured stiffness values. Once failed, the simulation suggested that for both bone metastases types, evolution of bone damage appears localized to the low bone modulus regions within the vertebral volume ( Figure 6 ). Additional studies are required to understand the impact of pathologic bone tissue structure and quality on the damage accumulation process and, ultimately, the failure of pathologic vertebrae.

4.4 Study limitations

This study has several limitations beyond the limitations of cadaveric spines.

4.4.1 Experimental protocol

Although our sample size of 10 human vertebrae is small, the vertebral samples were obtained from donors with known spine metastases containing a wide range of clinically observed bone metastases. Our study ( Stadelmann et al., 2020 ) used isolated vertebral bodies with the endplate and posterior elements removed to permit imaging in a small-bore high-resolution CT device and provide a standardized specimen geometry. Computational analysis ( Maquer et al., 2014 ) found that removing the vertebral endplates had a low impact on the vertebrae’s predicted strength and overall damage distribution. However, the degenerative state of the intervertebral disc ( Palanca et al., 2023 ) and the mechanical interaction of the posterior elements ( Alkalay et al., 2021 ) affect the strength, stiffness, mechanical instability and failure patterns ( Alkalay et al., 2021a ) of pathologic spines. Our current ongoing work aims to incorporate a newly designed computer control mechanical test system designed to allow time-lapse imaging under controlled loading conditions within the XtremeCT II (SCANCO Medical AG, Switzerland) to provide detailed measurements of the effects of bone metastasis on the spine deformation and failure.

Our study used a quasi-static monotonic mechanical loading protocol ( Stadelmann et al., 2020 ). However, monotonic loading has limited fidelity concerning physiological loading conditions. Importantly, the interaction of these loads with the poroelastic intervertebral disc joints ( Newell et al., 2020 ) and, through this interaction, the deformation of the vertebral endplates, will modulate and alter the stress/strain within the affected vertebrae ( Alkalay and Harrigan, 2016 ; Palanca et al., 2023 ), thereby affecting vertebral failure ( Jackman et al., 2016 ). Recapitulating dynamic testing conditions in vitro remains a significant challenge ( Costi et al., 2021 ). Hence, our study offers a narrow simulation of these effects.

As discussed in our previous study ( Stadelmann et al., 2020 ), we could not produce perfectly parallel sections during sample preparation. As a result, we could not produce “zero stress” conditions at the initiation of the test affecting the measurements of vertebral stiffness. Although this does not affect the resulting experimental strength or the model simulation of strength, our study reported experimental and simulated values for vertebral stiffness represent “apparent” vertebral stiffness ( Burghardt et al., 2007 ).

4.4.2 Modeling approach

Although the FE damage-based simulation strongly predicted the experimental measures of vertebral stiffness and strength, we did not directly validate the model’s damage predictions. We sought to predict vertebral strength and stiffness using our damage model as a modeling approach. The accuracy of the prediction validates the approach without explicit measurement of damage, which was not necessary for our study. Previous studies employed cyclic loading experiments to evaluate bone damage ( Zysset and Curnier, 1996 ), yielding loading boundary conditions significantly different from our simulations. We are evaluating controlled loading experiments within high-resolution CT, which, combined with image-based volumetric registration, allow computation of the pattern of internal bone strains to validate the FE damage predictions.

Our FE model employed an isotropic bone volume density parameter to compute the morphology-elasticity relationship. Although bone volume density accounts for 84%–87% of the variation in bone stiffness and strength ( Maquer et al., 2015 ), fabric anisotropy was found to explain a small but appreciable (∼10%) amount of the variation in the elastic properties of bone ( Maquer et al., 2015 ). Given the variation in the location, geometry, and character of metastases, including anisotropy in our constitutive model may further improve the predictability of our “optimization-based” models.

4.5 Clinical implications

Vertebral fractures can be devastating events for patients with spinal metastases. Once spinal metastases are identified, there are potential prophylactic therapies that may reduce PVF risk. All are invasive, and treating physicians seek to avoid them unless necessary. Attempts to generate individualized fracture risk predictions have had limited success thus far ( Fourney et al., 2011 ). Although our sample set was small, it includes vertebrae impacted by all metastatic lesion types demonstrating a wide range of vertebral stiffnesses and strengths, permitting characterization of the load at which they fail and, critically, localizing the vertebral region of initial failure. This information, representing an unmet clinical need, may guide localized interventions such as vertebral augmentation to the portions of the vertebral body at greatest risk upon further development and validation. In many patients, the affected vertebral bodies will have lost height due to failure that may be subclinical at presentation. In these cases, methods developed for predicting the fracture risk of intact bodies cannot be expected to succeed. Our approach holds promise for estimating residual strength after the initial fracture. This question of residual strength applies to both pathologic and osteoporotic fractures.

Data availability statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Ethics statement

Ethical approval was not required for the studies involving humans because Cadaver human samples are exempt based on NIH rules. The studies were conducted in accordance with the local legislation and institutional requirements. The human samples used in this study were acquired from all cadaveric spines were purchased with NIH grant funding from the Anatomical Gift Registry, Hanover, MD, United States. Written informed consent to participate in this study was not required from the participants or the participants’ legal guardians/next of kin in accordance with the national legislation and the institutional requirements.

Author contributions

ZS: Data curation, Formal Analysis, Investigation, Methodology, Software, Validation, Visualization, Writing–review and editing. MX: Data curation, Formal Analysis, Investigation, Methodology, Software, Validation, Visualization, Writing–original draft, Writing–review and editing. RR: Conceptualization, Data curation, Formal Analysis, Investigation, Methodology, Project administration, Software, Supervision, Validation, Writing–original draft, Writing–review and editing. MS: Data curation, Formal Analysis, Methodology, Software, Validation, Writing–review and editing. DH: Conceptualization, Formal Analysis, Funding acquisition, Resources, Writing–original draft, Writing–review and editing. RA: Conceptualization, Data curation, Formal Analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Supervision, Validation, Visualization, Writing–original draft, Writing–review and editing.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The author(s) declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbioe.2024.1424553/full#supplementary-material

hFE model simulated damage accumulation pattern at failure for the study’s remaining osteosclerotic, osteolytic, and mixed lesion vertebrae. For ease of visualization, the medial half of the vertebral body was made transparent to view the spatial evolution of damage within the vertebral body. The lateral half was rendered solid with blue, indicating no damage at a specific location. The marked differences in the damage pattern of the osteosclerotic vertebrae reflect the high variation observed for their load-displacement response. By comparison, the osteolytic vertebra is simulated to fail at a lower value of damage accumulation. The mixed lesion vertebra was simulated to fail at a high level of damage, which may contribute to its measured high strength and stiffness.

References

Alkalay R. N. Adamson R. Miropolsky A. Davis R. B. Groff M. L. Hackney D. B. ( 2021a ). Large lytic defects produce kinematic instability and loss of compressive strength in human spines: an in vitro study . J. Bone Jt. Surg. Am. 103 ( 10 ), 887 – 899 . 10.2106/jbjs.19.00419 PMC9167060 33755638

Alkalay R. N. Groff M. W. Stadelmann M. A. Buck F. M. Hoppe S. Theumann N. ( 2021b ). Improved estimates of strength and stiffness in pathologic vertebrae with bone metastases using CT-derived bone density compared with radiographic bone lesion quality classification . J. Neurosurg. Spine. 36 ( 1 ), 113 – 124 . 10.3171/2021.2.SPINE202027 34479191 PMC9210826

Alkalay R. N. Harrigan T. P. ( 2016 ). Mechanical assessment of the effects of metastatic lytic defect on the structural response of human thoracolumbar spine . J. Orthop. Res. 34 , 1808 – 1819 . 10.1002/jor.23154 26748564

Arthur Moore T. L. Gibson L. J. ( 2002 ). Microdamage accumulation in bovine trabecular bone in uniaxial compression . J. Biomech. Eng. 124 ( 1 ), 63 – 71 . 10.1115/1.1428745 11873773

Atkins A. Burke M. Samiezadeh S. Akens M. K. Hardisty M. Whyne C. M. ( 2019 ). Elevated microdamage spatially correlates with stress in metastatic vertebrae . Ann. Biomed. Eng. 47 ( 4 ), 980 – 989 . 10.1007/s10439-018-02188-8 30673956

Bailey S. Hackney D. Vashishth D. Alkalay R. N. ( 2020 ). The effects of metastatic lesion on the structural determinants of bone: current clinical and experimental approaches . Bone 138 , 115159 . 10.1016/j.bone.2019.115159 31759204 PMC7531290

Bailey S. Stadelmann M. A. Zysset P. K. Vashishth D. Alkalay R. N. ( 2022 ). Influence of metastatic bone lesion type and tumor origin on human vertebral bone architecture, matrix quality, and mechanical properties . J. Bone Min. Res. 37 ( 5 ), 896 – 907 . 10.1002/jbmr.4539 PMC9158727 35253282

Begg A. C. Stewart F. A. Vens C. ( 2011 ). Strategies to improve radiotherapy with targeted drugs . Nat. Rev. Cancer 11 ( 4 ), 239 – 253 . 10.1038/nrc3007 21430696

Burghardt A. J. Kazakia G. J. Majumdar S. ( 2007 ). A local adaptive threshold strategy for high resolution peripheral quantitative computed tomography of trabecular bone . Ann. Biomed. Eng. 35 ( 10 ), 1678 – 1686 . 10.1007/s10439-007-9344-4 17602299

Burke M. Akens M. Kiss A. Willett T. Whyne C. ( 2018 ). Mechanical behavior of metastatic vertebrae are influenced by tissue architecture, mineral content, and organic feature alterations . J. Orthop. Res. 36 ( 11 ), 3013 – 3022 . 10.1002/jor.24105 29978906

Burke M. Atkins A. Kiss A. Akens M. Yee A. Whyne C. ( 2017 ). The impact of metastasis on the mineral phase of vertebral bone tissue . J. Mech. Behav. Biomed. Mater 69 , 75 – 84 . 10.1016/j.jmbbm.2016.12.017 28040610

Burke M. Golaraei A. Atkins A. Akens M. Barzda V. Whyne C. ( 2017 ). Collagen fibril organization within rat vertebral bone modified with metastatic involvement . J. Struct. Biol. 199 ( 2 ), 153 – 164 . 10.1016/j.jsb.2017.06.008 28655593

Burke M. V. Atkins A. Akens M. Willett T. L. Whyne C. M. ( 2016 ). Osteolytic and mixed cancer metastasis modulates collagen and mineral parameters within rat vertebral bone matrix . J. Orthop. Res. 34 ( 12 ), 2126 – 2136 . 10.1002/jor.23248 27027407

Campbell G. M. Pena J. A. Giravent S. Thomsen F. Damm T. Gluer C. C. ( 2017 ). Assessment of bone fragility in patients with multiple myeloma using QCT-based finite element modeling . J. Bone Min. Res. 32 ( 1 ), 151 – 156 . 10.1002/jbmr.2924 27454865

Cavuoto R. Lenarda P. Tampieri A. Bigoni D. Paggi M. ( 2024 ). Phase-field modelling of failure in ceramics with multiscale porosity . Mater. Des. 238 , 112708 . 10.1016/j.matdes.2024.112708

Chevalier Y. Charlebois M. Pahr D. Varga P. Heini P. Schneider E. ( 2008 ). A patient-specific finite element methodology to predict damage accumulation in vertebral bodies under axial compression, sagittal flexion and combined loads . Comput. Methods Biomech. Biomed. Engin 11 ( 5 ), 477 – 487 . 10.1080/10255840802078022 18608338

Chevalier Y. Pahr D. Zysset P. K. ( 2009 ). The role of cortical shell and trabecular fabric in finite element analysis of the human vertebral body . J. Biomechanical Eng. 131 , 111003 – 111012 . 10.1115/1.3212097 20353254

Chevalier Y. Pahr D. Zysset P. K. ( 2009 ). The role of cortical shell and trabecular fabric in finite element analysis of the human vertebral body . J. Biomech. Eng. 131 ( 11 ), 111003 . 10.1115/1.3212097 20353254

Choudhari C. Chan K. Akens M. K. Whyne C. M. ( 2016 ). μFE models can represent microdamaged regions of healthy and metastatically involved whole vertebrae identified through histology and contrast enhanced μCT imaging . J. Biomech. 49 ( 7 ), 1103 – 1110 . 10.1016/j.jbiomech.2016.02.034 26947031

Costa M. C. Eltes P. Lazary A. Varga P. P. Viceconti M. Dall'Ara E. ( 2019 ). Biomechanical assessment of vertebrae with lytic metastases with subject-specific finite element models . J. Mech. Behav. Biomed. Mater. 98 , 268 – 290 . 10.1016/j.jmbbm.2019.06.027 31280054

Costi J. J. Ledet E. H. O'Connell G. D. ( 2021 ). Spine biomechanical testing methodologies: the controversy of consensus vs scientific evidence . JOR Spine 4 ( 1 ), e1138 . 10.1002/jsp2.1138 33778410 PMC7984003

Cronin K. A. Lake A. J. Scott S. Sherman R. L. Noone A. M. Howlader N. ( 2018 ). Annual report to the nation on the status of cancer, part I: national cancer statistics . Cancer 124 ( 13 ), 2785 – 2800 . 10.1002/cncr.31551 29786848 PMC6033186

Dall'Ara E. Luisier B. Schmidt R. Kainberger F. Zysset P. Pahr D. ( 2013 ). A nonlinear QCT-based finite element model validation study for the human femur tested in two configurations in vitro . Bone 52 , 27 – 38 . 10.1016/j.bone.2012.09.006 22985891

Dall'Ara E. Schmidt R. Pahr D. Varga P. Chevalier Y. Patsch J. ( 2010 ). A nonlinear finite element model validation study based on a novel experimental technique for inducing anterior wedge-shape fractures in human vertebral bodies in vitro . J. Biomech. 43 ( 12 ), 2374 – 2380 . 10.1016/j.jbiomech.2010.04.023 20462582

Ding M. Odgaard A. Hvid I. ( 1999 ). Accuracy of cancellous bone volume fraction measured by micro-CT scanning . J. Biomechanics 32 ( 3 ), 323 – 326 . 10.1016/s0021-9290(98)00176-6 10093033

Espinosa B. O. U. Quijada J. G. P. Kurkina E. Lukyanov O. ( 2023 ). Surrogate aerodynamic wing modeling based on a multilayer perceptron . Aerospace 10 ( 2 ), 149 . 10.3390/aerospace10020149

Faruqi S. Tseng C. L. Whyne C. Alghamdi M. Wilson J. Myrehaug S. ( 2018 ). Vertebral compression fracture after spine stereotactic body radiation therapy: a review of the pathophysiology and risk factors . Neurosurgery 83 ( 3 ), 314 – 322 . 10.1093/neuros/nyx493 29048517

Fourney D. R. Frangou E. M. Ryken T. C. Dipaola C. P. Shaffrey C. I. Berven S. H. ( 2011 ). Spinal instability neoplastic score: an analysis of reliability and validity from the spine oncology study group . J. Clin. Oncol. 29 ( 22 ), 3072 – 3077 . 10.1200/jco.2010.34.3897 21709187

Garcia D. Zysset P. K. Charlebois M. Curnier A. ( 2009 ). A three-dimensional elastic plastic damage constitutive law for bone tissue . Biomechanics Model. Mechanobiol. 8 , 149 – 165 . 10.1007/s10237-008-0125-2 18398628

Goulet R. W. Goldstein S. A. Ciarelli M. J. Kuhn J. L. Brown M. B. Feldkamp L. A. ( 1994 ). The relationship between the structural and orthogonal compressive properties of trabecular bone . J. Biomech. 27 ( 4 ), 375 – 389 . 10.1016/0021-9290(94)90014-0 8188719

Gudmundsdottir H. Jonsdottir B. Kristinsson S. Johannesson A. Goodenough D. Sigurdsson G. ( 1993 ). Vertebral bone density in Icelandic women using quantitative computed tomography without an external reference phantom . Osteoporos. Int. 3 ( 2 ), 84 – 89 . 10.1007/bf01623378 8453195

Jackman T. M. DelMonaco A. M. Morgan E. F. ( 2016 ). Accuracy of finite element analyses of CT scans in predictions of vertebral failure patterns under axial compression and anterior flexion . J. Biomechanics 49 , 267 – 275 . 10.1016/j.jbiomech.2015.12.004 PMC4955561 26792288

Johnson T. P. M. Socrate S. Boyce M. C. ( 2010 ). A viscoelastic, viscoplastic model of cortical bone valid at low and high strain rates . Acta Biomater. 6 , 4073 – 4080 . 10.1016/j.actbio.2010.04.017 20417735

Jones A. C. Wilcox R. K. ( 2007 ). Assessment of factors influencing finite element vertebral model predictions . J. Biomech. Eng. 129 ( 6 ), 898 – 903 . 10.1115/1.2800791 18067394

Maquer G. Musy S. N. Wandel J. Gross T. Zysset P. K. ( 2015 ). Bone volume fraction and fabric anisotropy are better determinants of trabecular bone stiffness than other morphological variables . J. Bone Mineral Res. 30 , 1000 – 1008 . 10.1002/jbmr.2437 25529534

Maquer G. Schwiedrzik J. Zysset P. K. ( 2014 ). Embedding of human vertebral bodies leads to higher ultimate load and altered damage localisation under axial compression . Comput. Methods Biomech. Biomed. Engin 17 ( 12 ), 1311 – 1322 . 10.1080/10255842.2012.744400 23237518

Mauch S. Noels L. Zhao Z. Radovitzky R. ( 2006 ). Lagrangian simulation of penetration environments via mesh healing and adaptive optimization . 25th U. S. Army Sci. Conf. , 1 – 3 .

McDougall J. A. Bansal A. Goulart B. H. McCune J. S. Karnopp A. Fedorenko C. ( 2016 ). The clinical and economic impacts of skeletal-related events among medicare enrollees with prostate cancer metastatic to bone . Oncologist 21 ( 3 ), 320 – 326 . 10.1634/theoncologist.2015-0327 26865591 PMC4786354

Molinari L. Falcinelli C. ( 2022 ). On the human vertebra computational modeling: a literature review . Meccanica 57 ( 3 ), 599 – 622 . 10.1007/s11012-021-01452-x

Nazarian A. Stechow D. V. Zurakowski D. Müller R. Snyder B. D. ( 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 . 10.1007/s00223-008-9174-x 18946628

Newell N. Rivera Tapia D. Rahman T. Lim S. O'Connell G. D. Holsgrove T. P. ( 2020 ). Influence of testing environment and loading rate on intervertebral disc compressive mechanics: an assessment of repeatability at three different laboratories . JOR Spine 3 ( 3 ), e21110 . 10.1002/jsp2.1110 33015585 PMC7524251

O'Brien L. Fitzmaurice G. ( 2004 ). Analysis of longitudinal multiple source binary data using generalized estimating equations . Appl. Stat. 53 ( 1 ), 177 – 193 . 10.1046/j.0035-9254.2003.05296.x

Oefelein M. G. Conrad V. R. Resnick M. I. ( 2002 ). Skeletal fractures negatively correlate with overall survival in men with prostate cancer . J. Urol. 168 ( 3 ), 1005 – 1007 . 10.1016/s0022-5347(05)64561-2 12187209

Oster G. Lamerato L. Glass A. G. Richert-Boe K. E. Lopez A. Chung K. ( 2013 ). Natural history of skeletal-related events in patients with breast, lung, or prostate cancer and metastases to bone: a 15-year study in two large US health systems . Support Care Cancer 21 ( 12 ), 3279 – 3286 . 10.1007/s00520-013-1887-3 23884473

Pahr D. H. Reisinger A. G. ( 2020 ). A review on recent advances in the constitutive modeling of bone tissue . Curr. Osteoporos. Rep. 18 , 696 – 704 . 10.1007/s11914-020-00631-1 33068252 PMC7732794

Pahr D. H. Schwiedrzik J. Dall'Ara E. Zysset P. K. ( 2014 ). Clinical versus pre-clinical FE models for vertebral body strength predictions . J. Mech. Behav. Biomed. Mater 33 , 76 – 83 . 10.1016/j.jmbbm.2012.11.018 23333770

Palanca M. Cavazzoni G. Dall'Ara E. ( 2023 ). The role of bone metastases on the mechanical competence of human vertebrae . Bone 173 , 116814 . 10.1016/j.bone.2023.116814 37257631

Pond G. R. Sonpavde G. de Wit R. Eisenberger M. A. Tannock I. F. Armstrong A. J. ( 2014 ). The prognostic importance of metastatic site in men with metastatic castration-resistant prostate cancer . Eur. Urol. 65 ( 1 ), 3 – 6 . 10.1016/j.eururo.2013.09.024 24120464

Preve D. Lenarda P. Bianchi D. Gizzi A. ( 2024 ). Phase field modelling and simulation of damage occurring in human vertebra after screws fixation procedure . Comput. Mech. 10.1007/s00466-024-02450-y

Radovitzky R. Seagraves A. Noels L. ( 2011 ). A scalable 3D fracture and fragmentation algorithm based on a hybrid, discontinuous Galerkin, cohesive element method . Comput. Methods Appl. Mech. Eng. 200 ( 1-4 ), 326 – 344 . 10.1016/j.cma.2010.08.014

Saad F. Lipton A. Cook R. Chen Y. M. Smith M. Coleman R. ( 2007 ). Pathologic fractures correlate with reduced survival in patients with malignant bone disease . Cancer 110 ( 8 ), 1860 – 1867 . 10.1002/cncr.22991 17763372

Schwiedrzik J. J. Wolfram U. Zysset P. K. ( 2013 ). A generalized anisotropic quadric yield criterion and its application to bone tissue at multiple length scales . Biomechanics Model. Mechanobiol. 12 , 1155 – 1168 . 10.1007/s10237-013-0472-5 23412886

Siegel R. DeSantis C. Virgo K. Stein K. Mariotto A. Smith T. ( 2012 ). Cancer treatment and survivorship statistics, 2012 . CA Cancer J. Clin. 62 ( 4 ), 220 – 241 . 10.3322/caac.21149 22700443

Stadelmann M. A. Schenk D. E. Maquer G. Lenherr C. Buck F. M. Bosshardt D. D. ( 2020 ). Conventional finite element models estimate the strength of metastatic human vertebrae despite alterations of the bone's tissue and structure . Bone 141 , 115598 . 10.1016/j.bone.2020.115598 32829037 PMC9206866

Taneichi H. Kaneda K. Takeda N. N. Abumi K. Satoh S. ( 1997 ). Risk factors and probability of vertebral body collapse in metastases of the thoracic and lumbar spine . Spine 22 ( 3 ), 239 – 245 . 10.1097/00007632-199702010-00002 9051884

Tschirhart C. E. Finkelstein J. A. Whyne C. M. ( 2006 ). Metastatic burst fracture risk assessment based on complex loading of the thoracic spine . Ann. Biomed. Eng. 34 , 494 – 505 . 10.1007/s10439-005-9063-7 16482411

Tschirhart C. E. Finkelstein J. A. Whyne C. M. ( 2007 ). Biomechanics of vertebral level, geometry, and transcortical tumors in the metastatic spine . J. Biomech. 40 ( 1 ), 46 – 54 . 10.1016/j.jbiomech.2005.11.014 16427058

Tschirhart C. E. Nagpurkar A. Whyne C. M. ( 2004 ). Effects of tumor location, shape and surface serration on burst fracture risk in the metastatic spine . J. Biomechanics 37 ( 5 ), 653 – 660 . 10.1016/j.jbiomech.2003.09.027 15046994

Rebecca L. Siegel M. P. H. Kimberly D. Miller M. P. H. Hannah E. Fuchs B. S. ( 2022 ). Cancer statistics, 2022 . 10.3322/caac.21708

van der Velden J. M. Versteeg A. L. Verkooijen H. M. Fisher C. G. Chow E. Oner F. C. ( 2017 ). Prospective evaluation of the relationship between mechanical stability and response to palliative radiotherapy for symptomatic spinal metastases . Oncologist 22 ( 8 ), 972 – 978 . 10.1634/theoncologist.2016-0356 28469043 PMC5553957

Vashishth D. Koontz J. Qiu S. J. Lundin-Cannon D. Yeni Y. N. Schaffler M. B. ( 2000 ). In vivo diffuse damage in human vertebral trabecular bone . Bone 26 ( 2 ), 147 – 152 . 10.1016/s8756-3282(99)00253-7 10678409

Werner B. Ovesy M. Zysset P. K. ( 2019 ). An explicit micro-FE approach to investigate the post-yield behaviour of trabecular bone under large deformations . Int. J. Numer. Method Biomed. Eng. 35 ( 5 ), e3188 . 10.1002/cnm.3188 30786166

Whyne C. M. Hu S. S. Lotz J. C. ( 2001 ). Parametric finite element analysis of verterbral bodies affected by tumors . J. Biomech. 34 ( 10 ), 1317 – 1324 . 10.1016/s0021-9290(01)00086-0 11522311

Wolfram U. Wilke H. J. Zysset P. K. ( 2011 ). Damage accumulation in vertebral trabecular bone depends on loading mode and direction . J. Biomech. 44 ( 6 ), 1164 – 1169 . 10.1016/j.jbiomech.2011.01.018 21295781

World Health Organization Collaborating Centre for Metabolic Bone Diseases ( 2023 ). FRAX® WHO fracture risk assessment tool .

Wu M. Y. Li C. J. Yiang G. T. Cheng Y. L. Tsai A. P. Hou Y. T. ( 2018 ). Molecular regulation of bone metastasis pathogenesis . Cell. Physiol. Biochem. 46 ( 4 ), 1423 – 1438 . 10.1159/000489184 29689559

Yao A. Sarkiss C. A. Ladner T. R. Jenkins A. L. 3rd ( 2017 ). Contemporary spinal oncology treatment paradigms and outcomes for metastatic tumors to the spine: a systematic review of breast, prostate, renal, and lung metastases . J. Clin. Neurosci. 41 , 11 – 23 . 10.1016/j.jocn.2017.04.004 28462790

Zysset P. K. ( 2003 ). A review of morphology-elasticity relationships in human trabecular bone: theories and experiments . J. Biomech. 36 ( 10 ), 1469 – 1485 . 10.1016/s0021-9290(03)00128-3 14499296

Zysset P. K. Curnier A. ( 1996 ). A 3D damage model for trabecular bone based on fabric tensors . J. Biomech. 29 ( 12 ), 1549 – 1558 . 10.1016/s0021-9290(96)80006-6 8945653

Nomenclature