A Viscoelastic Modeling for Failure Analysis of Human Vertebral Bone Undergoing Multi-Rate Compression

1. Introduction

Vertebral fractures are one of the most common bone injuries and are a major clinical problem, especially in older people and in patients with low bone density [1,2,3]. These fractures can lead to severe pain, changes in the shape of the spine, reduced mobility, and a greater risk of future fractures [2,3]. As a result, predicting the mechanical behavior and failure of vertebrae has become an important subject in spine biomechanics and injury research.

The mechanical behavior of trabecular bone and vertebral structures has been investigated using different experimental, analytical, and computational methods. Multiscale micromechanical models have been proposed to describe the anisotropic mechanical properties of vertebral trabecular bone by Haj-Ali et al. [4]. Green et al. [5] have used experimental studies, together with finite element analysis to examine stress–strain behavior and the start of damage in trabecular bone. In addition, models such as poroelastic formulations have been developed to study the effect of fluid flow inside the porous trabecular structure by Lim and Hong [6].

Finite element models have been widely used to investigate the mechanical behavior of spinal components and vertebrae under various loading conditions. In addition, viscoelastic formulations have been applied to examine load sharing within spinal structures at different loading rates [7] and to identify critical loading conditions through nonlinear analyses [8]. Experimental studies have also characterized the viscoelastic behavior of trabecular bone across multiple length scales using techniques such as dynamic mechanical testing and nanoindentation [9]. Many studies have demonstrated that the mechanical properties of trabecular bone are strongly dependent on apparent density. Several works have reported power-law relationships between elastic modulus and density, although these relationships vary with anatomical location and microstructural characteristics [10,11,12,13]. Wu et al. [14] have also shown that these properties can vary significantly and that density plays an important role in bone stiffness and strength. While density-dependent elastic behavior has been extensively studied, relatively few analytical models incorporate both density dependence and strain-rate sensitivity within a unified framework. In addition to density effects, trabecular bone exhibits pronounced strain-rate dependence. Under dynamic loading conditions such as falls, sports impacts, or motor vehicle accidents, its mechanical response can differ significantly from that observed under quasi-static loading. Previous research has shown that both stiffness and strength increase with increasing strain rate [11,12,13]. This behavior is primarily attributed to the porous microstructure of trabecular bone, which contains marrow and fluid and exhibits time-dependent mechanical effects [15]. Several studies have investigated the viscoelastic properties of trabecular bone using creep and stress relaxation tests. For example, Manda et al. [16] examined the relationship between viscoelastic behavior and bone volume fraction using creep tests on bovine trabecular bone, and later extended this work to include nonlinear viscoelastic effects [17]. However, most of these experiments were performed at low strain rates (about 0.01 s⁻¹), which may not fully represent traumatic conditions where strain rates are much higher. As a result, viscoelastic parameters obtained from low strain-rate testing may not completely describe the behavior of vertebral bone under dynamic loading conditions. Including strain-rate effects by using compression tests at different loading rates may improve prediction accuracy. Although many studies have investigated density-dependent elastic properties, only a limited number of studies have included density-dependent viscoelastic parameters in analytical models of vertebral mechanics. Developing models that consider both density changes and rate-dependent behavior is still an important challenge.

This study presents an analytical model to predict the mechanical response and fracture force of human vertebral specimens under different loading rates. It also investigates the effect of including viscoelastic behavior compared with a purely linear elastic model. The model incorporates density-dependent elastic properties together with a Kelvin–Voigt representation of trabecular bone. Unlike many previous studies in which viscoelastic parameters were derived primarily from low strain-rate creep or relaxation experiments, the present approach identifies these parameters directly from vertebral compression tests performed at different loading rates.

In this context, although finite element methods can provide detailed predictions of vertebral mechanics, they usually require complex image processing, specialized expertise, and high computational cost. Therefore, the analytical model presented in this study can provide a practical and computationally efficient approach for vertebral evaluation.

2. Materials and Methods

2.1. Vertebral Specimens and Experimental Data

Experimental data used in this study were obtained from a previously published biomechanical investigation on cadaveric spinal segments conducted by Rezaei et al. [18] at the Mayo Clinic. In the original study, 28 spinal specimens were tested under compression loading conditions. From these specimens, 10 intact vertebrae were selected for the present study. Three vertebrae were tested under quasi-static loading, and seven were tested at higher strain rates to simulate fracture conditions. The tests gave force–displacement curves, which were used to calculate stiffness and fracture force for each specimen. These results were used as reference data to calibrate the model in this study.

2.2. Image-Based Density and Geometric Measurements

Quantitative computed tomography (QCT) images were used to obtain geometric and density-related parameters for each vertebra. For each specimen, the vertebral body was divided into cortical and trabecular regions using several axial CT slices. The density and cross-sectional area of the cortical and trabecular regions were measured in each slice. These values were averaged along the vertebral height to obtain the mean cortical density, trabecular density, cortical cross-sectional area, and trabecular cross-sectional area for each specimen. The average values were used as input parameters in the analytical model. The specimen specifications are shown in Table 1.

2.3. Model Description

The vertebra was represented using an analytical model in which the vertebral body was approximated as a cylindrical structure composed of two mechanically distinct components: cortical bone and trabecular bone. A schematic illustration of the vertebral body, including the cortical shell and trabecular core, together with the corresponding equivalent mechanical model, is presented in Figure 1. The cortical region was modeled as a linear elastic material, whereas the trabecular region was described using a viscoelastic formulation. Under axial compression, cortical and trabecular bone experience the same overall deformation. Therefore, the two components were assumed to act mechanically in parallel, resulting in identical strains in both trabecular and cortical bone, i.e., ε = ε_t = ε_c. where ε denotes the overall axial strain of the vertebra, and ε_t and ε_c are the trabecular and cortical bone strains, respectively. The viscoelastic behavior of the trabecular bone was represented using a Kelvin–Voigt model consisting of an elastic spring and a viscous dashpot connected in parallel. The Kelvin–Voigt formulation has been widely used to characterize the time-dependent mechanical response of biological tissues [19]. Accordingly, the constitutive relation for trabecular bone can be written as

$${\sigma }_t=E_t+\mathsf{\eta }\dot{\epsilon }$$

(1)

where σ_t is the stress in the trabecular part, E_t is the trabecular elastic modulus, η is the viscous coefficient, and ε denotes the strain rate.

The cortical component was assumed to behave as a purely elastic material, described by σ_c = E_cε Since the cortical and trabecular components act in parallel, the total stress carried by the vertebra can be expressed as the sum of the stresses in each phase,

F = F_c + F_t

(2)

where F_c and F_t denote the forces carried by cortical and trabecular bone, respectively. These forces can be expressed as F_c = A_c σ_c and F_t = A_t σ_t where A_c and A_t represent the cross-sectional areas of the cortical shell and trabecular core. Substituting the constitutive relations into the force equilibrium equation yields

$$F={\mathrm{A}}_c{E}_c\epsilon +{\mathrm{A}}_t(E_t\epsilon +\mathsf{\eta }\dot{\epsilon })$$

(3)

2.4. Density-Dependent Material Formulation

To consider differences in bone density between specimens, density-dependent material relationships were used. The elastic modulus of the bone was defined as a power-law function of density. This type of relationship has been widely reported in experimental and analytical studies on trabecular bone [20,21,22].

E = aρ^b

(4)

where, E is the elastic modulus, ρ is the ash density, a and b are material constants obtained through parameter calibration. In order to account for density-dependent viscoelastic behavior, the Kelvin-Voigt model's viscous coefficient was also defined as a linear function of bone density. The well-established relationship of trabecular bone mechanical characteristics on density shown in the literature [23,24,25,26,27] is consistent with this idea. The relationship is stated as follows,

η = cρ + d

(5)

where c and d are constants determined through the optimization procedure.

This formulation enables the analytical model to capture both the density dependence of bone stiffness and the strain-rate sensitivity associated with the viscoelastic behavior of trabecular bone.

2.5. Fracture Force Estimation

The compressive fracture force of each vertebral specimen was estimated using the effective stiffness predicted by the model. The effective Young’s modulus, which includes both elastic and viscous effects, was then used to estimate the corresponding stress response. The combined contribution of the cortical and trabecular regions was taken into account when calculating the effective modulus of the vertebral column. The effective modulus can be written as follows using the Kelvin–Voigt formulation:

$$E_{\mathit{eff}}=\frac{{\mathrm{A}}_c{E}_c+{\mathrm{A}}_t(E_t+\mathsf{\eta }\frac{\dot{\epsilon }}{\epsilon })}{{\mathrm{A}}_{\mathit{total}}}$$

where the cross-sectional areas of cortical bone, trabecular bone, and their total area are represented by A_c, A_t, and A_total, respectively. Based on this, the effective modulus and the shape of the specimen were used to calculate the stiffness of each specimen. The compressive stress was calculated by assuming linear elastic behavior up to failure σ = E_effε where E_eff represents the effective modulus derived from the analytical model. To estimate the fracture condition, a critical compressive strain of 10% was assumed for vertebral failure [28,29,30]. This assumption has been widely adopted in previous experimental and computational studies of vertebral compression behavior [28,29,30]. The corresponding fracture force was then calculated as F_fracture = σ_fA_Total Where σ_f is the stress at the assumed failure strain.

2.6. Normalization of Specimen Geometry

Because vertebral specimens differ in size and geometry, the force–displacement data were normalized to remove geometric effects. The applied force was divided by the cross-sectional area of the vertebra, and the displacement was divided by the specimen height. This procedure yields a normalized stiffness that reflects intrinsic material behavior rather than geometric differences among specimens.

2.7. Parameter Calibration Identification

The experimental data used to calibrate the model were taken from a published study [18]. These data come from compression tests performed on human lumbar vertebral specimens under controlled loading conditions. The available dataset is somewhat sparse because of the inherent limits of experimental testing of biological tissues, such as restricted specimen availability, heterogeneity in bone quality, and testing complexity. However, the dataset reflects the general pattern of the mechanical response across the explored density range and offers a suitable basis for model calibration despite the small number of data points and inherent variability.

Optimization methods are widely used in structural and biomechanical engineering to identify parameters and improve design. For example, they have been used to find the best placement of piezoelectric patches to reduce stress in smart structures [31,32]. In the present study, an optimization procedure based on the Nelder–Mead algorithm was employed to identify the unknown material parameters of the proposed analytical model.

The unknown material constants a, b, c, and d, which define the density-dependent elastic and viscous relationships, were determined using the experimental data obtained from biomechanical testing. For each specimen, the analytical model was used to estimate the mechanical stiffness based on the specimen-specific geometric parameters and density measurements. The optimal parameters was obtained by minimizing the normalized sum of squared differences between the model predictions and the experimentally measured stiffness values across all specimens. The error function was defined as.

$$\underset{a,b,c,d}{\mathit{min}}{\displaystyle \underset{i=1}{\overset{n}{\Sigma }}}{\left(\frac{k_{\mathit{model},i}-k_{\mathit{exp},i}}{{\overline{k}}_{\mathit{exp},i}}\right)}^2$$

where k_model,i, k_exp,i denote the predicted and experimentally measured stiffness values for specimen i, respectively, ${\overline{k}}_{\mathit{exp},i}$ represents the mean experimental stiffness across all specimens, and n is the total number of specimens. The experimental stiffness was determined by calculating the slope of the force-displacement curve, whereas the projected stiffness was achieved using the effective modulus obtained from the viscoelastic model. The Nelder–Mead optimization algorithm was used to determine the parameter values that minimize this objective function.

2.8. Evaluation of Model Predictions

The predictive ability of the analytical model was evaluated by comparing the model results with the experimental measurements. First, relationships between density and both the elastic modulus and viscous coefficient were obtained using the calibrated parameters. These relationships were used to estimate the mechanical behavior of each specimen. Effective stiffness and compressive fracture force of each vertebra were calculated using the predicted material properties and specimen geometry. The coefficient of determination (R²) was used to evaluate model performance, and scatter plots were used to visually compare the predicted and experimental values.

3. Results

3.1. Calibrated Density-Dependent Material Parameters

The optimization procedure determined the material constants appearing in Eqs. (4) and (5), which describe the density dependence of the elastic modulus and viscous coefficient of trabecular bone. The calibrated values are summarized in Table 2.

The identified parameters indicate that both elastic stiffness and viscous resistance increase with increasing bone density. This trend suggests that denser trabecular bone provides greater load-bearing capacity together with stronger resistance to rate-dependent deformation.

3.2. Stiffness Prediction

The accuracy of the proposed model was evaluated by comparing the predicted normalized stiffness values with those obtained from experimental force–displacement curves [18].

The predicted stiffness was calculated using the effective modulus from the Kelvin–Voigt model, which includes both elastic and viscous effects. The model gives an R² value of 0.45 for all samples (Figure 2), showing moderate agreement with the experimental results. In comparison, the purely linear elastic model gives a slightly lower R² value of about 0.42.

Table 3 summarizes the stiffness and fracture force values from the experiment and the model, together with the accompanying errors.

3.3. Fracture Force Prediction

As shown in Figure 3, the model predicts fracture force and compares it with experimental data. Peak forces were obtained from compression tests reported in a previous study [15], with fracture force defined as the maximum recorded load.

Overall, the comparison between model predictions and experimental results yields a coefficient of determination of R² = 0.47 (Figure 3a), indicating moderate agreement. In comparison, the purely linear elastic model yields a lower R² value of 0.35, highlighting a more pronounced improvement when viscoelastic effects are included. When only dynamically tested specimens are considered, the correlation improves to R² = 0.60 (Figure 3b). This improvement suggests that the viscoelastic formulation more effectively captures vertebral behavior under dynamic loading conditions.

3.4. Model Accuracy and Error Analysis

Prediction errors were examined to identify specimens with unusually large deviations between model predictions and experimental measurements. Grubbs’ test [33] identified specimen 5186L3 as a statistical outlier in both stiffness and fracture force prediction errors. A significance level of α = 0.05 was used. This specimen showed the largest errors in the dataset, with 341% error in stiffness and 219% error in fracture force. After excluding this outlier, the coefficient of determination increased from R² = 0.45 to R² = 0.53 for stiffness prediction and from R² = 0.47 to R² = 0.60 for fracture force prediction. For dynamically tested specimens, fracture force prediction improved more substantially, with R² increasing from 0.60 to 0.88 after exclusion of the outlier specimen.

Specimen 5186T8, obtained from the same donor, also showed a relatively high stiffness prediction error; however, it was not identified as a statistical outlier and was therefore retained in the primary analysis. These findings suggest that donor-specific factors and microstructural characteristics not captured by density alone may influence the mechanical response of the vertebrae.

It is worth noting that the predicted elastic modulus and fracture force (Figure 4) for two specimens from the same donor (5186L3 and 5186T8) showed large positive errors. In both cases, the model clearly overestimated the results. This suggests that something other than density is affecting the mechanical response of these samples. Density is important, but it does not fully describe the trabecular structure. Other factors, such as anisotropy, connectivity, cortical thickness, possible microdamage, and specimen geometry, may also influence the dynamic behavior of the vertebra.

If only density is used, the mechanical properties may be overestimated, especially when the internal structure of the bone is degraded.

To better understand the source of these errors, the relationship between the error in elastic modulus and different density measures (cortical, trabecular, and apparent density) was studied.

No clear relationship was found between the prediction errors and any of the density measures (R² ≈ 0), as shown in Figure 5. The errors are scattered over the whole range of density values. This shows that density alone cannot explain the differences between the model and the experimental results. This also suggests that other factors, such as trabecular structure, anisotropy, and specimen geometry, play an important role in the mechanical behavior of vertebral bone.

4. Discussion

In this study, a density-dependent viscoelastic model was used to describe the mechanical behavior of human vertebral bone. The model relates the elastic response to bone density and also includes time-dependent behavior using a Kelvin–Voigt model. The predicted stiffness values showed moderate agreement with the experimental data, which is reasonable because trabecular bone has high natural variability.

Samples with similar average density may still behave differently because their internal structures are not the same. Factors such as trabecular orientation, connectivity, and local density variations can strongly affect the mechanical response. Some variation in the data may also be related to the assumptions used in the model. In addition, the presence of an outlier specimen suggests that the model is sensitive to specimen-specific features that are not fully represented by density-based relationships. Although the Kelvin–Voigt model can represent rate-dependent behavior, it does not account for nonlinear behavior or damage. These effects may become more significant at higher loads and may explain some of the differences between the predicted and experimental results.

The relatively small improvement of the present model compared with the purely linear elastic model in stiffness prediction may be related to the limited effect of viscosity within the studied strain-rate range. Considering the natural variability of cadaveric bone and the assumptions used in the model, this level of agreement still indicates that the model can capture the general trend of stiffness across the specimens. In particular, the viscous coefficient increased with density, suggesting that denser trabecular bone may resist rate-dependent deformation more effectively. This observation is consistent with previous studies reporting that bone stiffness and strength increase with strain rate.

Additionally, in the present model, this behavior is represented through the viscous component, allowing the density-dependent viscous coefficient to reflect the combined effects of density and strain rate. A comparison with a purely linear elastic formulation showed that the viscoelastic model provided modest improvement in stiffness prediction and more pronounced improvement in fracture force prediction. Model performance also improved for specimens tested under dynamic loading conditions, indicating that the viscoelastic formulation became more important at higher strain rates.

For fracture force prediction, the model relies on an assumed failure strain. Although this is a common method, it introduces uncertainty because the actual failure strain may vary between specimens. Overall, the results suggest that density alone cannot fully explain variations in mechanical response. Microstructural characteristics and specimen-specific differences also appear to play an important role. Nevertheless, the model provides a useful representation of vertebral behavior under loading conditions related to spinal injury.

Compared with patient-specific finite element approaches, the present analytical framework has lower computational cost, requires fewer CT-based inputs, and avoids complex meshing and solver procedures. These features may make the model more suitable for rapid evaluation and wider clinical use.

Several limitations should be noted. First, the sample size is limited, which may contribute to the moderate agreement between the model predictions and the experimental data. Increasing the number of specimens could improve the results. Second, the vertebra was modeled as a cylinder, while real vertebrae have more complex and nonuniform geometries that can influence their mechanical behavior. Finally, the Kelvin–Voigt model does not account for nonlinear behavior or damage mechanisms, which may become important at higher load levels.

5. Conclusions

In this study, a density-dependent viscoelastic model was used to describe the mechanical behavior of human vertebral bone. The model combines a power-law relationship between elastic modulus and density with a density-dependent viscous term in a Kelvin–Voigt model. The results showed moderate agreement with the experimental data and improved prediction of mechanical response, especially under dynamic loading conditions. By including viscoelastic behavior, the model could represent the rate-dependent behavior of vertebral bone. The findings showed that although bone density is an important factor in the mechanical behavior of vertebral bone, it is not sufficient to fully explain the observed behavior. Other factors, such as bone microstructure, may also affect the behavior. Including these factors in the model may improve prediction accuracy. Overall, the model provides a useful approach for describing vertebral bone behavior under loading conditions related to spinal injury. The proposed model may also be a practical and computationally efficient alternative to more complex finite element models for preliminary assessment of vertebral strength.