Through the hospital’s electronic medical record system, one OVCF patient with an L1 vertebral injury and a healthy volunteer without a history of spinal disease were selected as research participants. All participants underwent computed tomography (CT) scans of their entire lumbar spine for diagnostic or differential diagnostic purposes. The OVCF patient was 69 years old, with a height of 155 cm, a body mass of 55 kg, and a body mass index (BMI) of 22.89 kg/m². The patient had no history of medical or surgical disease, physical disability, or trauma. A T-score of −3.5 and T < −2.5 meet the World Health Organization (WHO) diagnostic criteria for osteoporosis. The healthy volunteer was 15 years old, with a height of 180 cm and a body mass of 100 kg. This study was approved by the Ethics Committee of Harrison International Peace Hospital (Approval No. 2024246–1), and written informed consent was obtained from all participants.

In this study, biomechanical mechanisms were investigated in depth through digital computer-aided modeling. Bone mineral density was measured and analyzed using a dual-energy X-ray absorptiometry (DXA) scanner (OsteoSys EXA-3000, South Korea, OsteoSys Co., Ltd.) at our hospital. A 64-slice high-resolution spiral CT scanner (Siemens Biograph mCT PET/CT, Germany, Siemens AG) was used to scan all segments of the thoracolumbar region, with a slice thickness of 0.625 mm, a tube voltage of 140 kV, and a current of 200 mA. A total of 507 images in sagittal, coronal, and axial planes were obtained. The CT data were converted to DICOM format and stored on a CD-ROM. The acquired two-dimensional CT images were then imported into Mimics 21.0 modeling software for the extraction and transformation of the original model contours. Geomagic 2017 software was used for feature modification, surface smoothing, surface fitting, and separation of cortical and cancellous bones. SolidWorks 2024 was utilized for structural generation and assembly, while Ansys Workbench 2022 R¹ was employed for assigning material properties, mechanical meshing, ligament construction, and biomechanical analysis.

The DICOM-format CT data from the healthy volunteer were extracted and imported into Mimics 21.0 software to establish the finite-element model of the normal T12–L2 thoracolumbar spine. The “Bone Window Rendering Thresholding” function was applied to segment the imported files at the bone tissue–soft tissue interface, with a threshold range set from 259 to 1906 HU. The thoracolumbar spine model was reconstructed layer by layer using the “Edit Masks” tool, allowing for precise segmentation of the vertebrae and other structural units. The final model was saved in STL format. Structural unit optimization—including smoothing, grinding, denoising, surface reconstruction, and solidification—was performed using Geomagic Warp 2017 software. Simulated intervertebral disks, endplates, and cartilage were generated in SolidWorks 2024, with the intervertebral disks further subdivided into the nucleus pulposus and annulus fibrosus.

The parameterization of the thoracolumbar segment model was done with the help of Ansys 2022 software. The material properties of the osteoporotic vertebrae were modeled by assuming that all structures were made of homogeneous linear elastic materials. The assignment of material properties to parts of the 3D model was based on those used in recent OVCF studies (Table 1) (Huang et al., 2023; Li et al., 2024; Pan et al., 2024; Yu et al., 2024) (Table 2) (Constant and Murley, 1987). Goel et al. (1993), Zhang et al. (2010), Huang et al. (2023) proposed that the modulus of elasticity of cortical bone, vertebral endplates, and posterior vertebral structures of the thoracolumbar segment of the spine was reduced by 33%, and the modulus of elasticity of the cancellous bone was reduced by 66%, and osteoporosis was simulated by decreasing the modulus of elasticity of each vertebral type by a certain amount. The model is meshed by controlling the mesh type and size to ensure that the computational accuracy meets the analysis requirements. The mesh size of the cartilage is set to 0.5 mm, and the mesh size of the rest of the model is set to 2 mm. The ligament portion is replaced with a spring set to stretch only (Che et al., 2022). The contact of articular cartilage with other surfaces was defined as frictional, with a friction coefficient of 0.1, and the contact of all other surfaces was set as bound (Figure 1).

V-type Fracture Model: The V-type fracture model was constructed in Geomagic 2017 based on vertebral data from the healthy volunteer. Beginning at the posterior vertebral edge vertex, a diagonal penetrating line was created using through-mode line selection, extending anteriorly until 27.5% compression was achieved. The model was then divided into resected and retained portions. After removal of the resected section, the superior surface of the vertebra was reconstructed using the fill command. Feature removal and spike elimination functions were subsequently applied to prevent stress concentration, resulting in the final simulated V-type fracture model.

Double-V-type Fracture Model: The double-V-type fracture model was constructed in Geomagic 2017 using vertebrae from the healthy volunteer. Through-mode line selection was applied at identical angles to both the superior and inferior surfaces from the posterior edge vertex, extending anteriorly until 26.6% compression was achieved. Following segmentation and removal of the resected portions, the superior surface was reconstructed using the fill command. Feature removal and spike elimination were subsequently performed to minimize stress concentration, resulting in a simulated double-V-type fracture model.

U-type Fracture Model: The U-type fracture model was constructed in Geomagic 2017 using vertebrae from the healthy volunteer. The brush tool was uniformly applied across the superior and inferior surfaces, with iterative offset adjustments used to create an inward depression until 27.4% compression was achieved. A sculpting knife was then used to blunt sharp edges and prevent stress increasers, followed by the application of quick smoothing and sandpaper functions to remove excess material, resulting in the finalized U-type fracture model.

Extracted Fracture Model: A biconcave fractured vertebra, confirmed through clinical assessment, was processed using Geomagic 2017. Sandpaper and feature removal tools were applied to achieve comprehensive surface blunting while preserving the native anatomical morphology. The model was then adjusted to 26.5% compression to generate the extracted fracture model (Figure 2).

Referring to the EVOSG typing system (Ismail et al., 1999), this patient’s fracture type was assessed morphologically as biconcave based on basic radiographic measurements, including the degree of vertebral wedging (localized kyphosis) and vertebral height loss (expressed in millimeters or as a percentage), in compression fractures (Ruiz Santiago et al., 2016). Referring to the Genant semi-quantitative (SQ) typing system (Genant et al., 1993), the degree of compression of the fractured vertebral body of the OVCF model was 26.5% and the degrees of compression of the other three models were adjusted to be close to that of the OVCF model. Consequently, a comparative analysis was conducted between moderate biconcave vertebral fractures and other morphologically distinct yet compression-matched moderate vertebral fractures. By maintaining identical compression ratios with varying fracture morphology, the variables were effectively isolated to ensure more targeted and reliable research outcomes. This is consistent with the fact that a reduction in vertebral height of >15% is the morphometric criterion required for an imaging diagnosis of a new vertebral fracture (Ruiz Santiago et al., 2016). The degree of vertebral body compression is obtained by taking the posterior margin of the vertebral body as the pre-compression and the bottom end of the compression as the post-compression. That is, the degree of vertebral body compression = [(height of the posterior margin of the vertebral body - height of the lowest end of the vertebral body) ÷ height of the posterior margin of the vertebral body] × 100%.

Future computational models should incorporate paraspinal musculature and intra-abdominal pressure to enhance the clinical relevance of vertebral fracture simulations, given their influence. In the current model, all nodes on the inferior surface of the L2 vertebra were constrained in the X, Y, and Z directions. The coordinate system of the thoracolumbar spine segment was used as a reference frame to assign mechanical parameters corresponding to different degrees of freedom. Loading conditions were then defined based on the X-, Y-, and Z-axes to simulate various motion states. A compressive load of 500 N was applied coaxially along the Z-axis to the coupling node on the superior surface of the T12 vertebral body, simulating the vertical loading experienced by the thoracolumbar spine in an upright posture (Pan et al., 2024). This load was uniformly distributed across the entire surface. To simulate flexion and extension, torques of +10 N/m and −10 N/m were applied about the Y-axis. Side-bending motions were modeled by applying torques of +10 N/m and −10 N/m about the X-axis. For rotational loading, torques of +10 N/m (left rotation) and −10 N/m (right rotation) were applied about the Z-axis. These boundary conditions were used to solve for the equivalent forces and displacements of the vertebral bodies under different physiological motion scenarios.

All data were analyzed using SPSS Statistics version 27.0. For continuous variables, normality and homogeneity of variance were first assessed based on the maximum displacement values across the four fracture models. When both assumptions were met, one-way analysis of variance (ANOVA) was performed, followed by post hoc comparisons using the Student–Newman–Keuls q (SNK-q) test. For data that did not satisfy normality assumptions, the nonparametric Wilcoxon rank-sum test (Mann–Whitney U test) was applied. Results are reported as the mean ± standard deviation, with statistical significance set at P < 0.05.

The current finite-element model was validated under combined loading conditions of pure moments and follower loads. Normal bone density was first assigned to the spinal materials, and the model’s mobility was evaluated in all six anatomical directions. No significant differences were observed when compared with the experimental results reported in previous studies (Schmoelz et al., 2010; Huang et al., 2023). Subsequently, osteoporotic bone density was assigned to the model for force analysis.

Comparison of T12 vertebral force clouds revealed significant variations in force distribution under different loading conditions (Figure 3). According to the stress distribution of the four fracture models in six directions, the maximum stress was primarily concentrated in the cortical bone, and from the stress distribution on the upper surface of the cortical bone, it could be seen that, in the forward flexion direction, the force of the four models was concentrated on the anterior edge of the vertebral body and then distributed along the posterior edge of the vertebral body, which was more obvious in the OVCF and U-type models; in the direction of the posterior extension, the forces of all four models were concentrated on the posterior margin of the vertebral body, with the OVCF model being the most obvious; in the direction of left lateral bending, the forces of all four models were concentrated on the left side of the vertebral body, with the OVCF and U-type models being the most obvious; in the direction of the right lateral curvature, the forces of all four models were concentrated on the right side of the vertebral body, and there was no significant difference among the four models; in the direction of left rotation, the OVCF and U-type models were distributed on the left side of the vertebral body, while the V-type and double-V-type models were distributed in the region where the vertebral body was connected to the pedicle; in the rightward direction, the force on the U-type vertebra was concentrated on the right edge of the vertebrae, and the distribution of the force on the other three types of vertebrae was not obvious.

Comparison of the four intervertebral disk force cloud models revealed significant variations in force distribution among the models under different loading conditions (Figures 4, 5). In the disks of T12 and L1, there was no significant difference in the distribution of force among the four models. In the disks of L1 and L2, the U-type model was the closest to the OVCF model in the direction of forward flexion. There was no significant difference in the patterns of stress distribution of the four models in any other direction.

Comparison of L1 vertebral body stress clouds revealed significant differences in stress distribution under varying loading conditions (Figure 6). According to the stress distribution of the four fracture models in the six directions, the maximum stress was concentrated in the cortical bone, and the stress distribution on the upper surface of the cortical bone could show that, in the forward-flexion direction, the U-type model was closer to the OVCF model, which were all concentrated in the anterior margin of the vertebral body, and there was no specificity in the stress distribution of the V-type and double-V-type models; in the left-rotation direction, the double-V-type model was closer to the OVCF model, and it was distributed in the left rotation direction; the double-V-type model was closer to the OVCF model and was distributed on the right side of the posterior vertebral body edge; in the right rotation direction, the U-type model was closer to the OVCF model and was distributed on the left side of the posterior vertebral body edge; in the remaining directions, the four models were in good agreement.

Comparison of L2 vertebral body stress clouds revealed significant differences in stress distribution under different loading conditions (Figure 7). According to the stress distribution of the four fracture models in six directions, the maximum stress was concentrated in the cortical bone, and from the stress distribution on the upper surface of the cortical bone, it could be seen that, in the direction of posterior extension, the V-type model was closer to the OVCF model, with a significant part concentrated on the posterior edge of the vertebral body and a small part in the anterior edge of the vertebral body; in the direction of the left lateral curvature, the U-type was closer to the OVCF model, and was concentrated in the posterior edge of the left side of the vertebral body; the remaining directions showed no significant difference.

A comparative analysis of maximum vertebral displacements across six loading directions (Table 3; Figures 8–10) revealed that the U-type fracture model exhibited displacements of 1.7876 mm (axial), 6.1564 mm (flexion), 4.6520 mm (left lateral bending), 6.2224 mm (right lateral bending), 3.4119 mm (left rotation), and 3.1601 mm (right rotation) at the L1 level. Morphological comparisons demonstrated the closest displacement matching between the U-type and the extracted models in left lateral bending, while the double-V-type and the extracted models showed maximal similarity during extension and right rotation. A statistical analysis confirmed no significant intergroup differences (p > 0.05). A multi-vertebral (T12-L2) assessment further identified the closest biomechanical agreement between the U-type and the extracted models in flexion, left lateral bending, and left rotation. Representative displacement patterns of the U-type model under multidirectional loading are illustrated in Figure 11.