Dependence of lumbar loads on spinopelvic sagittal alignment: An evaluation based on musculoskeletal modeling

Still little is known about how spinopelvic alignment affects spinal load distribution. Musculoskeletal modeling can potentially help to discover associations between spine alignment and risk factors of spinal disorders (e.g. disc herniation, vertebral fracture, spondylolisthesis, low back pain). The present study exploited the AnyBody full-body musculoskeletal model to assess the relation between lumbar loads and spinopelvic alignment in the sagittal plane. The model was evaluated in the standing position. The simulated postures were set using spinopelvic parameters gleaned from the literature and characterizing the healthy adult population. The parameters were: sagittal vertical axis, Roussouly lumbar type, sacral slope, and pelvic incidence. A total of 2772 configurations were simulated based on the following measurements: compression force and anterior shear at levels L4L5 and L5S1; multifidus, longissimus spinae, and rectus abdominis muscle forces. Changes in global sagittal alignment, lumbar typology, and sacral inclination, but not in pelvic incidence, were found to affect intervertebral loads in the lumbar spine and spinal muscle activation. Considering these changes would be advantageous for clinical evaluation, due to the recognized relation between altered loads and risk of disc herniation, vertebral fracture, spondylolisthesis, and low back pain. Musculoskeletal modeling proved to be a valuable biomechanical tool to non-invasively investigate the relation between internal loads and anatomical parameters.


Introduction
The human spine forms an S-shape in the sagittal plane, with convex curvature in the thoracic and sacral regions. The anatomical spinopelvic parameters, obtained by radiographic examination and clinically used to describe spine alignment, are related to four spine regions: cervical, thoracic, lumbar, and sacropelvic.
The most common parameter to assess global alignment is the sagittal vertical axis (SVA), which is defined as the horizontal offset from a plumb line dropped from the seventh cervical vertebra to the posterior-superior corner of the sacral endplate [1]. Being an index of sagittal imbalance, when the SVA magnitude exceeds the normal range, either toward the front or the back, the spine is considered malaligned. The lumbar spine is lordotic in shape (i.e., convex PLOS  simulations ( Fig 1A). By default, the body model represents the size and weight of an average European male (1.76 m, 75 kg). The model has been extensively validated for the assessment of lumbar loads and muscle activation during assigned postures and movements in physiological conditions [22][23][24]. Musculoskeletal modeling characterizes bones as rigid segments connected by joints, and muscles as tensile elements attached to segments and providing movements. By means of the inverse dynamic approach, muscle forces and intersegmental forces acting during the execution of specific imposed kinematics are computed by minimizing muscle recruitment activation [25,26]. The default polynomial optimization criterion was exploited. For the spine, the model characterizes the twelve thoracic vertebrae and ribcage as a single segment, the five lumbar vertebrae (from L1 to L5) as segments connected by spherical joints, and the sacrum and pelvis as rigidly connected segments. A total of 188 muscle fascicles and abdominal pressure are accounted for in the lumbar region.

Simulation process
The full-body model was evaluated in standing postures simulating spinopelvic alignment ( Fig  1) corresponding to changes in four anatomical parameters (SS, PI, RT, and SVA) in asymptomatic adults. Since RT implicitly accounts for changes in lumbar lordosis (LL), LL was not included as an additional parameter in this study. However, the corresponding LL values for the four different RTs are given in the results. The parameters values (Table 1) were obtained from the study by Hu et al. [27] who provided the mean values and the ranges of SS, PI, and SVA differentiated in the four RTs (Fig 2). Increments of 1˚were considered for SS and PI. Model posture was set by adjusting the anatomical parameters in the following sequence: SS, PI, RT, and SVA. By default, the local reference system for the sacrum segment is oriented parallel to the global reference system, but the SS is approximately 30˚ (Fig 3A). Accordingly, the changes in SS were simulated by rotating the sacrum segment to account for the SS default value. For example, the segment was rotated 10˚anteriorly to simulate a SS of 40˚. This movement also provided the corresponding rotations of the segment properties, such as joints and muscle insertion points, thus preserving the morphological characterization of the sacrum and the pelvis (the sacrum and pelvis being rigidly connected).
For the PI, changes in this angle were simulated by shifting, in the sagittal plane, the position of the hip joints (defined in the pelvis segment) that connect the pelvis to the femoral heads ( Fig 1D and Fig 3C). This procedure changed the attachment point between the femur and the pelvis, without altering the position of the muscle insertion points in the respective segment. After setting the SS and the PI, the four RTs (RT1, RT2, RT3, and RT4) were modeled to match the reference examples reported by Roussouly et al. [2]. Each RT was obtained by rotating the vertebral segments in the sagittal plane (Table 1). Proceeding from L5 to L1, each vertebral segment was rotated to shift the slope of the vertebral mesh from its original value in the default AnyBody model (Fig 1C) to that required in the specific RT (Fig 2). With this approach, morphological characterization (i.e., position of the intervertebral joints, vertebral centroid, and muscle insertion points) is kept unaltered for each vertebra (Fig 3A). Furthermore, the position of the center of mass of the vertebral segments was placed more anteriorly than the vertebral centroids, depending on the lumbar level, according to the computed tomography (CT) scan measurements reported by Pearsall et al. [28] (Fig 3A).
Finally, the three SVA values (dependent on RT, Table 1) corresponding to balanced alignment (SVAmed) and to backward and frontward imbalanced postures (SVAback and SVAfront) were modelled. In detail, the required SVA was achieved by setting the rotation of the thoracic segment with respect to the L1 vertebra ( Fig 1B). Neck rotation was set opposite to thoracic rotation in order to maintain horizontal gaze; arms were kept vertically aligned.
In all, 2772 spinopelvic configurations were simulated. The simulations were run in batch process using custom routines written in MATLAB (MathWorks Inc., Natick, MA, USA).

Model outputs
The following measurements were computed for each simulated configuration: intersegmental force at L4L5 (F L4L5 ) and at L5sacrum (F L5S1 ) joints; muscle forces of the multifidus (F MF ), longissimus spinae (F LS ), and rectus abdominis (F RA ), these muscles being involved in maintaining the spine aligned in upright posture. The axial compression and the anterior shear  components of F L4L5 (F c L4L5 and F s L4L5 ) were obtained by projecting F L4L5 on the axis passing through the upper and lower intervertebral joints of L5 (axial direction, caudally oriented) and on the orthogonal axis (anteriorly oriented, parallel to the upper endplate), respectively ( Fig  3B). The components of F L5S1 (F c L5S1 and F s L5S1 ) were obtained by projecting F L5S1 on the axis orthogonal to the direction of the default sacral inclination (axial direction, caudally oriented) and on the orthogonal axis (anteriorly oriented, parallel to the sacral endplate), respectively ( Fig 3B).

Intersegmental forces
The compression force at L4L5, F c L4L5 , in RT1 was larger in SVAfront alignment than in SVAmed and SVAback ( Fig 4A). A moderate negative linear relation was observed to depend on SS changes in SVAmed and SVAback ( Fig 4B). Changes in PI did not substantially affect the F c L4L5 values ( Fig 4C). This lack of relation with PI changes was, in general, evident in all the predicted compression and shear forces at both L4L5 and L5S1 (Figs 5 and 6). Generally, a balanced alignment posture (SVAmed) produced lower compression and shear forces (upper row in Figs 5 and 6). RT1 produced higher F c L4L5 , lower F s L4L5 , and higher F s L5S1 than the other RTs ( Table 2). Since the forces were mainly arranged as planar surfaces (Figs 5 and 6), which are intrinsically characterized by a non-Gaussian distribution, the descriptive values are reported as median and range in Tables 2 and 3. In RT1 the median F c L4L5 ranged from 543N to 790N, the median F s L4L5 from 27N to 35N, and the median F s L5S1 from 329N to 398N. The ranges for RT3 and RT4 were larger, with a median F s L4L5 from 57N to 76N and from 63N to 72N, respectively.
In order to assess the dependence of the intersegmental forces on the changes in the SS, we compared the force values computed at minimum and maximum SS (SS min and SS max ) at average PI (dependent on RT, see Table 1) (Figs 7 and 8).
At the L4L5 level, the compression force F c L4L5 was similar for SS min and SS max in all four RTs (upper row in Fig 7). The shear force F s L4L5 was lower in SS max in all four RTs (lower row in Fig 7). At the L5S1 level, the compression force F c L5S1 was moderately lower in SS max in all four RTs (upper row in Fig 8). The shear force F s L5S1 was larger in SS max in all four RTs (lower row in Fig 8).

Muscle forces
As noted for the intersegmental forces, we observed no relation between PI changes and muscle forces. The multifidus force, F MF , was lower in RT1 and RT2, and higher in SVAfront ( Fig  9, upper row). The median ranged from 16N to 41N in RT1 and from 17N to 24N in RT2 ( Table 3). The force values were higher in RT3 (from 25N to 42N) and RT4 (from 28N to 48N). The F MF was generally higher in SVAfront, and particularly in correspondence with SS max alignments (Fig 9, upper row).
The erector spinae force, F ES , was generally higher for SVAfront in all four RTs, with a mild positive linear relation with SS changes (Fig 9, central row). Overall, median values were lower in RT2 (range 25N to 194N, Table 3).
The rectus abdominis was activated only in the SVAback postures in all four RTs and mildly in SVAmed in RT4 (Fig 9, lower row). The muscle force, F RA , was increased from RT1 to RT4, and moderately lower in SS max . The median force values ranged from 26N to 175N, progressively increasing from RT1 to RT4 (Table 3).

Discussion
This section begins by discussing the results related to the anatomical parameters in the following order: sagittal vertical axis (SVA) interpreting global sagittal alignment, lumbar typology (RT), and spinopelvic parameters (SS and PI) (Figs 1 and 2). Although lumbar lordosis (LL) was not directly set in the simulated postures, its value was calculated (Figs 7-9). Since the orientations from L5 to L1 were defined as fixed in each RT (Table 1), in order to obtain an LL value dependent on postural change, the lordosis angle was computed between T12 and S1 ( Fig 1B). The inclination of T12 was defined as the slope of the vertebral mesh in the thoracic segment, and thus was related to SVA. The slope of S1 matched the SS. Further, LL was indirectly related to RT, as the changes in SS and SVA were dependent on lumbar typology (Table 1).
SVA-As expected, a balanced posture (SVAmed) produced lower compression forces in all four RTs than the backward or frontward imbalanced postures (SVAback and SVAfront) (Figs 5 and 6, upper row). The increase in the compression force at L4L5 and L5S1 at SVAfront was particularly pronounced in RT1 (Figs 7 and 8, upper rows). At L4L5, SVAfront also produced greater anterior shear in RT1 and RT2, whereas the shear was increased by SVAback in RT3 and RT4 (Fig 7, lower row). For muscle activation, SVAfront generally led to increased multifidus force, F MF , in all four RTs, whereas SVAback did not produce significant alterations (Fig  9, upper row). As expected, the erector spinae and rectus abdominis forces (F ES and F RA ) were mainly activated in the SVAfront and SVAback postures, respectively (Fig 9, central and lower  rows). Indeed, recruitment of the erector spinae (located dorsal to the spine) is necessary to counteract the shifting of trunk weight in postures imbalanced frontward and in balanced alignment as well, albeit with lower forces. Conversely, the rectus abdominis (located frontal in the lumbar region) is activated to counteract postures imbalanced backward. In particular, the rectus abdominis was not activated in either the SVAmed or the SVAfront posture (Table 3). These observations indicate that a frontward imbalanced posture is the one most able to increase loads on the lumbar spine. Our findings also underline the importance of achieving appropriate alignment in spinal reconstruction procedures, e.g., correcting thoracic hyperkyphosis, which is known to cause chronic pain and increase spinal load [29] or correcting hypo-lordosis in reconstruction surgery [30].
RT-Lumbar typology RT1 produced larger compression forces at L4L5 (Fig 7, upper row, and Table 2). Both RT1 and RT2 produced lower shear loads (Fig 7, lower row). Furthermore, RT1 was the typology most affected by SVAfront alignment in terms of alterations in compression force. Overall, lower shears were predicted in RT1 with corresponding large SS values (see SS max in lower row of Fig 7).
Conversely, RT1 produced larger shear forces at the L5S1 level. This finding can be explained by the result depicted in Fig 3B. Although the predicted intersegmental loads (F L5S1 and F L4L5 ) are similar in modulus and direction, the anterior component of F L5S1 (projection of the vector on the anterior axis) is larger than that of F L4L5 . The difference between the directions of the anterior axes of L5 and S1 is larger in RT1 (43˚at SS max ) than in the other RTs (range 23˚to 26˚) ( Table 1 and Fig 2). RT4 produced larger muscle forces in all the muscles, whereas RT2 produced lower F MF and F ES (Fig 7).
These results highlight the importance of considering lumbar typology when planning the treatment of disorders characterized by alterations in intervertebral lumbar loads. Axial load is associated with the risk of disk bulging and herniation [31][32][33] and the risk of vertebral fracture in patients with osteoporosis and reduced bone mineral density [34]. RT1 and RT4 warrant particular attention in osteoporosis and in other conditions at increased risk for disc herniation, such as ageing, obesity, and physically demanding work. Moreover, RT1 and RT2 generated lower activation of the multifidus muscle (Fig 9, upper row), the atrophy of which (producing lower forces) is known to be strongly associated with low back pain [35]. Generally, RT3 and RT4 produced increased anterior shear at L4L5 (and RT1 at L5S1), which can be a risk factor of anterior displacement in spondylolisthesis and spondylolysis [36][37][38].
SS and PI-While SS is defined as the slope of the sacral endplate in the sagittal plane (Fig 1), the AnyBody model accounts for bones as simple rigid segments (i.e., spherical or ellipsoidal masses), neglecting their surface shape and properties. As a consequence, the bony mesh surfaces (e.g., sacrum) displayed in the model view have to be considered for visual purposes only, without any mechanical meaning (Fig 3A). In the default standing position of the AnyBody model, the local reference system of the sacrum segment is oriented parallel to the global reference system, but SS visually results approximately 30˚ (Fig 3A). This orientation was originally set according to anatomic studies [23,24,39]. Because the segment properties (i.e., center of mass, joints, and muscle insertion points) are defined in the local reference system, we simulated the SS changes (along with changes in the corresponding segment properties) by rotating the sacrum segment to account for the SS default value. Greater SS produced moderately higher compression forces at L4L5 in RT2 and RT3 (Fig 7,  upper row) and more consistently affected anterior shear. Indeed, lower and higher shear forces at SS max were noted at L4L5 and L5S1, respectively, in all four RTs (Figs 7 and 8, lower , upper plot) and anterior shear (F s L4L5 , lower plot) at the L4L5 level, computed at minimum and maximum sacral slope (SS min and SS max ) for the four lumbar typologies (RT1, RT2, RT3, RT4). SS min and SS max were, respectively, 25å nd 35˚for both RT1 and RT2, 35˚and 45˚for RT3, and 45˚and 55˚for RT4 ( Table 1). The results correspond to the central pelvic incidence (PI) values for the RT (40˚in RT1 and RT2, 50˚in RT3, 60˚in RT4, see Table 1). The three sagittal vertical axis (SVA) conditions interpreting the balanced global alignment (SVAmed) and the backward and frontward imbalanced alignments (SVAback and SVAfront) are presented as green, blue, and red bars, respectively. The lumbar lordosis (LL) is reported as well.
https://doi.org/10.1371/journal.pone.0207997.g007 rows). In RT3 and RT4, which are the most common lumbar types in healthy adults [2], postures with greater pelvic retroversion (corresponding to SS min values) produced greater shear loads at L5S1. Accordingly, the degree of sacral inclination should be carefully considered in each RT and in conditions with risk factors of anterior displacement (e.g., spondylolisthesis and spondylolysis). As concerns muscle forces, SS max produced greater F MF in RT4 and moderately greater F ES in RT1, RT2, and RT3 (Fig 9). These results should be interpreted with caution, however, since there is demonstrated evidence for an association between multifidus atrophy and low back pain, while the evidence regarding the erector spinae is conflicting [40].
Surprisingly, changes in PI did not affect compression force at L4L5 (Fig 4C) or at L5S1. This lack of relation with PI was also generally predicted for anterior shear and muscle forces. It is commonly reported and accepted that PI is correlated with the degree of lumbar lordosis (LL). In other words, subjects with a large PI will have a large LL, while those with a small PI , upper plot) and anterior shear (F s L5S1 , lower plot) at the L5S1 level computed at the minimum and maximum sacral slope values (SS min and SS max ) in the four lumbar typologies (RT1, RT2, RT3, RT4). SS min and SS max are, respectively, 25˚and 35˚for both RT1 and RT2, 35˚and 45˚for RT3, 45˚and 55˚for RT4 ( Table 1). The results correspond to the central pelvic incidence (PI) values for the RT (40˚in RT1 and RT2, 50˚in RT3, 60˚in RT4, see Table 1). The three sagittal vertical axis (SVA) conditions interpreting the balanced global alignment (SVAmed) and the backward and frontward imbalanced alignments (SVAback and SVAfront) are presented as green, blue, and red bars, respectively. The lumbar lordosis (LL) is reported as well.
https://doi.org/10.1371/journal.pone.0207997.g008 SS min and SS max are, respectively, 25˚and 35˚for both RT1 and RT2, 35˚and 45˚for RT3, and 45˚and 55˚for RT4 ( Table 1). The results are taken in correspondence to the central pelvic incidence (PI) values for each RT (40˚in RT1 and RT2, 50˚in RT3, 60˚in RT4, Table 1). The three sagittal vertical axis (SVA) conditions interpreting the balanced global alignment (SVAmed) and the backward and frontward imbalanced alignment (SVAback and SVAfront), are presented as green, blue and red, bars, respectively. The lumbar lordosis (LL) is reported as well.
https://doi.org/10.1371/journal.pone.0207997.g009 have a small LL. A mismatch between PI and LL of less than 10˚is targeted in spinal corrective surgery to attain satisfactory spinopelvic alignment [13]. A larger LL is known to cause changes in lumbar loads [7,41,42], whereas PI changes were not found to affect the load distribution in the present study. This finding suggests that although PI and LL are correlated, only the latter parameter can be considered directly responsible for altering spinal load. For example, in the balanced posture (SVAmed), moderately greater compression forces at L4L5 were obtained in RT2, RT3, and RT4 with increased LL (Fig 7, 10˚difference between SS min and SS max ). Furthermore, the PI-LL mismatch has been found as an important parameter associated to increased load in the lumbar spine, specifically in relation to adjacent segment degeneration following fusion surgeries [20,30,43]. Unfortunately, the direct comparison with those studies is not feasible. Indeed, they exploited parameters from pathological subjects or from in vitro tests, while the present study has taken into account physiological ranges. Moreover, in the present study the calculation of LL (based on the slope of T12) is potentially affected by assuming thorax as single segment. Future investigations exploiting real patient data and specifically setting LL could elucidate the effects of PI-LL, but also distinguish between the individual contributions of LL and PI.
The present study has several limitations. The AnyBody full-body model represents only the size and weight of an average European male. The stiffness of the intervertebral discs was neglected. The facet joints are not modelled in the default AnyBody model. Nevertheless, they are not expected to have an impact on the loads at the motion segment in the upright position, but rather on the load sharing among the different structures in each motion segment, which is beyond the scope of the present work. The thoracic segment, which was used to obtain the position of C7 and to set SVA (Fig 1), characterizes the twelve thoracic vertebrae and ribcage as a single segment, without allowing the differentiation of further vertebral arrangements. This limitation can cause potential artefacts in the loads at the thoracolumbar level (T12L1), which should be assessed in future studies. The imbalanced postures are modelled by rotating the thorax with respect to L1, instead of being distributed along the different thoracic levels as would be expected physiologically. The present study focused specifically on investigating loads at lower lumbar levels; however, the lumped modelling assumptions for the thorax could have implications at T12L1, which could be carried over and consequently affect loading predictions in the lower lumbar joints as well. Nonetheless, the rigid thorax assumption was demonstrated to affect loading predictions at the upper but not at the lower lumbar levels [44].
Another limitation is the effects of the rotation between the thoracic and the lumbar region (T12L1) and between the lumbar region and the sacrum (L5S1) to obtain the spinopelvic configurations. Indeed, while the lumbar typology (RT) is defined by fixed orientations from L5 to L1 (Table 1), the inclinations of the thorax and the sacrum (in the same RT) varied according to the SVA and SS (obtained by rotating T12L1 and L5S1 joints, respectively). This modelling strategy was necessary, given the absence of subject-specific anatomical data reporting the distribution of rotations between the lumbar and the adjacent spine levels for the specific RT. In principle, however, it could produce large wedge angles at T12L1 and L5S1. In the simulations the extension angle between L5 and S1 ranged from -43˚to -13˚, which is comparable with that reported in asymptomatic subjects (range -35˚to -10˚) [45]. Concerning T12L1, since the AnyBody model characterizes the twelve thoracic vertebrae as a single segment, the relative angle between T12 and L1 was obtained by measuring the slope of the mesh of T12 in the thoracic segment (Fig 1B). The angle ranged from -23˚to +11˚, whereas in asymptomatic subjects it ranged from -7˚to +17˚ [44]. The larger value in extension (-23˚), although moderate, can affect the prediction of loads at T12L1. This aspect needs to be better clarified in future investigations, although it is related to the default choice for obtaining orientation of T12. Indeed, as mentioned above, the default AnyBody model does not allow to distribute thoracic rotation along different vertebral levels and constrains the user to orientate the mesh of T12 in the thoracic segment.
A further observation is the dependence of the spinopelvic parameters on each other. The sagittal parameters were varied according to the ranges obtained from the literature and independently one from the other. However, an inherited interdependence is expected in reality, e.g., frontally imbalanced (SVAfront) individuals may not present the lowest SS in a specific RT. Since the values published in the literature describe the distribution in the population, they do not allow to typify the interdependence between the parameters of single individuals. The simulations we evaluated explore all possible configurations among the parameters, including within those typifying interdependency. Nonetheless, real subject-specific data (e.g., from radiographic or CT images) are needed in future model developments to focus on alignments with interdependent parameters.
In conclusion, our findings indicate that changes in global sagittal alignment, lumbar typology, and sacral inclination, but not in pelvic incidence, can affect intervertebral loads in the lumbar spine and spinal muscle activation. Accounting for these variations would be advantageous for clinical evaluation, owing to the relation between altered loads and the risk of disc herniation, vertebral fracture, anterior displacement, and low back pain as well. Musculoskeletal modeling was found to be a valuable biomechanical tool to non-invasively investigate the relation between internal loads and spinopelvic parameters. In order to broaden the extent of the results, future developments will need to assess the relation between loads and anatomical parameters in other poses (e.g., trunk flexion-extension and bending) and dynamic conditions.