ARTICLE IN PRESS Journal of Biomechanics 43 (2010) 334–342
Contents lists available at ScienceDirect
Journal of Biomechanics journal homepage: www.elsevier.com/locate/jbiomech www.JBiomech.com
An accurate validation of a computational model of a human lumbosacral segment V. Moramarco a,, A. Pe´rez del Palomar b,c, C. Pappalettere a, M. Doblare´ b,c a
Department of Mechanical and Management Engineering (DIMeG), Politecnico di Bari, Italy ´n Institute of Engineering Research (I3A), Universidad de Zaragoza, Spain Group of Structural Mechanics and Materials Modelling (GEMM), Arago c ´n Biome´dica en Red. Bioingenierı´a Biomateriales y Nanomedicina (CIBER-BBN), Aragon Health Sciences Institute, Spain Centro de Investigacio b
a r t i c l e in f o
a b s t r a c t
Article history: Accepted 29 July 2009
Clinical studies have recently documented that there is sufficient evidence to suggest that abnormal motion may be an indicator of abnormal mechanics of the spine and, therefore, may be associated with some types of low-back pain. However, designating a motion as abnormal requires knowledge of normal motions. This work hence aims to develop an accurate computational model to simulate the biomechanical response of the whole lumbosacral spinal unit (L1–S1) under physiological loadings and constraint conditions. In order to meet this objective, computed tomography (CT) scanning protocols, finite element (FE) analysis and accurate constitutive modelling have been integrated. Then the ranges of motion (ROM) under flexion, extension and lateral bending moment were measured and compared with experimental data, finding an excellent agreement. In particular, the ability of the model to reproduce the relative rotation between each couple of vertebrae was proved. Finally, the shear stresses for the most extreme load cases were reported in order to predict which are the most risky conditions and where the maximum damage would be located. The results indicate that the greater values of the stresses were located at L4–S1 levels just in the interfaces between disc and vertebrae across the posterior and posterolateral zone. This result can be clinically correlated with the existence of damage exactly where the stresses were maximal in the proposed finite element model. & 2009 Elsevier Ltd. All rights reserved.
Keywords: Lumbar spine Disc degeneration Fibre reinforced material Collagen fibres Lumbar motion Finite element method
1. Introduction Low back pain is one of the most diffuse chronic pathologies and represents one of the highest direct and indirect cost for national welfare. It can affect for years the patient with the obvious disabling consequences and many times without a complete explanation of the causes (Devor and Tal, 2009). Some authors like Deyo and Weinstain (2001) explain that idiopathic low back pain may be originated by the mechanical degeneration of some rachial structure (discs, ligaments, facet joints, etc.). Moreover, other clinical studies (Dvorak et al., 1991; Parnianpour et al., 1988; Percy et al., 1984; Percy and Tibrewal, 1984) show a relation between the anomalous, larger or smaller than normal relative displacement, of each couple of vertebrae segment and low-back pain. In most cases, surgery is the only real solution to solve the problem. In spite of the more extended use of clinical imaging analysis like computed tomography (CT) scan or magnetic resonance imaging (MRI) that allows the clinician to see an
Corresponding author.
E-mail address:
[email protected] (V. Moramarco). 0021-9290/$ - see front matter & 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.jbiomech.2009.07.042
evidently degenerate situation, it is not possible to predict a specific behaviour of the rachis. Furthermore, the association between symptoms and imaging results is weak (White and Gordon, 1982). For these reasons, usually, orthopaedic therapy is essentially based on the experience of the surgeon who predicts the best solution for each patient. Using mathematical models and computer simulations could potentially be therefore an important tool to support clinical decisions. From the analyses of the motion behaviour of the patient’s lumbar segment we could be able to diagnose an abnormal situation and therefore try to design and establish better solutions. In the last years numerous FE models of lumbosacral segment have been created but most of them were focused on effects of a specific spinal prosthesis (Donezier and Ku, 2006; Rohlmann et al., 2005) or on a specific clinical problem (Boccaccio et al., 2008; Polikeit et al., 2004). On the contrary, not so many of them were specifically designed to simulate the behaviour of healthy lumbosacral segment. Some of these models like Shirazi-Adl (1994) and Ezquerro et al. (2004) were able to predict the rotation for a specific value of moment but they did not consider the extreme non-linearity of the spine behaviour. In some others, like Eberlein et al. (2004) and Zander et al. (2009), a very accurate modelling of each part of the entire lumbar segment was
ARTICLE IN PRESS V. Moramarco et al. / Journal of Biomechanics 43 (2010) 334–342
presented, but only the rotation of the entire segment with respect to the applied moment was compared with experimental results and no results about the relative rotation between each vertebrae pair under loading were reported. In Guan et al. (2006) the relative rotations are shown but the model was limited to the L4–S1 segment and the anisotropy of the intervertebral disc behaviour was not considered. The aim of the present study is to demonstrate that it is possible to create a FE model of the entire healthy lumbosacral segment that, not only is able to predict the relative rotations between each vertebrae pair under a specific load, but can simulate very accurately the non-linearity of these movements. For this purpose, a FE model of the complete L1–S1 segment of the human rachis, based on CT-scan imaging and experimental data, has been built in order to simulate the non-linear behaviour of the rachis measured in vitro by other authors (Guan et al., 2007; Panjabi et al., 1994). Finally, since some recent studies (Costi et al., 2007; Wilson et al., 2003) have shown how the shear effects are dominant in disc tissue failures, the maximum shear stresses in the discs have been studied for flexion, extension and lateral bending moments.
2. Materials and methods To reconstruct the FE model of the healthy lumbar segment, CT scan data of a subject with no current spinal pathology were used. The domain boundary of each s vertebra was created individually, using the built-in module MedCAD available in the software application Mimicss 10.0 by Materialise Inc. No smoothing operations were carried out on the surfaces in order to preserve the original geometry as well as possible. Then, the vertebrae were discretized using surface triangular 3-nodes elements and reassembled preserving the original positions. In all conducted FE simulations bones were modelled as rigid bodies; in fact, they are much stiffer than the soft tissues present in the model and their deformation can be consider negligible in comparison with those undergone by the latter. The intervertebral discs were modelled with deformable elements. The upper and the lower surfaces of each disc were defined using the surface of the vertebrae above and below respectively. The software I-Deas 9 (Altair Computing Inc.)
335
(Structural Dynamics Research Corporation, 1993) was used to generate these surfaces. Then, using Cubit 10.1 (Leland, 2001), the volume domain was generated and the annulus fibrosus and the nucleus pulposus were drawn so that volumetric ratio was 3:7 like reported in Goto et al. (2002). The domain was discretized with 8node hexahedral elements. Finally all disc meshes were adjusted to the related vertebrae surfaces in order to obtain a perfect continuity and to remove interferences between the bodies (Fig. 1). The ligaments, anterior longitudinal (ALL), posterior longitudinal (PLL), intertransverse (ITL), flavum (LF), capsular (JC), interspinous (ISL) and supraspinous (SSL) were modelled as tension only nonlinear truss elements. The mechanical properties, listed in Table 1, were based on those presented in literature (Pintar et al., 1992; Chazal et al., 1985; Goel et al., 1995; Goto et al., 2003). Finally, all parts were assembled in order to obtain the L1–S1 entire model (Fig. 2). The constitutive equation for modelling the mechanical behaviour of the intervertebral disc was implemented in a user defined UMAT subroutine for ABAQUS. A hyperelastic fibre-reinforced model with two families of fibres (Holzapfel, 2000) was used to simulate the anisotropic behaviour of this soft tissue (Eberlein et al., 2001; Pe´rez del Palomar et al., 2008). The strain energy function used for this hyperelastic material is given by
C ¼ C10 ðI 1 3Þ þ C01 ðI 2 3Þ þ C20 ðI 1 3Þ2 þ C02 ðI 2 3Þ2 K1 fexp½K2 ðI4 1Þ2 1g 2K2 K1 1 þ fexp½K2 ðI6 1Þ2 1g þ ðJ 1Þ2 D 2K2 þ C11 ðI 1 3ÞðI 2 3Þ þ
ð1Þ
where Cij are material constants related to the ground substance, K1 4 0 and K2 40 are the parameters which identify the exponential behaviour due to the presence of collagen fibres of both families (note that here to simplify the model, the same response was assumed for both families) and D identifies the tissue incompressibility modulus. The invariants are defined as I 1 ¼ trC ;
2
I 2 ¼ 12½ðtrC Þ2 trðC Þ
I 4 ¼ a0 C a0 ;
0
I6 ¼ b C b
0
ð2Þ
0
where a is the unitary vector defining the orientation of the first family of fibers, 0 b the direction of the second family both in the reference configuration and C is the modified right Green strain tensor: T
C ¼F F
with F ¼ ðJÞ1=3 F
being F the deformation gradient and J ¼ detðFÞ.
Fig. 1. L4–L5 intervertebral disc volume definition: vertebrae surfaces drawing (a), structurated mesh discretization (b) and bodies interferences removal (c,d).
ð3Þ
ARTICLE IN PRESS 336
V. Moramarco et al. / Journal of Biomechanics 43 (2010) 334–342
Each family of fibres was placed circumferentially around the nucleus with an angle orientation a of 7 303 with respect to the disc plain (Goel et al., 1995). The values of the aforementioned constants were set such the difference between the
displacements predicted by the numerical model of the rachis segment L4–L5 and those measured by Markolf and Morris (1974) were minimal (Fig. 3a). Moreover the magnitude of the found values for the material constants was in the range of
Table 1 Material properties, geometrical parameters and elements number of the lumbosacral spine ligaments (Chazal et al., 1985; Goel et al., 1995; Goto et al., 2003; Pintar et al., 1992). Ligament
E1 (MPa)
E2
e12 a
Elements
Area ðmm2 Þ
ALL PLL LF ITL SSL
7.8 1.0 1.5 10.0 3.0
20.0 2.0 1.9 59 5.0
0.12 0.11 0.062 0.18 0.20
5 5 3 4 3
32.4 5.2 84.2 1.8 25.2
Ligament
Spine level
Area ðmm2 Þ
Poisson’s ratio n
Elements
Stiffness k (N/mm)
JC
L1–L2 L2–L3 L3–L4 L4–L5 L5–S1 L1–L2 L2–L3 L3–L4 L4–L5 L5–S1
43.8
0.4
6
35.1
0.4
6
42:5 7 0:8 33:9 7 19:2 32:3 7 3:3 30:6 7 1:5 29:9 7 22:0 10:07 5:2 9:6 7 4:8 18:1 7 15:9 8:7 7 6:5 16:3 7 15:0
ISL
a
Note: e12 denotes the strain transition between the two bilinear moduli (E1 and E2 ).
Fig. 2. Complete finite element model of the spinal segment L1–S1: (a) front view and (b) lateral view.
ARTICLE IN PRESS V. Moramarco et al. / Journal of Biomechanics 43 (2010) 334–342 the available data in the literature (Ebara et al., 1996; Eberlein et al., 2001; Zander et al., 2009). The goodness of the found data was then verified comparing the posterior disc bulge in our model with the experimental values measured by Brown et al. (1957) (Fig. 3b). To define the nucleus pulposus it was considered that the behavior of this part was isotropic and almost incompressible, so it was modelled as a hyperelastic NeoHookean material with material constants C10 ¼ 0:16 MPa and D ¼ 0:024 MPa1 (Pe´rez del Palomar et al., 2008). Contact constrain was activated between the facets of the posterior processes using GAP elements. Table 2 sums up all mechanical and element type properties of the model.
337
In order to validate this numerical model we reproduced in silico some in vitro experimental results available in literature, where the behaviour of complete lumbosacral-spine specimens, dissected from all non-ligamentous soft tissue, were analysed under specific loads. For this purpose, the lower face of the vertebra S1 was fully constrained. The loads were applied on the most proximal vertebra in our case L1. Four different finite element analyses were performed using ABAQUS 6.5 (Fig. 4). Firstly, according to Guan et al. (2007), an incremental 4 Nm pure flexion moment was applied on the model so that the nonlinear response of the spine under low loads was tested. Then, to simulate the tests proposed by Panjabi et al. (1994), an axial compressive pre-loading of 100 N was applied on the model, such that the rotation of the vertebra L4 was minimized, and two loads cases, respectively, a flexion and an extension moment both of 10 Nm, were applied incrementally to prove the nonlinear behaviour of the model in the entire flexion– extension range of the movement. Finally, in the fourth simulation carried out the response of the model under a lateral bending moment of 10 Nm was also studied.
3. Results First the validation of the model was proved showing its capacity to reproduce the nonlinear behaviour of a real spine under different loading conditions. In particular, in Fig. 5 the model response to an increasing pure flexion moment of 4 Nm is shown. This was compared with the values measured by Guan et al. (2007). The trends of the numerical and experimental curves showed a very good agreement with the whole range of motion (ROM) predicted by the FE model falling within the experimental standard deviation interval. The maximum found error between experimental data and our numerical results for a flexion moment of 4 Nm was smaller than 0:33 (L2–L3: exp. 2:83 , FE 2:53 ; L3–L4: exp. 3:43 , FE 3:53 ; L4–L5: exp. 4:43 , FE 4:33 ). Also in agreement with the experimental measures, the highest ROM was found in the segment L5–S1 (exp.: 6:83 , FE: 4:43 ) and the smallest in L1–L2 (exp.: 2:53 , FE: 1:63 ). Moreover, in the same segments, the highest disagreement with respect to the mean value found was 30% for the maximum applied moment.
Fig. 3. Comparison between the numerical model and the experimental data: (a) compressive axial force the L4–L5 relative displacement and (b) posterior disc bulge versus the compressive axial force.
Fig. 4. Schematic representation of load and boundary conditions applied to the numerical model. In Guan et al. (2007) set-up no axial load was applied.
Table 2 Material properties and element types (Hibbit and Sorensen, 2004) in the finite element model. Component
Element type
Elements
Bone Nucleus
R3D3 C3D8
73016 10106
Annulus
C3D8
19421
Ligaments Posterior processes contact
T3D2 GAP
180 30
Material constants
C10 ¼ 0:16 MPa D ¼ 0:024MPa 1 C10 ¼ 0:1 MPa C20 ¼ 2:5MPa K1 ¼ 1:8MPa K2 ¼ 11 D ¼ 0:306MPa 1 Table 1
Notes Rigid body Hyperelastic Neo-Hookean (Pe´rez del Palomar et al., 2008) Hyperelastic fibre reinforced material ða ¼ 7 303 Þ (Eberlein et al., 2001)
Tension only truss
ARTICLE IN PRESS 338
V. Moramarco et al. / Journal of Biomechanics 43 (2010) 334–342
In Fig. 6, instead, the reaction under a flexion–extension of 10 Nm moments was studied and related with the results presented by Panjabi et al. (1994). In this case too, the model
displays a good concordance with the experimental curves always remaining within the standard deviation interval, except for L5–S1 segment in extension and L2–L3 in flexion. In flexion the
Fig. 5. Flexion results for L1–S1 and comparison with Guan et al. (2007).
ARTICLE IN PRESS V. Moramarco et al. / Journal of Biomechanics 43 (2010) 334–342
Fig. 6. Comparison between the relative rotation angle under flexion–extension moment predicted by the FE model and that measured by Panjabi et al. (1994).
339
ARTICLE IN PRESS 340
V. Moramarco et al. / Journal of Biomechanics 43 (2010) 334–342
differences between experimental results and numerical prediction for 10 Nm were in the range from 0:53 (L5–S1 level) to 1:53 (L1–L2 level), while in extension they were in the interval from 0:053 (L3–L4 level) to 2:23 (L5–S1 level). According to experimental data, the ROM for the maximum flexion moment of the segments L4–L5 and L5–S1, 8:63 and 7:73 respectively, resulted substantially higher than for the other segments. Under extension moment the model fits very well the experimental curves, except for the L5–S1 segment, as it will be explained in Discussion section. The maximum ROM in extension is almost similar in all segments ð 33 Þ and quite lower than in flexion, as expected. As final validation, the ROM for a 10 Nm flexo-extension (Fig. 7) and lateral bending moments (Fig. 8) is shown and compared with experimental (Panjabi et al., 1994) and previous numerical (Shirazi-Adl, 1994) results. In both cases it is remarkable how the present model has, differently from previous finite element models, the same distribution of the ROM between the segments with the exception of the L5–S1 segment. This can be attributed to several specific causes that will be discussed later. In subjects not afflicted with bone tissue related pathologies, disc injury within physiological load represents one of the main causes of the alteration of the spine behaviour and, therefore, pain. In particular, the interface between intervertebral discs and vertebral surfaces represents the most critical region for the onset of spinal diseases. Moreover, as mentioned in Introduction section, some recent studies (Costi et al., 2007; Wilson et al., 2003) have shown how the shear effects are dominant in disc tissue failures. For this reasons, it was considered interesting to plot, for the more stressed disc, in Panjabi configuration (Panjabi et al., 1994), the maps of the maximum shear stress for flexion, extension and lateral bending moments (Figs. 9–11). These results show that, in every loading configuration, the most critical disc/vertebra interface is located between L5–S1 disc and S1 vertebra. Furthermore, as expected, the most risky load condition corresponds to the flexion moment. In this case, in fact, the shear stress had a maximum value of 30 MPa. In lateral bending, on the other hand, both upper and lower interfaces present high values (20 MPa) of shear stress, at the inner face of the spine. On the contrary, the extension condition, as expected, seems to be less problematic than the others with a maximum shear stress of 4 MPa.
Fig. 7. Comparison between the relative rotation angle under a 10 Nm flexion– extension moment predicted by the FE model and that measured by Panjabi et al. (1994) and predicted by Shirazi-Adl (1994).
4. Discussion The goal of this work was to demonstrate that starting from a CT scan of an apparently healthy human lumbosacral segment it is possible to create a FE model able to mimic the real behaviour of the lumbosacral spine. To demonstrate the goodness of this model some in vitro tests were reproduced in silico under similar conditions (no muscles were considered; just the lower vertebra fully constrained; loads applied on the upper vertebral body).
Fig. 8. Comparison between the relative rotation angle under a 10 Nm lateral moment predicted by the FE model and that measured by Panjabi et al. (1994) and predicted by Shirazi-Adl (1994).
Fig. 9. Maximum shear stress distribution under a flexion moment for disc L5–S1.
ARTICLE IN PRESS V. Moramarco et al. / Journal of Biomechanics 43 (2010) 334–342
Fig. 10. Maximum shear stress distribution under a extension moment for disc L5–S1.
Fig. 11. Maximum shear stress distribution under a lateral bending moment for disc L5–S1.
The annulus fibrosus was simulated as a fibre reinforced material. The related coefficients were set in order that the L4-L5 mono-segment axial displacement under a static increasing load in accordance with that reported by Markolf and Morris (1974). Then, to prove and consolidate the selected constants, it was shown that the model predicted the same posterior bulge measured by Brown et al. (1957), proving the biaxiality of the annulus material behaviour. The precedent models of the lumbosacral spine presented by Eberlein et al. (2004) and by Zander et al. (2009) compared only the kinematic response of the whole lumbar spine under loads with experimental data. In Shirazi-Adl (1994), instead, only the behaviour for a specific load value was considered (e.g. a 10 Nm flexion moment). In this work the capacity of the model to predict the relative rotation between each couple of vertebrae under
341
different loads was shown also including their intrinsic nonlinearity. The results obtained under the same loading conditions showed a very good agreement with Guan et al. (2006) experimental data. They were almost perfect for the L4–L5 segment, with a maximal error of 3% compared with the mean value of the experimental results. These results compared with those shown by Guan et al. (2007), where the flexion behaviour appeared almost linear, would prove the importance of using a nonlinear material behaviour for ligaments and a nonlinear anisotropic material behaviour for the annulus fibrosus. In segments L1–L2, the most cephalic, and L5–S1, the most caudal, the concordance was not so good. This can be explained by the fact that both loads and constraints were applied in those segments. Another explanation could be related with the ligaments behaviour. We have to take into account that ligaments are rather different from upper to lower levels with different area and mechanical properties. However, as a first approach and because of the lack of data in the literature, we used the same definition at all levels, with exception of JC and ISL ligaments. These parameters represent good mean values for L2–L5 ligaments but not for L1–L2 and L5–S1. In comparison with Panjabi et al. (1994) experimental values, we can observe that the FE model presents some interesting similarities. The spinal behaviour in extension and in flexion is asymmetric in agreement with literature. The average highest ROM for each segment in flexion ð6:33 Þ is usually much bigger than in extension ð3:13 Þ. Furthermore while in extension the maximal rotation value has little variation between upper (3:53 for L5–S1 level) and lower (2:63 for L3–L4) segments, in flexion this range is much bigger from 4:23 , in L1–L2, to 8:63 , in L4–L5. However, as mentioned before, the results for L5–S1 are not the ones expected. This is related with the fact that the posterior processes are closer than in other normal spines which implies a premature contact reducing the maximal reachable rotation as can be seen in Fig. 6. This effect is pointed out, also, in the results reported in Fig. 7, where the relative rotation angle between each couple of vertebrae under a 10 Nm flexion-extension moment of the model was compared with experimental (Panjabi et al., 1994) and previous numerical (Shirazi-Adl, 1994) data. It is important to remark that the present model improves the previous developed ones since it is not only able to predict a ROM within the range of variation of experimental data for each spinal level, but also presents the same distribution between the different segments (excluding L5–S1 due to the anomalous behaviour in extension). This is also the case for lateral bending. In this case, we also noticed an unusual behaviour of the L5–S1 in comparison with the experimental results. A possible explanation may be, in our opinion, the variation in flexion-extension behaviour, in the ligaments behaviour, and in particular in the intertransverse ligament (ITL) that, in this movement, is the most important. After validating the model regarding its kinematical behaviour, the stress distribution in the discs for the different loading patterns was analysed. One of the most common pathologies of the intervertebral disc is herniation. Several studies have shown how the maximum shear stresses may be related with the initiation of soft tissue failure (Costi et al., 2007). For this reason, we also studied the response of our spinal model regarding the maximum shear stresses under different physiological loads. We observed that the most serious condition corresponds to flexion and to lateral bending. In both cases the most stressed discs were L5–S1 and L4–L5; but while in flexion moment the peak value was found in the posterior part, for lateral bending it was located in the posterolateral part of the disc. These results are supported by clinical studies (Latorraca and Forni Niccolai Gamba, 2004) that show that the greatest proportion of lumbar hernias (about 90%)
ARTICLE IN PRESS 342
V. Moramarco et al. / Journal of Biomechanics 43 (2010) 334–342
are located in these segments (L5–S1 and L4–L5) and in particular in the posterior and posterolateral parts of the disc. Other studies (Lotz et al., 1998) also prove that the disc failure starts, usually, at the interface between disc and vertebra. This can be also explained by our results since the most stressed areas appeared at those interfaces mainly due to the hard discontinuity between these two materials. In conclusion, the present model can be considered as a useful tool to predict the stresses undergone by the discs under different loading conditions, and may be considered as a first step for future patient-specific models to predict the appearance and evolution of different pathologies and for preoperative planning and implant design.
Conflict of interest statement There is not any financial and personal relationship with other people or organizations that could inappropriately influence this work. References Boccaccio, A., Vena, P., Gastaldi, D., Franzoso, G., Pietrabissa, R., Pappalettere, C., 2008. Finite element analysis of cancellous bone failure in vertebrae body of healthy and osteoporotic subjects. Proceedings of the Institution of Mechanical Engineers [H] 222, 1023–1036. Brown, T., Hansen, R.J., Yorra, A.J., 1957. Some mechanical test on the lumbosacral spine with particular reference to intervertebral disc. Journal of Bone and Joint Surgery 39, 1135–1164. Chazal, J., Tanguy, A., Bouges, M., Gaurel, G., Escande, G., Guillot, M., Vanneuville, G., 1985. Biomechanical properties of spinal ligaments and an histological study of the supraspinal ligament in traction. Journal of Biomechanics 18, 167–176. Costi, J.J., Stokes, I.A., Gardner Morse, M., Laible, J.P., Scoffone, H.M., Iatridis, J.C., 2007. Direct measurement of intervertebral disc maximum shear strain in six degrees of freedom: motion that place disc tissue at risk injury. Journal of Biomechanics 40, 2457–2466. Donezier, G., Ku, D.N., 2006. Biomechanical comparison between fusion of two vertebrae and implantation of an artificial intervertebral disc. Journal of Biomechanics 39, 766–775. Deyo, R.A., Weinstain, J.N., 2001. Low back pain. New England Journal of Medicine 344, 363–370. Devor, M., Tal, M., 2009. What causes low back pain?. Pain 142, 11–12. Dvorak, J., Panjabi, M.M., Chang, D.G., Theiler, R., Grob, D., 1991. Functional radiographic diagnosis of the lumbar spine. Flexion–extension and lateral bending. Spine 16, 562–571. Ebara, S., Iatridis, J.C., Setton, L.A., Foster, R.J., Mow, V.C., Weidenbaum, M., 1996. Tensile properties of nondegenerate human lumbar annulus fibrosus. Spine 21, 452–461. Eberlein, R., Holzapfel, G.A., Schulze-Bauer, C.A.J., 2001. An anisotropic model for annulus tissue and enhanced finite element analysis of intact lumbar disc bodies. Computer Methods in Biomechanics and Biomedical Engineering 4, 209–229. ¨ Eberlein, R., Holzapfel, G.A., Frohlich, M., 2004. Multi-segment FEA of the human lumbar spine including the heterogeneity of the annulus fibrosus. Computational Mechanics 34, 147–165. Ezquerro, F., Simo´n, A., Prado, M., Pe´rez, A., 2004. Combination of finite element modelling and optimization for the study of lumbar spine biomechanics considering the 3D thorax–pelvis orientation. Medical Engineering & Physics 26, 11–22.
Goel, V.K., Monroe, T., Gilbertson, L.G., Brinckmann, P., 1995. Interlaminar shear stresses and laminae separation in a disc. Spine 20, 243–246. Goto, K., Tajima, N., Chosa, E., Totoribe, K., Kuroki, H., Arizumi, Y., Arai, T., 2002. Mechanical analysis of the lumbar vertebrae in the three-dimensional finite element method moments in which intradiscal pressure in the nucleus pulposus was used to establish the model. Journal of Orthopaedic Science 7, 243–246. Goto, K., Tajima, N., Chosa, E., Totoribe, K., Shinichiro, S., Kuroki, H., 2003. Effect of lumbar spinal fusion on the other lumbar intervertebral levels (three dimensional finite element analysis). Journal of Orthopaedic Science 8, 577–584. Guan, Y., Yoganandam, N., Zhang, J., Pintar, F.A., Cusick, J.F., Wolfla, C.E., Maiman, D.J., 2006. Validation of a clinical finite element model of the human lumbosacral spine. Medical and Biological Engineering and Computing 44, 633–641. Guan, Y., Yoganandam, N., Moore, J., Pintar, F.A., Zhang, J., Maiman, D.J., Laud, P., 2007. Moment-rotation responses of the lumbosacral spinal column. Journal of Biomechanics 40, 1975–1980. Hibbit, K., Sorensen, Inc., 2004. ABAQUS Analysis User’s Manual, v. 6.5. HKS inc., Pawtucket, RI, USA. Holzapfel, G.A., 2000. Nonlinear Solid Mechanics. Wiley, New York. Latorraca, A., Forni Niccolai Gamba, C., 2004. Analysis on 149 consecutive cases of intervertebral disc prolapse operated with microendoscopic (MetrX) technique. Reumatismo 56, 31–35. Leland, R., 2001. Cubit mesh generation environment. Version 10.0, vol. 1 edition, CUBIT Development Team, Sandia National Laboratory, Albuquerque. Lotz, J.C., Colliou, O.K., Chin, J.R., Duncan, N.A., Liendenberg, E., 1998. Compressioninduced degeneration of the intervertebral disc: an in vivo mouse model and finite-element study. Spine 23, 2493–2506. Markolf, K.L., Morris, J.M., 1974. The structural components of the vertebral disc. A study of their contribution to the ability of the disc to withstand compressive force. Journal of Bone and Joint Surgery 56, 675–687. Panjabi, M.M., Oxland, T.R., Yamamoto, I., Crisco, J.J., 1994. Mechanical behavior of the human lumbar and lumbosacral spine as shown by three-dimensional load–displacement curves. American Journal of Bone and Joint Surgery 76, 413–424. Parnianpour, M., Nordin, M., Kahanovitz, N., Frankel, V., 1988. The triaxial coupling of torque generation of trunk muscles during isometric exertions and the effect of fatiguing isoinertial movements on the motor output and movement patterns. Spine 13, 982–992. Percy, M., Portek, I., Shepherd, J., 1984. Three-dimensional X-ray analysis of normal movement in the lumbar spine. Spine 9, 294–297. Percy, M., Tibrewal, S.B., 1984. Axial rotation and lateral bending in the normal lumbar spine measured by three-dimensional radiography. Spine 9, 582–587. Pe´rez del Palomar, A., Calvo, B., Doblare´, M., 2008. An accurate finite elements model of the cervical spine under quasi-static loading. Journal of Biomechanics 41, 523–531. Pintar, F.A., Yoganandam, N., Myers, T., Elhagediab, A., Sances Jr., A., 1992. Biomechanical properties of human lumbar spine ligaments. Journal of Biomechanics 25, 1351–1356. Polikeit, A., Nolte, P., Ferguson, S.J., 2004. Simulated influence of osteoporosis and disc degeneration on the load transfer in a lumbar functional spinal unit. Journal of Biomechanics 37, 1061–1069. Rohlmann, A., Zander, T., Bergmann, G., 2005. Effect of total disc replacement with ProDisc on the biomechanical behavior of the lumbar spine. Spine 30, 738–743. Shirazi-Adl, A., 1994. Biomechanics of the lumbar spine in sagittal/lateral moments. Spine 19, 2407–2414. Structural Dynamics Research Corporation, Sandia National Laboratory, Albuquerque, 1993. I-DEAS Manuals: I-DEAS Master Series, Eastman dr., Milford, Ohio 45150 edition. White III, A.A., Gordon, S.L., 1982. Synopsis: workshop on idiopathic low-back pain. Spine 7, 141–149. Wilson, W., van Rietbergen, B., van DonKelaar, C.C., Huiskes, R., 2003. Pathways of load induced cartilage damage causing cartilage degeneration in knee after meniscectomy. Journal of Biomechanics 36, 845–851. Zander, T., Rholmann, A., Bergmann, G., 2009. Influence of different artificial disc kinematics on spine biomechanics. Clinical Biomechanics 24, 135–142.
