An Experimentally Validated Mathematical Model of Axial Spinal Compression to Visualize Vertebral Compression Fracture

Summary

Here, we provide a predictive model of VCFs to support surgical intervention. Biomechanical fracture patterns of human and porcine cadaveric vertebrae under compression loading were used to verify an FEA model, providing a resource for better visualizing vertebral fracture and supporting using porcine samples in spine research.

Abstract

The finite element analysis (FEA) model established in this study predicts the mechanical behavior of a vertebral body under pure compression loading. Four human and four porcine cadaveric spines were used in the region from T9-T12 (human) and T12-T14 (porcine) to derive biomechanical failure data under pure axial compression. By implementing the data from these axial crush experiments and combining them with computed tomography (CT)-derived three-dimensional (3D) reconstructions of vertebrae, this new model mathematically predicts the behavior of the spine in vertebral compression fractures (VCFs). The development of an accurate mathematical model has the potential to aid surgical intervention of VCFs during vertebroplasty, kyphoplasty, and potentially additional techniques by accurately representing biomechanical behavior. It can be adapted for osteoporosis, fusion, and implant biomechanics cases to improve surgical accuracy. Understanding the vertebral fracture pattern is essential in surgical pre-planning and understanding prior and future risks. The presented method may serve as an additional resource for studying and treating VCF and further supports the use of porcine specimens in spine research.

Introduction

Vertebral Compression Fractures (VCFs) are biomechanical failures of the anterior vertebral body in response to axial compressive load^1,2. There are 1-1.5 million patients experiencing VCFs in the US each year, with notable risk factors including osteoporosis, advanced age, and female sex^3,4,5. Patients with VCFs can experience pain, disability, altered pulmonary or respiratory function, secondary vertebral fracture, and increased mortality risk^6,7. Preferred clinical management includes physical therapy, rehabilitation, and pain management¹. Surgical intervention for VCF treatment includes vertebroplasty and kyphoplasty, the Osseo-Fix Spinal Fracture Reduction System, internal bracing, and anterior and posterior decompression and stabilization^2,8,9. Clinical diagnosis of VCF may involve magnetic resonance imaging (MRI), CT, positron emission tomography (PET), and single-photon emission computed tomography (SPECT)¹⁰-however, additional visualization and simulation of vertebral fracture may circumvent limitations of these imaging modalities. This paper aims to outline a finite element method of VCF, validated by physical compression experiments, that allows for accurate modeling and visualization of fracture biomechanics.

Finite element analysis (FEA) is a widely used engineering method to simulate and mathematically predict an object's response to physical conditions¹¹. While numerous studies of VCF utilize FEA^12,13,14, many of these studies include CT-derived modeling from living patients and would benefit from a side-by-side validation of the vertebra's response in actual physical compression fracture. This study outlines not only the FEA process but also the protocol for the physical simulation of axial compression fracture.

FEA has been used in surgical pre-planning in a variety of clinical settings, including colorectal repair, oral and maxillofacial procedures, and orthopedic reconstruction^15,16,17. FEA provides a virtual 3D recreation of a site or organ, which may facilitate anatomic understanding of a unique patient's surgical needs. A finite element model can recreate the size, shape, orientation, and behavior of a physical object in response to loads or external pressures, which is useful in optimizing surgical techniques.

This study aims to describe this additional model to help optimize the visualization of VCF, which may augment treatment by vertebroplasty and balloon kyphoplasty. The method described may be beneficial in standardizing the research performed in biomechanical fracture experiments and offers advantages over traditional protocols related to cost, time, and sample selection. This method also provides support for previously constructed mathematical models of compression fractures that utilized medical imaging and finite element analysis^12,13,14,18. Supplemental visualization of the VCF morphology can provide a tool to optimize treatment. This study proposes a mathematical model that predicts VCF under axial compression loading, supported by physical experimental testing, in human and porcine specimens.

Protocol

All experiments were performed in accordance with relevant human research ethics and animal care guidelines at the University of Illinois College of Medicine.

1. Sample preparation

NOTE: VCFs are extensively described for the thoracolumbar region, as most wedge compression fractures occur in the mid-thoracic area¹⁹. For this reason, thoracic vertebrae T9-T12 from humans and T12-14 from porcine were used for these experiments (Table 1).

1. Obtain four cadaveric human spines.
2. Obtain four commercially available cadaveric porcine spines.
3. Identify vertebrae T9-12 and T12-14 for the human and porcine spines, respectively. Ensure there are a total of sixteen human vertebrae and twelve porcine vertebrae.
4. Using a round-tip bone saw and dissection scalpel, dissect the identified vertebrae and remove soft tissue, muscle ligaments, and intervertebral discs.
5. Perform CT scans for each spinal vertebra with the specifications listed in Table 2.
6. Using dissection tools, remove posterior elements of the vertebrae (spinous process, transverse process, facet joints, lamina, and pedicles), following the curvature between the pedicle and vertebral body, such that only the body remains (Supplementary Figure 1).
7. Polish the superior and inferior surfaces of each vertebral body to ensure a planar surface for compression. Remove soft tissue debris and remaining tissue on the vertebral superior and inferior surfaces with a standard scalpel and forceps.
8. Obtain 15 and 12 cortical cores from the upper and lower facets of the vertebrae for human and porcine specimens, respectively. Obtain cortical cores using an oscillating bone saw to remove a 1 cm x 1 cm sample of cortical bone and use a standard scalpel to remove the remaining cancellous bone.
9. Perform density test on samples by weighing each cortical core within a volume-measured pipe to determine scaling (Supplementary Figure 2).
10. In preparation for compression experiments, use the CT images to create a 3D-printed plastic polymer guide of the vertebral body lower profile (Supplementary Figure 3).
11. Using a laser cutter machine, create four holes to define sagittal and coronal planes through the center of mass.

2. Creating a 3D model mesh from radiographic imaging (Supplementary Figure 4)

1. Open the MIMICS (image processing software) software.
2. Import the CT images previously obtained in step 1.5. Click File > New Project Wizard > Import DICOM File.
3. Create a new Threshold by selecting New under the Masks tab and apply the default settings selected for Bone (CT).
4. Rename the new mask appropriately (e.g., Human Spine 1, T9).
5. Click Control + Select the new mask, select Calculate 3D > Optimal Quality > Calculate. The 3D model will appear in the lower right quadrant.
6. Select File > Export > 3-matic > Select 3D Object, which contains the named file.
7. Select Fix > Smooth > Select 3D object > Apply at Smooth Factor 0.5.
8. Select Fix > Fix Wizard > Follow Advice. Continue to click Follow Advice until the # detected prompt appears green. This step optimizes edges, shells, and overlapping triangles.
9. Select Design > Create Analytical Primitive > Create Datum Plane.
10. Select points and Apply for Points 1-3 of the Datum Plane. The average cranial and caudal surfaces defined the transverse plane parallel with the superior and inferior endplates. The sagittal plane is defined as perpendicular to the transverse plane, passing through the center of mass from the anterior to posterior direction of the vertebrae following the spinous process.
11. In the lefthand toolbar, change the surface mesh from Smooth Shaded to Filled with Triangle Edges.
12. Select Remesh > Uniform Remesh > Set Triangle Edge Length to 1.0 mm > Apply.
13. Select Remesh > Create Volume Mesh > Apply.
14. Save *.INP file.

3. Assigning material properties to the volume mesh (Supplementary Figure 5)

1. Export the file back into the image processing software.
2. Convert CT grayscale values to real density values using mineral calibration as previously outlined, using density values obtained in step 1.9^20,21.
   NOTE: The resulting values for the Young modulus (E), the yield strain (εy), and the Poisson Coefficient (ν) are shown below. The Young modulus is calculated as a function of the apparent density and assumed ratio ρ_app/ρ_ash equal to 0.6.
   
3. Under FEA, select Material Assignment.
4. Assign values for Young's modulus, yield strain, and Poisson Coefficient as determined above. The minimum grayscale value limits the density to 0.01 g/cm³ to avoid unobtainable negative density values.
5. Click Export > Abaqus > Select File> Add > Check boxes for Single Output File and Create Assembly > Select OK.
6. Open Abaqus CAE > Right click > Properties > Start folder.
7. Right click > New > Rich Text Format > Rename and convert from .txt to .m
8. Open the new file in MATLAB > in Editor, assign post-yield behavior as perfectly plastic, as described by the following equation, with constraint (σ):
   
9. Click Run to perform the function.

4. Finite element analysis

1. Select Module > Step > Create > Dynamic, Explicit > Continue > Ok.
2. Select Module > Load > Boundary Condition Manager > Step: Initial > Displacement/Rotation > Continue.
3. Select the inferior boundary of the mesh and select Done.
4. Select Module > Mesh > Select Part > Assign Element Type > Done > 3D Stress > Select Yes for Element deletion.
5. Select Module > Load > Boundary Condition Manager > Create > Step-1 > ensure category is Mechanical > Continue.
6. Select the superior boundary of the mesh and select Done.
7. Select Module > Interaction > Reference Point > Assign to Mesh.
8. Select Constraint Manager > Create > Rigid Body > Tie and select superior boundary of mesh > Assign Reference Point > Done. Find the controlled point location using the follower load technique, where the follower load is tangent to the curve of the spine, ensuring that the controlled point is located at 10% of the width of the vertebral body along the sagittal plane from the center of mass projection from the upper plate (Supplementary Figure 6)²².
9. Select Module > Load > Boundary Condition Manager > Create > assign load at the reference point.
10. Within Edit Boundary Condition, select Create Amplitude > Continue > assign displacement rate of 3 mm/min.
11. Select Module > Property > Material Manager > Create > Mechanical > Elasticity > Elastic > assign input data (Table 3).
12. Select Mechanical > Plasticity > Plastic.
13. Select Module > Step > Field Output Requests Manager > Created > Failure/Fracture > select Compressive Damage and State/Field/User/Time > STATUS.
14. Select Module > Job > Job Manager > Results > Visualization > DAMAGE > Frame Selector > increase load using Step Slider.
15. Document the point of failure previously defined.

5. Experimental testing

1. Configure the upper loading plate to the MTS connector, allowing free rotation of up to 5°.
2. Place the vertebral body on the lower compression plate.
3. Trace the silhouette of the vertebral body to ensure that the center of loading is aligned with 10% of the width of the vertebral body from the center of mass in the sagittal direction.
4. Place the vertebral body lower profile onto the outline, oriented by the 3D-printed polymer piece described in steps 1.10 and 1.11.
5. Orient data tracking equipment to capture motion changes and fix the digital marker on the loading plate.
6. Obtain pre-compression photos of the vertebral bodies in all views. Orient the video capture device.
7. Using the materials test systems (MTS) machine, apply an axial load of a quasistatic displacement rate of 3 mm/min, where failure is set as ⅓ of the vertebra height. The MTS machine automatically outputs the load (N) for each experiment.
8. Use the coordinates of three digitized features on the scale to calculate the displacement from the center marker.
9. Use the coordinates of four digitized points of the sagittal and coronal planes to calculate the upper surface (Supplementary Figure 7).
10. Obtain post-compression photos of the vertebral bodies in all views.

6. Data analysis

1. The MTS collects the load and extension values to compute the load-displacement curve. From both FEA and experimental compression, utilize stiffness and strength values to perform statistical analysis.
   NOTE: Stiffness is the slope in its linear trajectory, and strength is the highest point of elastic response. Finally, the failure pattern is defined by identifying any fracture along the compression load from finite elements with a nonzero plastic strain.
2. Perform all statistical analysis using appropriate data analysis software.

Results

From the 12 porcine vertebrae, the load-displacement curves obtained for all specimens are depicted in Figure 1, where all experimental data curves showed similar patterns. Results show an average vertebral stiffness and strength of 9.54 ±1.1 kN/mm (range 6.73-13.63 kN/mm) and 10.2 ± 0.86 kN (range 8.1-13 kN). Video analysis of the experimental compression tests demonstrated a common failure pattern across specimens. The failure pattern consisted of three fracture lines. The first two fractu...

Discussion

While CT-derived FEA is a well-established tool for modeling biomechanical behavior, few studies apply mathematical modeling in conjunction with physical validation in the context of VCFs. This study demonstrates an accurate predictive mathematical model of thoracolumbar VCF, supported by experimental compression of porcine and human cadaveric spines. The model incorporates strength, stiffness, fracture pattern, load mechanics, and material properties to mathematically predict vertebra response to axial load. The morphol...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
3-MATIC	Materialise	v15.0
Abaqus/CAE	Simulia	n/a
Analytical Balance	Sigma Aldrich	OH30122634
Cadaveric porcine spines			Local butcher shop, Chicago, IL
CT machine	Siemens	SOMATOM Definition AS
Epredia Shandon Round Tip Bone Saw, Blade: 8 in. x 2.5 in. (20.3 in. 6.4 cm), standard, 11.5 in. (29.2 cm)	Fischer Scientific	40018
Fisherbrand Fine Precision Medium Tipped Tweezers/Forceps	Fischer Scientific	12-000-157
Integra Miltex Sterile Standard Scalpels	Fischer Scientific	12-460-451
Laser cutter machine/3D printer	Fusion3	F400HFR
Mimics Software	Materialise	v2.0
NDI Optotrak Certus	NextGen Ergonomics	Certus Products
Orthopedic Bone Saw	OrthopedicDrills	OTS-1
SPSS Statistics	IBM	v27
Absorbent Underpads with Waterproof Moisture Barrier, 58.4 x 61 cm, 410 mL, Blue, 200/cs	VWR INTERNATIONAL INC	56616-032