Effect of Passive Support of the Spinal Muscles on the Biomechanics of a Lumbar Finite Element Model

Finite element (FE) modeling of the passive ligamentous spine is widely used to assess various biomechanical behaviors. Currently, FE models that incorporate the vertebrae, ligaments, and the personalized geometry of the bony spine may be used in conjunction with external loads from the muscles. However, while the muscles place a load (moment) on the spine and support it simultaneously, the effect of the passive support from the adjacent spinal muscles has not been considered. This study thus aims to investigate the effect of passive support from the psoas major, quadratus lumborum, and erector muscles on the range of motion (RoM) and intradiscal pressure (IDP) of the lumbar spine. Various L2-sacrum spinal models that differed only in their muscle properties were constructed and loaded with a pure moment (2.5–15.0 Nm) alone or combined with a compressive (440 or 1000 N) follower load. The RoM and IDP of the model that excluded the effect of muscles closely matched previous FE results under the corresponding load conditions. When the muscles (40–160 kPa) were included in the FE model, the RoM at L2 was reduced by up to 6.57% under a pure moment (10 Nm). The IDP was reduced by up to 6.45% under flexion and 6.84% under extension. It was also found that the erector muscles had a greater effect than the psoas major and quadratus muscles.


Introduction
Lower back pain is a major medical issue that affects about 80% of all adults worldwide [1] and is the second most common symptom for patients visiting a physician for a general illness [2]. Due to the high prevalence, there have been many studies on the characteristics of the spine, demonstrating that determining the physiological load on the lumbar spine is crucial for understanding lumbar spine disease and for effective lumbar spine surgery and injury prevention [3,4].
Computational modeling of the passive ligamentous spine is widely used to assess various biomechanical behaviors [5][6][7][8][9][10]. Currently, finite element (FE) modeling of the spine is frequently used to assess deformation and stress distribution within the joint structure of the spinal system [11][12][13][14][15]. To establish valid spinal models, accurate mechanical properties are required for each component and their representations.
Recently, studies have investigated the effects of ligament properties on FE model responses. For example, when the strength of the ligaments varied considerably within the model, momentary sharing and internal forces within the disc and joints were affected [21][22][23][24]. Discs have also been modeled with both linear [25,26] and nonlinear properties [27][28][29][30]. However, the modeling of the muscles has often been ignored, being replaced by corresponding external loads on the spinal bone. In a recent study, more realistic muscle load conditions were assessed using a combined musculoskeletal and FE model of the lumbosacral spine [31]. However, the effect of the passive support from adjacent spinal muscles was not considered.
This study thus aims to investigate the effect of passive support from the psoas major, quadratus lumborum, and erector muscles on the range of motion (RoM) and intradiscal pressure (IDP) of the lumbar spine. We also investigate the effects of the stiffness of each component. We constructed a lumbar spine model that incorporated these muscle groups by adding computed tomography (CT)-based muscle elements to a conventional passive ligamentous spine model verified using data from the literature [32]. The effects of the passive support from the adjacent spinal muscles on the RoM and IDP of the lumbar spine were tested under various load conditions by comparing the results of the models with and without muscle parts.

FE Model without Muscles
The geometry of the L2-sacrum vertebrae was obtained from the lumbar spine of a 34-year-old woman reconstructed from CT scans (thickness = 2 mm). The model was reconstructed using Mimics (Materialise Technologies, Belgium) and a mesh was generated with Hypermesh (Altair Hyper-works 14.0, Troy, MI, USA). Figure 1 shows the spinal model for the L2-sacrum vertebrae. The lumbar bones consisted of a cortical bone and a cancellous bone with 4-node tetrahedral elements. The thickness of the cortical bone was 1 mm, and the disc was divided into the nucleus pulposus at the center and the annulus fibrosus surrounding it. Since the disc nucleus and annulus could not be clearly distinguished in the spinal CT images, the cross-section of the disc nucleus was manually constructed, to have approximately 30-50% of the space occupied by the annulus [33]. Seven ligaments (ALL, PLL, LF, SSL, ISL, ITL, and CL) were included in the model. Table 1 shows the number of elements for each component of the model. Ligaments were modeled as elastic bodies with volume. All interactions between different bony parts were considered as frictionless contact while the ligaments were attached to the bone with fixed conditions. Appl. Sci. 2020, 10, x FOR PEER REVIEW  2 of 11 Recently, studies have investigated the effects of ligament properties on FE model responses. For example, when the strength of the ligaments varied considerably within the model, momentary sharing and internal forces within the disc and joints were affected [21][22][23][24]. Discs have also been modeled with both linear [25,26] and nonlinear properties [27][28][29][30]. However, the modeling of the muscles has often been ignored, being replaced by corresponding external loads on the spinal bone. In a recent study, more realistic muscle load conditions were assessed using a combined musculoskeletal and FE model of the lumbosacral spine [31]. However, the effect of the passive support from adjacent spinal muscles was not considered.
This study thus aims to investigate the effect of passive support from the psoas major, quadratus lumborum, and erector muscles on the range of motion (RoM) and intradiscal pressure (IDP) of the lumbar spine. We also investigate the effects of the stiffness of each component. We constructed a lumbar spine model that incorporated these muscle groups by adding computed tomography (CT)based muscle elements to a conventional passive ligamentous spine model verified using data from the literature [32]. The effects of the passive support from the adjacent spinal muscles on the RoM and IDP of the lumbar spine were tested under various load conditions by comparing the results of the models with and without muscle parts.

FE Model without Muscles
The geometry of the L2-sacrum vertebrae was obtained from the lumbar spine of a 34-year-old woman reconstructed from CT scans (thickness = 2 mm). The model was reconstructed using Mimics (Materialise Technologies, Belgium) and a mesh was generated with Hypermesh (Altair Hyperworks 14.0, Troy, MI, USA). Figure 1 shows the spinal model for the L2-sacrum vertebrae. The lumbar bones consisted of a cortical bone and a cancellous bone with 4-node tetrahedral elements. The thickness of the cortical bone was 1 mm, and the disc was divided into the nucleus pulposus at the center and the annulus fibrosus surrounding it. Since the disc nucleus and annulus could not be clearly distinguished in the spinal CT images, the cross-section of the disc nucleus was manually constructed, to have approximately 30-50% of the space occupied by the annulus [33]. Seven ligaments (ALL, PLL, LF, SSL, ISL, ITL, and CL) were included in the model. Table 1 shows the number of elements for each component of the model. Ligaments were modeled as elastic bodies with volume. All interactions between different bony parts were considered as frictionless contact while the ligaments were attached to the bone with fixed conditions.   The density of the mesh for the soft tissue sections expected to be more deformed was increased and optimized. When the density of the constructed mesh was increased 1.5-fold, the difference in the maximum value of the stress of the model was less than 5%. The constructed mesh was selected for the analysis as shown in Table 1. The material properties are listed in Table 2. For all materials, a linear elastic behavior was considered.

Muscle Components
We divided the muscles associated with the lumbar vertebrae into two groups based on the transverse process: (1) The psoas and quadratus muscles, and (2) the erector muscle adjacent to the lumbar spinous process. To investigate the effects of each muscle group, we constructed four versions of the model: (1) Bone without muscle, (2) bone with psoas major and quadratus lumborum muscles (P&Q), (3) bone with erector spinae muscle (E), and (4) bone with psoas major, quadratus lumborum, and erector spinae muscles models (P&Q+E) (Figure 2). Since the intensity of the muscles in the CT images was similar to that of organs and blood vessels, the muscles were manually selected from the CT data and then stacked to create the models. All muscles were meshed using 3-node trials at 1 mm. The volume meshes were meshed elements with 4-node tetrahedral solid elements. Psoas major and quadratus muscles were defined by one volume. This simplified muscle volume is adjacent to the lumbar spine, and the erector muscle is adjacent to the lumbar spinous process. Each of these adjacent muscles was anchored to the lumbar bone, and all the interactions between different parts were modeled as frictionless contact conditions.
To evaluate the influence of the stiffness parameters for each component, Young's modulus in the range of 40-160 kPa was evaluated [29,35]. Poisson's ratio was set at 0.49 for all muscles to represent the quasi-incompressibility of muscle tissue [29,35].
Appl. Sci. 2020, 10, x FOR PEER REVIEW 4 of 11 muscles was anchored to the lumbar bone, and all the interactions between different parts were modeled as frictionless contact conditions. To evaluate the influence of the stiffness parameters for each component, Young's modulus in the range of 40-160 kPa was evaluated [29,35]. Poisson's ratio was set at 0.49 for all muscles to represent the quasi-incompressibility of muscle tissue [29,35].

Simulated Tasks
The actual behavior of the human body begins with the force generated by the muscles. However, the application of the movement of these muscles to analysis needs to be simplified because the conditions are difficult to establish. Therefore, the active force produced by the muscles analyzed in this study was expressed as the loading conditions. To evaluate the effect of the muscles on the RoM, we considered various loading conditions used in previous works [36]. We first considered six pure moments of flexion and extension (2.5, 5.0, 7.5, 10.0, 12.5, and 15.0 N·m) applied to the top of L2 while the sacrum was fixed. The moment was also combined with a preload to evaluate the effects on IDP in the nucleus of the disc. For the combined loads with a moment of 15 N·m, a follower preload (1000 N for flexion and 440 N for extension) was used based on a previous study [23]. All analyses were conducted in Abaqus (Dassault Systemes Simulia Corp., Johnston, RI, USA).

Simulated Tasks
The actual behavior of the human body begins with the force generated by the muscles. However, the application of the movement of these muscles to analysis needs to be simplified because the conditions are difficult to establish. Therefore, the active force produced by the muscles analyzed in this study was expressed as the loading conditions. To evaluate the effect of the muscles on the RoM, we considered various loading conditions used in previous works [36]. We first considered six pure moments of flexion and extension (2.5, 5.0, 7.5, 10.0, 12.5, and 15.0 N·m) applied to the top of L2 while the sacrum was fixed. The moment was also combined with a preload to evaluate the effects on IDP in the nucleus of the disc. For the combined loads with a moment of 15 N·m, a follower preload (1000 N for flexion and 440 N for extension) was used based on a previous study [23]. All analyses were conducted in Abaqus (Dassault Systemes Simulia Corp., Johnston, RI, USA).

RoM
The RoM for the L2 top plate models was compared in relation to the moment ( Table 3). The models produced a wide RoM in terms of flexion and extension depending on the muscle properties. In all cases, the P&Q + E model was the most restricted in terms of the two movements. At a muscle modulus of 160 kPa and a moment of 10 Nm, the motion modulus for the P&Q + E model was 6.57% lower than that of the muscle-free model.

IDP under Combined Loading
The IDP of various models under a combined loading were compared ( Table 4). The predicted IDPs were highest in the model for the bone without muscle at all levels, whereas the P&Q + E muscle model predicted the lowest IDPs. For the nucleus 5S muscle at 160 kPa, the P&Q + E and muscle-free models had an IDP difference of about 109 kPa (6.45%). There was a difference of about 24.1 kPa (3.49%) at nucleus 23, 36.7 kPa (4.43%) at nucleus 34, and 56.7 kPa (5.20%) at nucleus 45.
The IDPs predicted by the proposed models under extension were also compared, with the highest values observed in the model for the bone without muscle and the lowest observed in the P&Q + E model (Figure 3). At 160 kPa, the P&Q + E and muscle-free models differed by about 1.52 kPa (2.21%) at nucleus 23, 2.98 kPa (5.31%) at nucleus 34, and 6.49 kPa (25.5%) at nucleus 45. However, the IDP at nucleus 5S exhibited an opposite trend to that of the other nuclei.

Model Validation
Two FE models without muscles were tested and compared with the results reported by Jiang et al. [25]. For validation, the same components and material properties were considered for comparison. RoM was compared using a model consisting of only the L4-sacrum, disc, and ligaments under a preload of 150 N combined with 7.5 Nm on the top of the L4 planes, whereas the nucleus IDP was compared using a model of L3-L4 under a preload of 500 N combined with 7.5 Nm on the top of the L3 plane. Flexion, extension, lateral bending, and axial rotation were considered as boundary conditions. In the Figure 4a, the RoMs for the two versions of the FE model without muscles were compared to previously reported numerical data [37][38][39][40][41]. Under flexion, a maximum error of about 1.9 • was calculated for all sections, compared to 2.4 • under extension, and the average error for all cases was about 0.6 • . The global spine motion exhibited a similar trend.
In the Figure 4b, the IDP for L3-4 in terms of flexion, extension, lateral bending, and axial rotation was compared to data reported by El-Rich et al. [42] and Schmidt et al. [43]. For flexion, lateral bending, and axial rotation, the proposed model predicted slightly lower IDPs compared to the existing models and predicted significantly lower IDPs for extension. Jiang et al. [25] reported that the IDP for vertebral motion can be affected by the inherent sacral slope in each individual. Thus, though the IDPs predicted for extension were not similar, the difference was most likely due to differences in the shape of the model: The current model predicted 3-fold less IDP for extension.
rotation was compared to data reported by El-Rich et al. [42] and Schmidt et al. [43]. For flexion, lateral bending, and axial rotation, the proposed model predicted slightly lower IDPs compared to the existing models and predicted significantly lower IDPs for extension. Jiang et al. [25] reported that the IDP for vertebral motion can be affected by the inherent sacral slope in each individual. Thus, though the IDPs predicted for extension were not similar, the difference was most likely due to differences in the shape of the model: The current model predicted 3-fold less IDP for extension.

Effect of Muscle Volume on RoM
The proposed FE model accurately predicted spinal force, and the incorporation of the surrounding muscles had no effects on spinal behavior. The version of the model with all muscles included had a RoM difference of 0.4° at most compared to the model without muscle.

Effect of Muscle Volume on RoM
The proposed FE model accurately predicted spinal force, and the incorporation of the surrounding muscles had no effects on spinal behavior. The version of the model with all muscles included had a RoM difference of 0.4 • at most compared to the model without muscle.
It was found that the higher the Young's modulus of the muscle was, the lower the vertebral RoM became. Additionally, the higher the stiffness of the muscle was, the stronger the resisting force of the vertebrae became, leading to increased stiffness in the entire spine. In addition, when comparing the effects of each muscle, it was found that the erector muscle always had a stronger effect on the vertebral RoM for the various moments. This is because the erector muscle is larger than the other two muscle groups and muscles on the side adjacent to the spinous process have a greater influence on the spinal transverse process.

Effect of Muscle Volume on IDP
The IDP of the muscle model was relatively low in most cases, with a maximum difference of 7%. The muscles are responsible for supporting the spine, thus dispersing the pressure on the spinal disc. As the modulus of the spine increased, the IDP decreased ( Table 4). The overall pressure distribution in the vertebrae revealed that the largest load was on the bottom disc, which is most affected by the load of the upper body near the caudal bone.
Under extension, nucleus 5S predicted IDP trends that differed from all of the other nuclei. This may be because the pressure is lower at nucleus 5S, or it may be an inherent characteristic of the model used. In future research, it is necessary to further analyze this observation by comparing the vertebrae of various case studies using FE analysis.

Limitations and Future Works
FEM analysis studies involving muscles and ligaments are lacking compared to bones. It is more accurate to assume the muscle is a hyper-elastic body while for the limited range of motions, the linear material assumption also provides the reasonable accuracy [30]; therefore, it is frequently used in previous works. However, the behavior of nonlinear models is largely different from that of linear models when the deformation of the model is large. Therefore, the results of the current study would not be directly applicable for the model with large motions.
In this study, the calculation was also simplified with the assumption that the contact between different parts were correctly attached under frictionless contact conditions. For a more realistic interpretation, we should consider more realistic contact conditions with the experimentally obtained friction coefficients. Moreover, the effects of the intra-abdominal pressure or volume of the muscles, other surrounding tissues were not investigated.
Furthermore, the assessment of the passive support process was performed without taking into account the line of action based on the point of contact, as the direction of the supporting force would vary depending on the shape in which the muscle was compressed. In addition, applying anisotropic properties with different physical properties according to the direction of the muscle will result in a more precise model, which may better express the direction of the passive bearing force. The inclusion of the line of action and anisotropic properties in the passive support of muscles will be an essential part of future research. Moreover, there is a difference in the volume of muscle in the supine and standing positions, and there may be an effect, and the differences in the size of the muscles are more prominent from person to person.
This study only considered the effects of additional stiffness from the muscles, not the mass of the muscles. Research on the effects of the muscles in the spinal dynamic analysis is needed. Lastly, the future studies with experimental validations would be of great value.

Conclusions
In this study, a passive muscle model was constructed for the analysis of in vivo vertebral biomechanics. The results showed that the incorporation of muscle resulted in no more than a 7% difference in vertebral biomechanics. Thus, the RoM and IDP calculated from models without muscles, which are frequently used, may be at most 7% overestimated. The 7% differences among the models may not be significant as there are from 10-20% differences in the results of the previous results: The effects of the different geometry and the analysis conditions have more prominent effects on the results of the IDP. The amount of additional support from the muscles may not be significant and the effect throughout the model was uniform. Therefore, the passive muscles model should only be needed when a very accurate spinal musculoskeletal mechanics is required.