As society ages, the number of osteoporosis patients has been increasing, turning bone fragility fractures into a major societal problem. Fractures can cause the loss of motor function, lead to other diseases, decrease quality of life, and cause elderly people to be bedridden―these are serious problems for aging societies. After osteoporosis is diagnosed, it is important to prevent fractures, or discover them early should they occur. Yet at present there is no established diagnostic method that allows the risk of fracture due to bone destruction to be evaluated quantitatively. Vertebral fractures are common in osteoporosis. While the basic treatment is conservative, balloon kyphoplasty (BKP) is often performed in cases with persistent pain. In this procedure, a balloon is inserted percutaneously into the fracture area to raise the vertebra, and then the area is filled with bone cement. However, fractures are said to occur frequently in the adjacent vertebrae after BKP. In the present study, we used the finite element method (FEM) to investigate the stress profiles of vertebrae in cases of vertebral fracture that had undergone BKP. Our purpose was to clarify the mechanism of fractures that occur in the adjacent vertebrae after BKP.
2. Materials and Methods
2.1. Patient-Specific FE Modeling
The patients provided written, informed consent prior to enrollment, and the study protocol was approved by the Ethics Committee of Juntendo University. The subjects were two patients (first case: 74-year-old woman, bone mineral density (BMD) T-score 69%; second case: 88-year-old woman, BMD T-score 64%) who underwent BKP for osteoporotic vertebral fractures (L1) at our hospital. They were consecutively recruited and studied from April 2015 to august 2015.Computed tomographic (CT) digital imaging and communications in medicine (DICOM) data was collected at 1-mm intervals after BKP was performed at L1 (Figure 1). The bone analysis software program Me- chanical Finder (Research Center of Computational Mechanics Co., Japan)  was then
Figure 1. Construction of 3D FE model. Bone geometrical features were extracted from CT DICOM by using Mechanical Finder (MF) software  .
used to construct heterogeneous 3-dimensional (3D) finite element (FE) models (T11- L3) (Figure 2). In each CT image, the regions of interest (ROI) were the outer edges of the cortical bone, which were used to create a numerical model of the anatomical structure of the vertebrae. This was also performed for the areas of bone cement HV-R in the postoperative models. The FE models were constructed from tetrahedral elements and shell elements. Tetrahedral elements with a mean size of 1 mm were used for the bone cement, inside of cortical bone, and cancellous bone. Shell elements 0.4 mm thick were used for the outside of cortical bone. On average, there were 1,500,454 and 1,566,517 tetrahedral solid elements and 191,620 and 211,230 triangular shell elements for first case model and second case model, respectively. To reflect the heterogeneity of actual bone structure, Young’s modulus and compressive yield stresses were estimated from CT image Hounsfield unit (HU) values using Keyak et al.’s experimental formula  . Based on reports by Keyak et al.  and Reilly et al.  , 0.4 was used for Poisson’s ratio. Poisson’s ratios for intervertebral discs and facet joint cartilage were 0.45 and 0.2, respectively. Young’s modulus for intervertebral discs and facet joint cartilage were 8.4 MPa and 11 MPa, respectively (Table 1). For bone cement, Young’s modulus was 3.7 MPa and Poisson’s ratio was 0.37 (Table 2)  .
The FE models were loaded with a compressive force of 1000N in four rotational/ moment loadings on the superior surface of the T11 intervertebral disc. The models were made to simulate the four physiological motions of the spine: flexion, extension,
Figure 2. Finite element models: (a) First case; (b) Second case .
Table 1. Facet joint cartilages and intervertebral disc material properties.
Table 2. Bone cement material properties  .
lateral bending, and axial rotation. The inferior side of the L3 intervertebral disc was rigidly fixed. The loading details are listed and depicted in Table 3 and Figure 3, respectively  . We referenced the boundary conditions used by Zhang et al.  , and loaded the vertebrae with 85% and the posterior components with 15% of the load. FEM analyses of the stress profiles after BKP were compared. Among the results, we focused on stress distributions on the vertebrae designations.
Figure 4 shows the Young’s modulus distributions. Young’s modulus were lower in the second case than in the first case for all vertebrae. We measured the midsagittal and midlateral plane distance the halfway point of the midsagittal plane distance of the superior and inferior endplates. Maximum Drucker-Prager stresses after BKP in the first case and second case were, respectively, 22.7 MPa and 32.1 MPa for compression, 2.7 MPa and 6.4MPa for flexion, 0.4 MPa and 1.7 MPa for extension, 11.3 MPa and 5.9 MPa for lateral bending, and 9.2 MPa and 11.3 MPa for axial rotation (Figure 5). Figure 6 shows the strain energy distributions after BKP for compression, flexion, extension, lateral bending, and axial rotation loading. Strain energy density was concentrated on the superior and inferior vertebrae (T12, L2) under all the conditions. Average strain energy densities in the first case and second case were as follows. For compression the respective values for T12 were 31.65 KJ/m3 and 101.32 KJ/m3, for L1 11.28 KJ/m3 and 13.98 KJ/m3, and for L2 18.59 KJ/m3 and 37.07 KJ/m3. For flexion the respective values
Figure 3. Loads and boundary conditions.
Table 3. Loading conditions  .
(a) (b)(c) (d)
Figure 4. (a) Young’s modulus distribution; (b) Anterior; (c) Central; (d) Posterior.
Figure 5. Maximum Drucker-Prager stress. Compression: first case 22.7 MPa, second case 32.1 MPa, flexion; first case 2.7 MPa, second case 6.4 MPa, extension; first case 0.4 MPa, second case 1.7 MPa, lateral bending; first case 11.3 MPa, second case 5.9 MPa, axial rotation; first case 9.2 MPa, second case 11.3 MPa.
Figure 6. Strain energy density distribution.
for T12 were 0.90 KJ/m3 and 2.55 KJ/m3, for L1 0.39 KJ/m3 and 1.15 KJ/m3, and for L2 0.59 KJ/m3 and 1.85 KJ/m3. For extension the respective values for T12 were 0.01 KJ/m3 and 0.14 KJ/m3, for L1 0.01 KJ/m3 and 0.07 KJ/m3, and for L2 0.01 KJ/m3 and 0.10 KJ/m3. For lateral bending the respective values for T12 were 0.78 KJ/m3 and 1.84 KJ/m3, for L1 0.80 KJ/m3 and 0.84 KJ/m3, and for L2 0.35 KJ/m3 and 1.11 KJ/m3. For axial rotation the respective values for T12 were 1.78 KJ/m3 and 0.78 KJ/m3, for L1 0.49 KJ/m3 and 0.80 KJ/m3, and for L2 1.50 KJ/m3 and 0.35 KJ/m3 (Figure 7). Note, the strain energy density of L1 was reduced and the strain energy densities of the adjacent vertebrae were increased.
Recent advancements in computational mechanics technology have enabled mechanical analyses using FEM that reflect the complex structural morphology and material characteristics of bone  . These techniques are being applied to FEM analyses of vertebral fractures. For example, 3D bone models have been constructed from DICOM data obtained from quantitative CT images. FEM has then been used to perform stress analyses to quantify strengths under external force from different directions and of different sizes  -  . Recent attempts that take individual patients’ bone strengths into consideration have used CT data to understand the pathology of and assess therapeutic effects in osteoporosis   . These techniques have also been used to evaluate surgeries using instrumentation  . One important issue for bone FEM has been how to treat the distribution of bone density in the modeling stage. Keyake et al. proposed an experimental formula that uses concentration distributions from CT images to reflect bone density  . They verified the accuracy of this formula. This formula also strongly
(a) (b)(c) (d)
Figure 7. Average strain energy density: (a) Compression; (b) T12; (c) L1; (d) L2 Compression; T12 first case 31.65 KJ/m3, second case 101.32 KJ/m3, L1 first case 11.28 KJ/m3, second case 13.98 KJ/m3, L2 first case 18.59 KJ/m3, second case 37.07 KJ/m3.
correlates with the results of vertical compression load tests on lumbar vertebral samples from fresh cadavers and FEM by using CT data  .
BKP is a superior surgical method that can quickly alleviate persistent pain due to vertebral fractures and can be safely performed in elderly people due to its low level of invasiveness. However, there have been numerous reports of intraoperative and postoperative complications with BKP. Fractures of the adjacent vertebrae after BKP are a particularly significant problem. Reported incidence rates of the adjacent vertebral fractures vary from 9.0% to 50%  . Lindsay et al. reported that the incidence rate of new secondary vertebral fractures was 19.2% from natural causes  . This indicates that BKP could be a contributing factor in increasing the incidence of vertebral fractures, especially considering that these fractures often occur within a few months of surgery.
In the present study, the CT data were from a 74-year-old woman and an 88-year-old woman, respectively. Bone density and Young’s modulus were lower in the second patient. In addition, maximum Drucker-Prager stresses were higher in the second case than in the first for compression, flexion, extension, lateral bending, and axial rotation. These findings indicate that the second patient was at higher risk of a secondary vertebral fracture. Elevated stress concentration in a vertebra increases the likelihood of damage to that vertebra, particularly in elderly people with reduced bone concentration  .
The strain energy density distributions after BKP showed increased concentrations in the adjacent vertebrae (T12, L2) for all conditions, which indicates that stress concentration increased in the vertebrae centered around the vertebra that was filled with bone cement. Bone cement increased rigidity in the filled vertebra, while elevating stress concentration in the adjacent vertebrae   .
In both cases, the average strain energy density in L1 was markedly reduced. We believe this was caused by BKP increasing the rigidity of L1 and constraining the deformation  . Reduced bone density of L1 due to stress shielding is also suggested. We believe that reduced strain energy density at L1 caused strain energy density to rise in the adjacent vertebrae, and thus elevated structural instability. Clinically, the vast majority of the adjacent vertebral fractures after BKP occur in the vertebrae immediately above and below, which is consistent with our results. The results of the present study indicate that BKP is more of a contributing factor in the adjacent vertebral fractures than are natural causes. However, BKP is a beneficial therapy that can quickly relieve pain and improve the activities of daily living. Actually, this study patients relieved pain after BKP. Further investigation into spinal alignment, vertebral height, optimal amount of bone cement, and other factors are needed to help prevent fractures in the adjacent vertebrae. The limitations of this study are a small sample size and we don’t have experimental validations. Validation study is very difficult, but we hope to try cadaver study.
We used the finite element method to quantitatively evaluate changes in stress concentration in vertebrae after balloon kyphoplasty. Our results suggest that adjacent vertebral fractures after balloon kyphoplasty are due not only to bone fragility, but also increased rigidity in the vertebrae filled with bone cement, which raises the stress concentration on the adjacent vertebrae and increases the likelihood of fracture.
Written informed consent was obtained from the patient for publication of this paper and for the use of any ac-companying images.
Conflict of Interest
The authors report no conflict of interest concerning the findings specified in this paper.
 Tawara, D., Sakamoto, J., Murakami, H., Kawahara, N., Oda, J. and Tomita, K. (2010) Mechanical Evaluation by Patient-Specific Finite Element Analysis Demonstrates Therapeutic Effects for Osteoporotic Vertebrae. Journal of the Mechanical Behavior of Biomedical Materials, 3, 31-40.
 Keyak, J.H., Rossia, S.A., Jones, K.A. and Skinner, H.B. (1998) Prediction of Femoral Fracture Load Using Automated Finite Element Modeling. Journal of Biomechanics, 31, 125-133.
 Kurtz, S.M., Villarragaa, M.L., Zhao, K. and Edidin, A.A. (2005) Static and Fatigue Mechanical Behavior of Bone Cement with Elevated Barium Sulfate Content for Treatment of Vertebral Compression Fractures. Biomaterials, 17, 3699-3712.
 Han, K.S., Rohlmann, A., Zander, T. and Taylor, W.R. (2013) Lumbar Spinal Loads Vary with Body Height and Weight. Medical Engineering and Physics, 35, 969-977.
 Zhang, L., Yang, G., Wu, L. and Yu, B. (2010) The Biomechanical Effects of Osteoporosis Vertebral Augmentation with Cancellous Bone Granules or Bone Cement on Treated and Non-Treated Vertebral Bodies: A Finite Element Evaluation. Clinical Biomechanics, 25, 166-172.
 Brekelmans, W.A., Poort, H.W. and Slooff, T.J. (1972) A New Method to Analyze the Mechanical Behavior of Skeletal Parts. Acta Orthopaedica Scandinavica, 43, 301-317.
 Bauer, J.S., Sidorenko, I., Mueller, D., Baum, T., Issever, A.S., Eckstein, F., Rummeny, E.J., Link, T.M. and Reath, C.W. (2014) Prediction of Bone Strength by CT and MDCT-Based Finite-Element-Models: How Much Spatial Resolution Is Needed. European Journal of Radiology, 83, e36-e42.
 Julienne, E.B., Floor, M.L., Van Rietbergen, B., Ito, K. and Huiskes, R. (2009) Comparison of Bone Loss Induced by Ovariectomy and Neurectomy in Rats Analyzed by In Vivo Micro-CT. Journal of Orthopaedic Research, 27, 1521-1527.
 Mazlan, M.H., Todo, M., Takano, H. and Yonezawa, I. (2014) Finite Element Analysis of Osteoporotic Vertebrae with First Lumbar (L1) Vertebral Compression Fracture. International Journal of Applied Physics and Mathematics, 4, 267-274.
 Tawara, D., Sakamoto, J., Murakami, H., Kawahara, N. and Tomita, K. (2011) Patient Specific Finite Element Analyses Detect Significant Mechanical Therapeutic Effects on Osteo-porotic Vertebrae during a Three-Year Treatment. Journal of Biomechanical Science and Engineering, 6, 248-261.
 Mazlan, M.H., Todo, M., Takano, H. and Yonezawa, I. (2015) Effect of Cage Insertion Orientation on Stress Profiles and Subsidence Phenomenon in Posterior Lumbar Interbody Fusion. Journal of Medical and Bioengineering, 5, 93-97.
 Keyak, J.H. (2001) Improved Prediction of Proximal Femoral Fracture Load Using Nonlinear Finite Element Models. Medical Engineering and Physics, 23, 165-173.
 Silva, M.J., Wang, C., Keaveny, T.M. and Hayes, W.C. (1994) Direct and Computed Tomography Thickness Measurements of the Human, Lumbar Vertebral Shell and Endplate. Bone, 15, 409-414.
 Taylor, R.S., Fritzell, P. and Taylor, R.J. (2012) Balloon Kyphoplasty in the Management of Vertebral Compression Fractures: An Updated Systematic Review and Meta-Analysis. European Spine Journal, 16, 1085-1100.
 Lindsay, R., Silverman, S.L. and Cooper, C. (2001) Risk of New Vertebral Fracture in the Year Following a Fracture. The Journal of the American Medical Association, 285, 320-323.
 Luo, J., Adams, M.A. and Dolan, P. (2010) Vertebroplasty and Kyphoplasty Can Restore Normal Spine Mechanics Following Osteoporotic Vertebral Fracture. Journal of Osteoporosis, [Epub].
 Chiang, C.K. and Wang, J.L. (2006) Strain Energy Density Variations of Osteoporotic Spine Column after Bone Cement Augmentation—An In Vitro Porcine Biomechanical Model. Journal of Biomechanics, 39, S150.