The Effect of Muscle Direction on the Predictions of Finite Element Model of Human Lumbar Spine

The normal physiological loads from muscles experienced by the spine are largely unknown due to a lack of data. The aim of this study is to investigate the effects of varying muscle directions on the outcomes predicted from finite element models of human lumbar spine. A nonlinear finite element model of L3–L5 was employed. The force of the erector spinae muscle, the force of the rectus abdominis muscle, follower loads, and upper body weight were applied. The model was fixed in a neural standing position and the direction of the force of the erector spinae muscle and rectus abdominis muscle was varied in three directions. The intradiscal pressure, reaction moments, and intervertebral rotations were calculated. The intradiscal pressure of L4-L5 was 0.56–0.57 MPa, which agrees with the in vivo pressure of 0.5 MPa from the literatures. The models with the erector spinae muscle loaded in anterior-oblique direction showed the smallest reaction moments (less than 0.6 Nm) and intervertebral rotations of L3-L4 and L4-L5 (less than 0.2 degrees). In comparison with loading in the vertical direction and posterior-oblique direction, the erector spinae muscle loaded in the anterior-oblique direction required lower external force or moment to keep the lumbar spine in the neutral position.


Introduction
Finite element (FE) analysis has been used extensively as a powerful computational tool for resolving difficult or complex clinical and biomechanical situations [1]. The quality of the outcome is strongly influenced by the anatomical geometry [2], material properties [1], and boundary and loading conditions [3]. The anatomical geometry is typically generated from medical images of a representative subject and the material properties are often measured from cadaveric specimens. However, given the highly invasive nature of recording in vivo data, the normal physiological loading in the spine is still largely unknown due to a lack of data on muscle forces [4].
Various assumptions have been made in experimental and numerical studies. Compressive force or/and pure moment was used earlier [4]. Recently, the concept of using a single compressive follower load was introduced to represent body weight and muscles [3,5,6]. The follower load is a compressive load applied along a follower load path that approximates a tangent to the curve of the lumbar spine. For posture, the musculature is one of the most important components for preserving spinal stability [7,8]. Many studies have attempted to incorporate muscle forces into FE simulations in order to improve the accuracy of the model. Kong and Goel [9] developed an optimization-based force model with experimental input to predict muscle forces in a single-joint FE model. Some authors estimated muscle forces using a kinematics-driven algorithm [10,11] and an EMG-assisted optimization method [10,12]. Calisse et al. [13] and Zander et al. [14] recorded in vivo data from internal spinal fixators to predict muscle forces. Khurelbaatar et al. [15] developed a spinal FE model with trunk muscles and estimated muscle forces using an optimization technique.
However, it was difficult to share the above-mentioned muscle forces for other researchers due to their private 2 BioMed Research International  [20] Ground substance of annulus fibrosis Hyperelastic, neo-Hookean, C 10 = 0.3448, D 1 = 0.3 [22] Cartilage of facet joint Soft contact [23,24] Ligaments Nonlinear [20,25] optimization techniques. We have previously carried out sensitivity studies when transferring muscle forces from a shared inverse static model to load an FE model [16,17]. The findings showed that many parameters have to be considered when combining two biomechanical models and artificial errors may be induced in the model which may have a significant impact on the outcome. Rohlmann et al. [26] considered four groups of loads: the force of the erector spinae (ES) muscle, the force of the rectus abdominis (RA) muscle, follower loads, and upper body weight. The magnitude of loads and action points were estimated from representative studies and data, but all muscle forces were applied in a vertical direction. This is an oversimplification of the true loading condition in the body. In our opinion, more realistic muscle directions may lead to more physiological loads. FE analysis regarding mechanism for some diseases and biomechanical efficacy for surgery may be beneficial. Therefore, the aim of this study is to investigate the effects of muscle direction on the simulation outcomes during FE modeling of lumbar spine.

Finite Element Model of L3-L5.
A nonlinear FE model of L3-L5 [27] was employed for this study. The vertebral bodies of the lumbar spine were modeled as isotropic cortical shells with a thickness of 0.5 mm [28,29], also with a transverse isotropic cancellous core and isotropic posterior bony structures. Cartilaginous endplates were assigned a thickness of 0.8 mm [30,31]. The intervertebral discs were comprised of an incompressible nucleus pulposus surrounded by a composite annulus fibrosus. The annuli fibrosi were modeled as fiberreinforced hyperelastic composites. Rebar elements of two times seven layers were used to represent the fibers and the stiffness was increased as the fibers radiated out from the center [21]. A 0.25 mm thick cartilaginous layer was also added to each facet articular surface, and a gap of 0.5 mm was created between the curved facet joints [23]. All seven ligaments of the lumbar spine were represented by tensiononly spring elements with nonlinear material properties [20,25]. The material properties of the various tissues in the FE model are listed in Table 1.

Validation of Mobility.
In order to validate the mobility of the model, 7.5 Nm moments were applied to the top of L3 to simulate flexion, extension, right lateral bending, and right axial rotation, while all six degrees of freedom of the inferior endplate of L5 were rigidly fixed. The intervertebral rotation (IVR) of L3-L4 and L4-L5 was calculated and compared with in vitro experiment data [18].

Loading and Evaluation.
Four groups of loads [26] were applied to the models in this study. The global back muscle forces were represented by the force of ES with a magnitude of 170 N acting 40 mm dorsal to the related disc center. Similarly, the force of RA was used to represent the global abdominal muscle forces. The magnitude was set at 20 N and acted 153 mm ventral to the related disc center. The action line of ES and RA was varied in three directions: anterior-oblique direction, vertical direction, and posterioroblique direction. The angles between the anterior-oblique direction and vertical direction and between the posterioroblique direction and vertical direction were 15 degrees. Nine  conditions. The inferior endplate of L5 was rigidly fixed in six degrees of freedom. The intradiscal pressure (IDP), reaction moments, and intervertebral rotations (IVRs) were calculated. The FE program ABAQUS 6.13 (Dassault Systèmes, Versailles, France) was used for all simulations.

Validation.
The IVRs were recorded after applying moments of 7.5 Nm in flexion-extension, lateral bending, and axial rotation. The results for L3-L4 and L4-L5 were within the range of in vitro experimental data [18] (Figure 2).

Intradiscal Pressure (IDP).
The IDP recorded for L4-L5 was within 0.56-0.57 MPa for all nine loading conditions, which was similar to the 0.5 MPa in vivo pressure reported by Wilke et al. [32]. Therefore, a truly physiological IDP was achieved by the combination of the axial spinal load and upper body weight. Varying the direction of the muscle forces only had a minor influence on IDP for all nine loading conditions.

Reaction Moment.
Reaction moments acting against the boundary conditions which maintained the FE model in a neutral position were calculated and compared among the nine models (Figure 3). The three models (loading condition 7-9) with ES acting in an anterior-oblique direction showed smaller reaction moments. Overall, loading condition 9 produced the smallest reaction moment, while loading condition 1 produced the greatest. A greater reaction moment can be interpreted to mean that a greater external load is required to poise muscle forces and upper body weight.

Intervertebral Rotation.
In these loading conditions, the L3-L4 spinal functional unit showed flexion, while the L4-L5 showed extension. In loading condition 1-3, the combination

Discussion
The musculature plays an important role in maintaining the stability of the spine, but muscle forces are often difficult to measure due to the invasive methods required for in vivo measurements. Attempts have been made to incorporate muscle forces into FE models [9-15, 26, 33], but most of these forces were estimated based on locally developed algorithm/data. A generally practical method can be helpful for applications of muscle forces and driving private FE models. Muscle forces combined with follower loads and upper body weight investigated by Rohlmann et al. [26] may play an important role. A nonlinear FE model was employed in the present study in order to investigate how the direction of application of muscle forces impacts the motion of the lumbar spine. The magnitudes of all forces as well as the direction of the follower loads and upper body weight were kept constant. The physiological conditions are variable, so a sensitivity study was carried on. The action line of the ES and RA muscle forces was varied in three directions: anterioroblique direction, vertical direction, and posterior-oblique direction. In total, nine loading conditions were simulated through a combination of the different directions of ES and RA. Motion and loading data was recorded for further analysis and to determine the impact of muscle direction on the spinal motion segment.
The model was validated by confirming that the IVRs of the L3-L4 and L4-L5 spinal functional units were within the range of in vitro experimental data [18]. There were large variations for intervertebral rotations of the human body. Thus, the FE predictions in our study represent a part of populations instead of all. In addition, the calculated IDP of L4-L5 was 0.56-0.57 MPa, which was similar to the 0.5 MPa value reported in an in vivo study by Wilke et al. [32]. Thus, the spinal compressive loads in this study created by a combination of muscle forces and upper body weight are in line with physiological loads experienced under similar conditions, at least with quite a part of populations, since there are large variations in human bodies. Varying the direction of the muscle forces only had a minor influence on IDP when the spinal position is unchanged.
As the IDP and IVR are sensitive to posture [32], the parameters calculated for nine loading conditions cannot be compared in a unified way without restriction. Therefore, this study rigidly fixed L3 in a neutral standing position. The reaction moments to the boundary conditions reflect the degree of deviation to the neutral standing position. A greater reaction moment means a greater external load is needed to poise the muscle forces and upper body weight. The magnitudes of the reaction moments in loading conditions 7, 8, and 9 were relatively small ( Figure 3). ES acting in an anterioroblique direction needed a lower external force or moment to keep the spinal FE model in a neutral position, in comparison to when acting in a vertical direction or posterior-oblique direction. This finding was in agreement with anatomical observations and other biomechanical models [34]. The direction of RA also affected the reaction moment but had less of an influence than ES. The biomechanical environment of the abdomen is more complicated than shown in this model, because both of abdominal muscles and abdominal pressure take part in maintaining stability.
In physiological environment, upper body weight and muscle forces are well balanced. Due to the L3 being fixed in the FE model, the overall IVR was 0 degrees. However, the IVR for each spinal motion unit did vary with different loads. In Figures 4 and 5, the loading conditions with ES acting in an anterior-oblique direction showed smaller IVRs for L3-L4 and L4-L5. A satisfied equilibrium was achieved by the back muscle, abdominal muscle, follower load, and body weight. This is in accordance with the one from reaction moment. Similarly, the variations caused by varying the direction of RA were rather small.
There are some limitations in this study that should be noted. Firstly, the geometry and material properties were constant across all models. This may not be representative of the entire population. Secondly, only two spinal functional units were involved and only effect in sagittal direction was investigated. Thirdly, the effect of local muscles was simulated by a follower load of 200 N, which is a compressive load applied along a follower load path, thus minimizing the shear forces. A more detailed model including additional muscles would be beneficial in further studies. Although several simplifications were made, the results drawn from the comparison among the loading conditions are still valid since the direction of the muscle forces was changed in a parameterized way and the configurations were the same for all loading conditions. Thus, the results of this study can be viewed as a comparative analysis that provides information for applying muscle forces.

Conclusion
The spinal compressive loads delivered by muscle forces and upper body weight in this study are comparable to physiological conditions. This study also found that the muscle forces and upper body weight are capable of maintaining equilibrium in the case of erector spinae muscle in anterioroblique direction. The erector spinae muscle acting in the anterior-oblique direction needed a lower external force or moment to keep the spinal FE model in a neutral position, in comparison to the vertical direction and posterior-oblique direction.