Optimizing Cement Distribution in the Fractured Area for Osteoporotic Vertebral Compression Fractures Through Biomechanical Effects Analysis Based on a Three Dimentional Lumbar Finite Element Modeling

Lin-qiang Ye Guangzhou University of Traditional Chinese Medicine First A liated Hospital De Liang Guangzhou University of Traditional Chinese Medicine First A liated Hospital Zhen Li Dongguan Hospital of Traditional Chinese Medicine Rui Weng Guangzhou University of Chinese Medicine Xue-cheng Huang Southern Medical University Yong-hong Feng Dongguan Hospital of Traditional Chinese Medicine Jian Zou Dongguan Hospital of Traditional Chinese Medicine Xiao-bing Jiang Guangzhou University of Traditional Chinese Medicine First A liated Hospital Guo-liang Lu (  dglugl@163.com ) Department of Spinal Surgery, Dongguan Hospital of Traditional Chinese Medicine, Dongguan, Guangdong, People’s Republic of China;


Introduction
Osteoporotic vertebral compression fractures (OVCFs) are very common in the elderly, with an estimated 1.4 million new fractures occurring every year worldwide [1]. Until recently, symptomatic OVCFs were treated commonly with conservative methods including bed rest, analgesics, braces, and physical therapy. However, percutaneous vertebral augmentation (PVA), such as percutaneous vertebraplasty and percutaneous kyphoplasty, has been introduced as an alternative treatment option [2,3]. Biomechanical studies showed signi cant increases in the stiffness and strength of individual fractured vertebra after PVA [4,5]. Apart from rapid pain relief, another immediate effect of PVA was an increase of anterior vertebral height, which reduced kyphosis in patients [6][7][8]. The realigned spinal column and regained height in the augmented vertebra may decrease pulmonary and gastro-intestinal complications and early morbidity related to compression fractures [9].
Our previous biomechanical study demonstrated that various cement distributions relative to the fractured area could produce different impacts on biomechanical behaviors of OVCFs treated with PVA [10]. Biomechanical behaviors of su cient versus insu cient cement distributions in the fractured area and symmetrical versus asymmetrical cement distributions around the fractured area were compared.
Insu cient cement distribution in the fractured area led to signi cant increase in maximum displacement of augmented vertebra, which might explain unrelieved pain after PVA in treatment of symptomatic OVCFs. While insu cient cement distribution in the fractured area increased maximum von Mises stress in cortical and cancellous bone of augmented vertebra signi cantly, asymmetrical cement distribution around the fractured area increased maximum von Mises stress in cancellous bone signi cantly, both of which might be the reason for recollapse of augmented vertebra. These ndings were veri ed by related clinical studies [11,12].
While cement distributes su ciently in the fractured area and relatively symmetrically around the fractured area, three types of cement mass location in the vertebral body are commonly seen according to our clinical practice when performing bipedicular PVA for OVCFs, including anterolateral (AL), anteromedial (AM) and posterolateral (PL). To the best of our knowledge, however, few biomechanical studies exist that show the differences of biomechanical behaviors among these three types of cement distribution so far. In order to further optimize cement distribution in the fractured area for OVCFs treated with PVA, this study was performed to investigate biomechanical effects of AL, AM and PL in the fractured area on OVCF on the basis of our previous study [10].

Development of normal L1-L5 nite element (FE) model
A normal three-dimensional digital anatomical FE model of L1-L5 was created based on CT images of a L1-L5 motion segment without injury or radiographic evidence of degeneration. The image data were from a 26-year-old male who requested CT scan for health examination at our hospital. He was explained about the research purpose and signed the consent form. Slice thickness was 0.625 mm. The slice images were preserved in a computer and then imported to Mimics 19.0 software (Materialise Inc., Leuven, Belgium) for generation of the 3-dimensional geometries of L1-L5 vertebra. Then, a smoothing process was performed in Geomagic Studio 2013 software (Geomagic Inc., Research Triangle Park, NC, USA) to remove spikes and holes on the surface of the vertebral geometries. The geometries of intervertebral discs which were di cult to separate from the CT images, were created using SolidWorks 2017 sofeware (SolidWorks Inp., Concord, MA, USA). Each vertebra consisted of a cancellous core surrounded by a cortical shell layer with thickness of 0.5 mm [13]. Cartilaginous endplates were simulated with thickness of 0.5 mm at both ends of each vertebra [14]. The nucleus pulpous occupied about 43% of the total disc [15]. Facet joints were modeled with a cartilage layer with thickness of 0.2 mm and a gap of about 0.6 mm [16]. Element types of cancellous bone, posterior bony element, cartilaginous endplates, annulus and nucleus pulpous were de ned as solid elements. Cortical bone and facet joint cartilage were de ned as shell elements. In order to mimic collagen bers of annulus brosus, four layers of collagen bers were embeded radially in the annular ground substances. For each layer two bundles of bers were modeled with orientation of about ± 30 ° with respect to the horizontal plane using truss elements [17]. A ber cross-sectional area of 0.1mm 2 was assumed [18]. Ligaments were modeled with truss elements. The assigned material properties were assumed to be homogeneous, linear, and isotropic and truss elements were tension-only. Tied constraints were used to ensure disc and ligament attachments to the vertebra and prevent any relative movement during the simulations. Frictionless Surface-based, nite sliding contact was de ned for facet joint cartilages [19]. Material properties were taken from the literatures and listed in Table 1 [18,20,21]. Abaqus 6.12 (Abaqus, Inc., Providence, RI, USA) was used for numerical analysis. A mesh convergence study was conducted and considered to be convergent when the prediction results obtained by two consecutive mesh resolutions had differences smaller than 5% [22,23]. The nal mesh consisted of 8042 solid elements, 177037 shell elements, 3840 truss elements and a total of 285469 nodes (Fig. 1).

Simulation Of Ovcf
Referring to simulation method reported by Chiang et al. [26], the model (Pre-augmented) was constructed with the following steps to simulate OVCF on L3. A cleft was horizontally penetrated into the vertebral body by 20 mm through the center of the anterior cortical shell. The size of the cleft was approximately 20 mm, 30 mm and 2 mm in depth, width and height, respectively (Fig. 2).

Simulation Of Cement Augmentation
For bone cement material, PMMA was applied because it represents the most common clinically used material. Its material property was listed in Table 1 [21]. Two cement cylinders with the same volume were vertically implanted around the fractured area symmetrically to mimic bipedicular PVA with su cient cement distribution in the fractured area and symmetrical cement distribution around the fractured area. The volume of each cement cylinder was approximately 3 ml. For simulation of AL, both of cement cylinders were embedded in the anterior column and in proximity to the lateral border of the vertebral body. Subsequently, the two cement cylinders were relocated medially adjacent to midline of the vertebral body to mimic AM. Being different from AM, the two cement cylinders of AL were relocated posteriorly with partial cement cylinders locating in middle column to simulate PL. Eventually, we got four different models for test, including the fractured model before PVA and three augmented models (Fig. 2).

Loads And Boundary Conditions
A pure moment of 7.5 N m was applied independently for exion/extension, left/right lateral bending, and left/right axial rotation. The inferior endplate of L5 vertebra was xed in all degrees of freedom and Loads were applied in L1 superior endplate. Abaqus 6.12 (Abaqus, Inc., Providence, RI, USA) was used for numerical analysis. The magnitudes and distributions of von Mises stress in cortical and cancellous bone and maximum displacement of L3 vertebral body were recorded. Von Mises stress has been proposed as a parameter of failure criteria for bone [21,27] and maximum displacement a parameter of stability [28].

Results
Validation of the normal L1-L5 FE model Under 7.5 N m pure moments, the ROMs of exion, extension, lateral bendings and axial rotations in each spinal level predicted by the FE model were within the ranges reported by Panjiabi et al. (Fig. 3a-d) [24]. L4-L5 IDPs under compressive forces of 300N and 1000N were in agreement with the results of the in vitro study (Fig. 3e) [25].
Distributions and magnitudes of von Mises stress in cortical bone of L3 in OVCF models before and after PVA The results of distributions of von Mises stress under exion in cortical bone of L3 in OVCF models before and after PVA were shown in Fig. 4a. Compared with distribution before PVA, it was unchanged after PVA, which was still concentrated at posterior unfractured area. Similar results could be seen in extension, lateral bendings and axial rotations.
Regarding magnitude of maximum von Mises stress in cortical bone of L3, it decreased signi cantly after PVA under all loading conditions with that in AL decreasing the most signi cantly. Comparing AM with PL, it was lower in AM under exion, extension and axial rotations while it was lower in PL under lateral bendings (Fig. 5a).

Distributions and magnitudes of von Mises stress in cancellous bone of L3 in OVCF models before and after PVA
The results of distributions of von Mises stress under exion in cancellous bone of L3 for OVCF models before and after PVA were shown in Fig. 4b, and demonstrated that von Mises stress was concentrated at posterior unfractured area before PVA, but transferred to be concentrated at cancellous bone surrounding bone cement after PVA. However, differences were not detected among augmented models, in which von Mises stress was concentrated below and above fractured area symmetrically. Similar results could be seen in extension, lateral bendings and axial rotations.
Regarding magnitude of maximum von Mises stress in the cancellous bone of L3, it increased after PVA under exion and lateral bendings while decreased under extension (with exception of PL) and axial rotations. It was lowest in AM under all the loading conditions among AL, AM and PL. Comparing AL with PL, it was lower in AL under exion and extension while it was lower in PL under lateral bendings, the differences under axial rotations were not obvious (Fig. 5b).
Maximum displacement of L3 for OVCF models before and after PVA Regarding maximum displacement of L3, it decreased signi cantly after PVA under all loading conditions with that in AL decreasing the most signi cantly. Comparing AM with PL, it was lower in AM under exion, extension and axial rotations while it was lower in PL under lateral bendings (Fig. 5c).

Discussion
This nite element model was developed on the data collected from a spinal CT scan of a healthy lumbar spine and biomechanical material properties re ecting the pathological characteristics of vertebral osteoporosis. A validated ve-vertebra segment was constructed, but not individual vertebra, because a ve-vertebra segment model with intervertebral discs and facet joints might not only highly simulated the motion and load transfer of lumbar spine when OVCF was simulated in the middle vertebra, but also avoid loading positions and boundary conditions being directly connected to the target vertebra of L3, which may have an in uence on biomechanical behaviors of L3 [29]. The validation test proved that the constructed three-dimensional nite element model could accurately simulate physiological activity and load transfer at the lumbar spine and, therefore, could be a valuable tool for later research. Although stress measurements of cortical and cancellous bone were not validated experimentally, this could not be an in uential factor in our study because the parameters in this study were compared in terms of relative differences among fractured models before and after PVA with different distributing patterns of cement.
The anterior cortex of compressed vertebra is usually injured and not continuous before union. In the current study, we attempted to re-create a compression fracture model where a cleft was horizontally penetrated into L3 vertebral body through the center of the anterior cortical shell. This failure model may be attributed to the type A1.2 fracture, that is, the wedge impaction fracture [30]. For bipedicular augmentation, two vertically orientated cement cylinders were implanted in the vertebral body. Similar methods to simulate PVA had been reported by Polikeit et al.[31] and Liebschner et al. [32], and they also agreed that this shape was comparable with the cement distribution seen on radiographs of treated patients. We acknowledge that using cylinders to simulate PVA is done for ease of computer simulation, and not a direct re ection of the reality of PVA. Because the simulation of PVA in the current study can be easily created through nite element analysis method and ensure the repeatability of the study. We are aware that different shapes of cement might occur in the vertebral bodies. Thus, preliminary experiments, in which another shape of cement (cement cake) with equivalent volume was used to simulate PVA, had been done to analyze how cement shape affect ndings in our previous study. We found that different stress and displacement can be produced using cement cakes to simulate PVA, but the same conclusion can be drawn using cement cakes to simulate PVA like cement cylinders. In addition, because our previous study demonstrated that su cient cement distribution in the fractured area and symmetrical cement distribution around the fractured area may provide better structural support, as compared with insu cient and asymmetrical cement distribution [10]. In this study, the two cement cylinders were implanted around the fractured area symmetrically to mimic bipedicular PVA with su cient cement distribution in the fractured area and symmetrical cement distribution around the fractured area. Using this simulation method may relatively accurately re ect different biomechanical effects of AL, AM and PL in the fractured area on OVCFs.
Our current study is focus on distributions and magnitudes of von Mises stress in cortical and cancellous bone and maximum displacement of L3 for OVCF models before and after PVA. Distribution of von Mises stress in cortical bone was unchanged while that in cancellous bone was transferred to be concentrated symmetrically at cancellous bone surrounding cement after PVA, which is consistent with our previous study [10] and a nite element analysis conducted by Polikeit et al. [21]. This is likely due to the fact that increase in vertebral stability of OVCF after PVA depends on interdigitation between cement and cancellous bone. In addition, it is known that back pain caused by OVCFs is mostly likely to be related to periosteal nerves aggravated by micromotion at the fractured area and recollapse is due to reduced strength of OVCFs. Thus, stabilization of the fractured area and restoration of vertebral strength are the most acceptable mechanisms for pain relief and preventing recollapse after PVA in treatment of OVCFs.
Interestingly, we discovered that maximum displacement of L3 decreased signi cantly after PVA under all loading conditions with that in AL decreasing the most signi cantly. It indicates that vertebral stability can be better restored by AL. Furthermore, maximum von Mises stress in cortical bone of L3 had the same trend as maximum displacement of L3, which implys that AL can better restore vertebral strength in terms of stress in cortical bone. However, AM created the lowest maximum von Mises stress in cancellous bone of L3, which suggests that AM can better restore vertebral strength in terms of stress in cancellous bone. Given that both cortical and cancellous bone contribute to vertebral strength, we believe that cement distribution between AL and AM can balance stress in cortical and cancellous bone, better restoring vertebral strength. Meanwhile, since difference of maximum displacement of L3 between AL and AM was not obvious, su cient vertebral stability can be provided with cement distribution between AL and AM. Some limitations of this study need to be mentioned. Firstly, according to our previous study, location of fractured area of OVCFs could be classi ed into 3 types including superior, middle, and inferior [11]. Some OVCFs involve one part while others may involve two or three parts. The current model only focused on OVCFs with middle fractured area. Secondly, wedge-like vertebra simulating compression fracture was not created. Only the idealized status of anterior vertebral height correction was simulated in this study, however, PVA only partially restore the vertebral shape in most cases. Third, simulation of cement augmentation with cement cylinders is not a direct re ection of the reality of PVA. Above shortcomings may limit the application of the results to real world clinical scenarios. Clinical study evaluating the ndings from this study would be expected in the future.

Conclusions
While cement distributes su ciently in the fractured area and relatively symmetrically around the fractured area, various cement mass location in the vertebral body including AL, AM and PL have differences in restoring vertebral stability and strength of OVCFs due to the fact that different stress and displacement of corresponding augmented models were produced. Cement distribution between AL and AM may balance stress in cortical and cancellous bone, better restoring vertebral strength, meanwhile, providing su cient vertebral stability.  OVCF models before and after PVA.