Finite element investigation osteoporotic lumbar L 1 vertebra buckling in a presence of torsional load

Olga CHABAROVA*, Vidmantas ALEKNA**, Rimantas KAČIANAUSKAS***, Oleg ARDATOV**** *Department of Strength of Materials and Engineering Mechanics, Vilnius Gediminas Technical University, Vilnius, Lithuania **Faculty of Medicine, Vilnius University, Vilnius, Lithuania ***Institute of Mechanics, Vilnius Gediminas Technical University, Vilnius, Lithuania, E-mail: rimantas.kacianauskas@vgtu.lt ****Department of Biomechanics, Vilnius Gediminas Technical University, Vilnius, Lithuania


Introduction
Osteoporosis is a disease, which is characterized by a low bone mass and bone fractures.It is one of the most widespread health problems which affects male and female; with an aging population its incidence is increasing.One out of three women, and one out of six men of age 50 and older, suffer from fractures caused by osteoporosis [1], [2].Vertebral compression fractures can cause a pain and can permanently alter the shape and strength of the spine [3], [4].Currently, osteoporotic fractures represent major cause of disability, loss of quality of life and even death amongst an elderly population [1].
From a mechanical perspective, vertebral fractures can occur due to several reasons.Evaluation of mechanical behavior of individual vertebrae is of critical importance in understanding the overall spine biomechanics.Research of mechanical properties during osteoporotic degradation demonstrate that vertebral elasticity strongly correlates with bone density decrease [5], [6], [7], [8], [9].These effects lead to visible geometrical alterations, for instance, vertebral bone shell or soft disc diminution, which directly impacts stability of spinal structure [10].Generally, mathematical models describing a mechanical instability phenomenon are time history type, and they are characterized by geometrical and/or nonlinearities.
Concept of spinal instability is determined by considering the severity of the spinal damage and adopting a management strategy, as according to various authors, unstable spinal fractures require surgical intervention.Some authors [11] suggest that the osteoporotic loss of vertebral stability are caused by cracks in the trabecular layer.After the resorption of transversal trabeculae only vertical trabeculae are maintained.Due to the lack of the transversal bonds the vertical trabeculae collapse and lose stability under the impact of the axial force.The studies of spinal torsion state analysis [12], [13] predicted fractures applying mechanical criteria of plasticity and decomposition, but the correlation between the loss of spinal stability and vertebral fractures occurrences was not included.Polikeit et al. [14] suggested imitating osteoporosis by adjusting only the properties of the necessary material.McDonald et al. [15] examined the microstructural model of the vertebra of the lumbar region to investigate the microarchitectural changes effects of the osteoporotic trabecular bone on the vertebral instability.
Recently, computational biomechanical modelling of the spine has been widely used and shown great promise in providing valuable information in treatment opportunities.The Finite Element Method (FEM) is one of the most commonly used and perspective modelling methods for the modelling of mechanical structures.It opened new opportunities for the modelling of the spine instability with the method of the finite element calculation.The modelling of the digital spine instability applying FEM was defined by Loughenbury et al, Johansen et al, Kim et al et str.[1], [16], [17].
The FE models applied to mechanical analysis of the spine range from submicron-scale assessment of vertebral bone to the analysis of the whole spine and surrounding structures.Developments of the entire human spine-scale FE models and their simulations are rather limited [17], [18].However, vertebrae are connected via intervertebral discs (IVD) and posterior elements having a relatively complex geometry.Numerical results obtained by applying a three-dimensional FE model of the human body [18] illustrate that maximal stresses are distributed along the entire skeleton, but fractures and severe deformations are concentrated in a local fragment of three vertebrae.Eight well-established FE models of the lumbar spine (L1-L5) created by different research centers were compared with in vitro and in vivo measurements [19].The graphical interpretation of the model, using original parameters of the assessed spine, provides a three-dimensional (3D) view.This digital model of the spine can then be manipulated: moved, rotated, various loads and force applied and changes to the form of the spine observed at the time of the physiological movements.
Imai [20] designed Finite Element (FE) models of the trabecular bone of the vertebrae, using 1 mm, 2 mm and 3 mm sized linear tetrahedron elements.In the above mentioned studies, all authors looked into the influence of the trabecular bone changes for vertebral stability, which emerged from compression force.These studies did not analyze the impact of torsion to the vertebral stability.
This study aimed to examine contribution of torsional load to buckling of the osteoporotic lumbar L1 human vertebral body which was axially compressed by applying the finite element method.The changes of the properties of the bone material observed in osteoporosis were modelled.Furthermore, the effect of these changes on load carrying capacity of vertebra was considered.Also a comparison between healthy and degenerated properties is presented.

Spine problem formulation and basic data analysis
The human spine is composed of bones, intervertebral discs, and ligaments.Bones that form the spine are called vertebrae.The spine can be considered as a mechanical structure constituting of relatively stiff bodies (vertebrae) coupled with relatively flexible links (discs).For many reasons the spine can be divided into the several zones, while one of them, the lumbar zone, is composed by L1-L5 vertebrae.Compression fractures of the spine usually occur at the first vertebra of the lumbar spine (L1).
Vertebral body is of cylindrical shape and it is the strongest element that consists the spine.It vertebral body plays an important role in the spineit manages varieties of loads, which occur in daily life.The vertebral body itself mainly consists of cancellous (spongy) bone and of dense material (outer thin cortical shell) (Fig. 1, b).The course part of the vertebral body is of an arc form with growths (Fig. 1, a).Vertebral body and intervertebral discs withstand approximately 75% of the total spines compression and other 25% by vertebral process [21].
Vertebra's models with geometrical sizes are generally defined by patient-specific data.After thorough analysis [22] an approximate lumbar L1 vertebrae body model containing several shape parameters was developed (Fig. 2).The height of the quasi-cylinder is approximately equal to 30 mm, while cross-sectional diameter is roughly 40 mm.
The thickness of the shell is equal to 0.5 mm.
Dense cortical shell is modelled as homogeneous isotropic elastic-perfectly plastic continuum, where the yield limit characterized by yield stress σy, while material model defined by the stress-strain diagram is given in Fig. 3. Trabecular phase is modelled as inhomogeneous elastic orthotropic continuum.Thereby transversal elasticity modulus Exx=Eyy is assumed to be fraction of longitudinal modulus Ezz, thus: The values of material parameters are shown in Table 1.
Osteoporotic Degradation of Vertebra is characterized by reduction of volume density ρ from 300 kg/m 3 up to 100 kg/m 3 leading to degradation of Young's modulus.Trabecular bone elasticity modulus is calculated according to formula given in [23] as follows: where ρ is real density.
Vertebral cortical layer and trabecular bone material physical properties are given in the Table 1.

Mathematical model
To achieve our goals formulated above, the nonlinear analysis problem is formulated.Traditionally, buckling is considered as buckling of cortical shell characterized by elastic instability due to the critical out-of-plane deformation of a structure reached under action of axial (in-plane) load.Critical state is characterized by the value of critical load at (bifurcation point).In a presence of continuum regions the critical instability may be characterized as asymptotic limit state having unlimited deformation of continuum part.For evaluation of strength limit, additional model for perfect plasticity is required.
In summary all these sub-problems may be resolved by applying the universal nonlinear model while continuum problem is replaced by a discrete finite element model.The nonlinear loading-path-dependent equilibrium is characterized by a set of nonlinear algebraic equations.Incremental formulation of this model is defined at time instant t as follows: here KG is the global nonlinear stiffness matrix comprising contribution of the finite displacements and plasticity and depending on current values of the displacement vector u(t), Δu and ΔF are increments of displacement and external load vectors.Actually, in order to find bifurcation point and to trace descending loading branch during instability of loading prescribed displacements are specified.Thus, external axial loading in time t is controlled by the instantaneous contribution of the vertical displacement uz(t), while torque ωz(t).The FE model was generated in a following manner.Vertebral cortical bone volume was discretized by shell finite elements.The shell element applied is a four-node element with six degrees of freedom at each node: translations in the x, y, and z directions, and rotations about the x, y, and z-axes.Element is associated with plasticity.It is suitable for analysing thin to moderately-thick shell structures.
A cancellous bone, endplate, vertebral process models were meshed with volumetric finite elements.This type solid element is a higher order 3-D 20-node solid element that exhibits quadratic displacement behavior.The element supports plasticity, large deflection, and large strain capabilities.
Cortical shell is bonded to trabecular bone nodeby-node connecting translational degrees of freedom.
The FE mesh of cortical shell contains 11882 nodes and 11665 shell elements.The cancellous bone, endplate, vertebral process were described by 3D mesh of 153663 nodes and 94839 solid elements.The meshed model is presented in Fig. 4.

Results obtained by numerical analysis
To evaluate contribution of the buckling type instability two basic FE models were generated.The first model contains material properties of the healthy trabecular phase, while the second model takes into consideration properties of the trabecular phase degenerated due to osteoporotic degradation.
Each of the models under action of the purely axial and the combined axial-torsional monotonic loadings was considered, respectively.The time history is defined for time interval 0 ≤ t ≤ 1.For of axial loading it is controlled by the specified displacement of the upper endplate uz(t)=uz,max, limited by maximal value uz,max=1 mm, while torque is controlled by rotational displacement (twist angle) ωz(t) limited by maximal value ωz,max=5º.
The stress-strain state of the body is obtained by solving Eq. ( 3) regarding the given loads and boundary conditions, and different selections of stresses, strains and displacements may be explored for illustration of results.
Physical nature of different models is qualitatively illustrated by deformed shapes of cortical shell is shown in Fig. 5.The displacement values defined by millimetres are illustrated in unified colour scale.It is obvious, that axial loading (Fig. 5, a, c) results into symmetrically deformed shapes while torsion (Fig. 5, b, d) is characterized by nosymmetric shift.Contribution of the osteoporotic degradation may be evaluated by significant differences in colour scale.For the case of axial load, local deformation in the frontal part at the vertical symmetry plane is observed.The torsion dramatically changes deformation character yielding increase deformation along the inclined sharing plane.
Quantitative differences between various models are illustrated in Figs.6-8, where time histories of the selected displacement and force parameters during entire loading period for 0 ≤ t ≤ 1 are plotted.The first curve denoted as (Curve a) illustrates behavior of healthy vertebrae under axial loading.Curve b illustrates behavior of the healthy vertebrae under combined axial and torsional loading.Analogously, Curve c and Curve d illustrates behavior of degenerated vertebrae.
Variation of transversal displacement of selected point A is shown in Fig. 6 while variation of von Mises stress in selected point A can be seen in Fig. 8.The point A location is chosen in advance to illustrate typical situation of cortical shell behavior.Variations of maximal carrying load are illustrated in Fig. 7. To understand the safety margin contour plot of the von Mises on the front view of cortical shall is presented in Fig. 8. Von Mises stress distributions on A points, while affected by marginal force are illustrated in Fig. 9.

Discussion and conclusions
Discussion on the obtained results is aimed to discover the role of buckling for the load bearing capacity of L1 vertebra body.Here, the lower bound of the safety margin is identified as parameter reflecting facture risk including not only classical strengths criterion but also deformation criterion reflecting large deformations leading to local instabilities.
On the basis of numerical results (Curves a) obtained for healthy vertebra body under axial load, it was found that the load carrying capacity of vertebra is characterized by strength criterion.Time histories of the von Mises stress σ (Fig. 8) show that strength criterion σy=64 MPa for the case of purely axial compression is satisfied atn time instant tax,health=0.36.Therewith, local displacement ux(tax,health)=0.066mm (Fig. 6) is relatively small, while the load at time instant tax,health (Fig. 7) may be considered as limit load Fax,health=F(tax,health)=15.7 kN.Distribution of von Mises stress (Fig. 9, a) clearly confirm this statement.
On the other hand, this load value indicates only local damage at point A, but increase of load under total load the structure is still able to resist external loading until the value of Fmax,ax,health=F(t2,ax,health)=34.5 kN obtained at time instant t2,ax,health=0.86.Distribution of stress at the time instant of failure is given in Fig. 9, b.
Behaviour of healthy vertebra under axial load combined with torsion (Curves b) may be characterised in the same manner.However, the main conclusion sounds that presence of torque reduces the load carrying capacity and increases fracture risk.
Analogously, to previous sample results also illustrate reduction of the total load carrying capacity.The maximal limit load Fmax,tors,health is reached at time instant t2,tors,health=0.51.thus Fmax,tors, health=F(t2,tors,health)=20.9 kN.
x y z Fig. 9 Distribution of the von Mises stress σ (MPa) on the cortical shell: a, bhealthy vertebra under axial load at lower limit load time instances tax,health=0.36and under total load; c, dhealthy vertebra under axial load combined with torque healthy vertebra under axial load at lower limit load time instances tax,health=0.3and under total load; e,fdegenerated vertebra under axial load at lower limit load time instances tax,health=0.78and under total load; g, hdegenerated vertebra under axial load combined with torque healthy vertebra under axial load at lower limit load time instances tax,health=0.46and under total load The other hand, this load value indicates only local damage at point A, but increase under total load the structure is still able to resist external loading.Distribution of stress at the time instant of failure is given in Fig. 9, d.
Compared to available experimental data [24], [25], manifesting carrying load varying between 10 and 15.9 kN specified for young man, the quite reliable safety margin varying from 26% up to 90%.
By considering osteoporotic degenerated lumbar, the main point is switched to displacement-based criteria.When compressing a degenerated vertebra, time history of the transversal displacement at point A (Fig. 6), exhibits unlimited character illustrating unstable deformation behaviour.Osteoporotic degradation reduces the resultant axial stiffness, therefore, occurrence of the yield limit is achieved later at time instant t2,ax,degen=0.78 (Fig. 8).It is obtained with enormous value of the transversal displacement ux(t2,ax,degen)=0.475mm (Fig. 6) indicating 428% increase when compared with an analogous points of the healthy vertebra.
As a result, the corresponding load F2,ax,degen=F(t2,ax,degen)=7.9 kN (Fig. 7) may be considered as limit load.From the strength point of view, the smooth character Curve F(t) means structure remains to withstand higher loads F(t> t2,ax,degen) > F2,ax,degen.It is obvious, that real lower bound of carrying loud should be restricted by introducing displacement constraints.Thus, the F2,ax,degen may be regarded just post-critical load.
Regarding combined axial-torsional loading the behaviour of vertebra body defined by Curves d is different.Variation of displacement (Fig. 6) clearly shows presence of the critical point at time instant ttors,degen=0.46,consequently, load bearing capacity is predefined by the critical buckling load Fmax,tors,degen=F(ttors,degen)=4.5 kN (Fig. 7).This load ischaractrised by elastic state, σ(ttors,degen)=61 MPa < 64 MPa.For limitation of transversal displacement even lower values may be required.
In summary, it could be concluded that buckling or more general instabilities play decisive role in deformation evaluation of damage of osteoporotic degradation vertebra, and local deformation criteria including buckling should be applied along with strength criteria.

Fig. 1
Fig. 1 a -A schematic view and bactual photograph of a lumbar vertebral body cut at the mid-sagittal plane, showing in particular the cortical shell, the vertebral endplates, and the porous core made of a network of trabeculae

Fig. 2
Fig. 2 Frontal section view of the model

Fig. 4
Fig. 4 External view of the finite elements model Calculation experiment was accomplished by using finite element method and software ANSYS.The FE model was generated in a following manner.Vertebral cortical bone volume was discretized by shell finite elements.The shell element applied is a four-node element with six degrees of freedom at each node: translations in the x, y, and z directions, and rotations about the x, y, and z-axes.Element is associated with plasticity.It is suitable for analysing thin to moderately-thick shell structures.A cancellous bone, endplate, vertebral process models were meshed with volumetric finite elements.This type solid element is a higher order 3-D 20-node solid element that exhibits quadratic displacement behavior.The element supports plasticity, large deflection, and large strain capabilities.Cortical shell is bonded to trabecular bone nodeby-node connecting translational degrees of freedom.The FE mesh of cortical shell contains 11882 nodes and 11665 shell elements.The cancellous bone, endplate, vertebral process were described by 3D mesh of 153663 nodes and 94839 solid elements.The meshed model is presented in Fig.4.

Fig. 5 Fig. 6 FFig. 7 Fig. 8
Fig. 5 The front view of deformed shapes of cortical shell and contour plot of total deformation after loading at time instant t=1: ahealthy vertebra under axial load; bhealthy vertebra under axial load combined with torque; cdegenerated vertebra under axial load; ddegenerated vertebra under axial load combined with torque

Table 1
Bone properties for healthy and osteoporotic body mass index (BMI) values