Sensitivity of muscle and intervertebral disc force computations to variations in muscle attachment sites

Abstract

The current paper aims at assessing the sensitivity of muscle and intervertebral disc force computations against potential errors in modeling muscle attachment sites. We perturbed each attachment location in a complete and coherent musculoskeletal model of the human spine and quantified the changes in muscle and disc forces during standing upright, flexion, lateral bending, and axial rotation of the trunk. Although the majority of the muscles caused minor changes (less than 5%) in the disc forces, certain muscle groups, for example, quadratus lumborum, altered the shear and compressive forces as high as 353% and 17%, respectively. Furthermore, percent changes were higher in the shear forces than in the compressive forces. Our analyses identified certain muscles in the rib cage (intercostales interni and intercostales externi) and lumbar spine (quadratus lumborum and longissimus thoracis) as being more influential for computing muscle and disc forces. Furthermore, the disc forces at the L4/L5 joint were the most sensitive against muscle attachment sites, followed by T6/T7 and T12/L1 joints. Presented findings suggest that modeling muscle attachment sites based on solely anatomical illustrations might lead to erroneous evaluation of internal forces and promote using anatomical datasets where these locations were accurately measured. When developing a personalized model of the spine, certain care should also be paid especially for the muscles indicated in this work.

Keywords:

• Muscle force; musculoskeletal model; spine loads; sensitivity; muscle attachment

1. Introduction

Computational models of the musculoskeletal system provide great opportunities to improve existing diagnosis and treatment practices for musculoskeletal and neuromuscular disorders. In the spine, detailed assessments of internal loads in the spinal tissues and muscle coordination strategies during trunk movement can help to achieve these clinical objectives (Hug and Tucker 2017). For example, Bruno et al. (2017) elaborated on the higher rates of vertebral fractures in the middle thoracic and thoracolumbar regions by computing vertebral compressive forces and the associated risk factors for a series of daily activities. The credibility of such biomechanical analyses depends largely on how well the models represent the true skeletal and muscular anatomy in all its complexity. Since muscle lines-of-action determine torque contributions at the joints, modeling muscle attachment sites at the bones–that is muscle origin, insertion, and via points–requires special attention for accurate assessment of joint and muscle forces.

In the past years, attachments of major spinal muscles were reported mostly through anatomical drawings that were created in detailed dissection studies on the lumbar spine (Bogduk, Macintosh, et al. 1992; Bogduk, Pearcy, et al. 1992; Bogduk et al. 1998; Macintosh and Bogduk 1991; Macintosh et al. 1986; Delp et al. 2001; Phillips et al. 2008; Brown et al. 2011) and cervical spine (Kamibayashi and Richmond 1998; Borst et al. 2011). These illustrations were prepared by visually inspecting the tendinous attachments of muscle bundles at the bones and estimating average locations across the cadavers studied rather than measuring the positions. Subsequently, the majority of generic musculoskeletal models adopted these illustrations for defining muscle lines-of-action (De Zee et al. 2007; Gagnon et al. 2011; Christophy et al. 2012; Ignasiak, Dendorfer, et al. 2016). Other studies introduced novel techniques to adjust muscle paths and cross-sectional areas of major trunk muscles based on two-dimensional medical images (Anderson et al. 2012; Bruno et al. 2015; Eskandari et al. 2017). Recently, we completed an entire anatomical dataset measured from one single human spine (Bayoglu et al. 2017a, 2017b). In this dataset, the coordinates of muscle attachment sites (quantified by optical tracking technology), spinal geometry, and comprehensive morphological parameters of 321 muscle-tendon elements were reported. Subsequently, a complete and coherent musculoskeletal model for the entire human spine (the Twente Spine Model) was developed based on this dataset (Bayoglu et al. 2019). With musculoskeletal models in general, and thus with ours, it is an important question that how sensitive the outcomes are to errors in the input data.

Previously, several authors investigated the sensitivity of their models against the locations of intervertebral center of rotation (Zander et al. 2016), spine rhythms and abdominal pressure (Arshad et al. 2016), ligament stiffness (Putzer, Auer, et al. 2016), and vertebral morphology (Putzer, Ehrlich, et al. 2016). Muscle recruitment during trunk exertions is dependent on muscle lines-of-action, however, there is only one study which aimed to investigate the sensitivity of muscle and intervertebral disc force computations against muscle attachment sites in the spine. Nussbaum et al. (1995) varied lines-of-action of some lumbar muscles within anatomically feasible ranges to reflect inter-individual differences and studied the effects on predicted muscle and disc forces. Unfortunately, their model included only a few muscle groups, and the findings were limited to a single disc level in the lumbar spine. Thus, the aim of the present study is to assess the effects of variations in muscle attachment sites on the computed muscle and intervertebral disc forces in a musculoskeletal model of the spine. We focused on potential measurement errors in muscle origin, insertion, and via point locations and subsequent implementation in musculoskeletal models. Specifically, each attachment location was systematically perturbed, and the global reaction of the muscular system to the corresponding perturbation was calculated. We also analyzed the changes in compressive and shear forces at the T6/T7, T12/L1, and L4/L5 discs.

2. Materials and methods

The Twente Spine Model was used to assess the sensitivity of muscle and intervertebral disc joint reaction force computations to muscle attachment locations (Bayoglu et al. 2019). This model was based on a complete and coherent anatomical dataset measured from one male cadaver (height: 154 cm, mass: 51 kg) (Bayoglu et al. 2017a, 2017b) and was made right to left symmetrical (see Figure 1). It included 60 bony segments (cervical, thoracic, and lumbar vertebrae, ribs, skull, sternum, hyoid, thyrohyoid, clavicles, scapulas, humeri, sacrum, and pelvis) and 351 muscle-tendon force actuators (per body side, see Table 1) and had 193 degrees-of-freedom (DoF). Intervertebral articulations were modeled as spherical, costo-vertebral and costo-transverse articulations as compound revolute, and the costo-sternal articulations as six DoF joints. Further, the lines-of-action of each muscle-tendon actuator (muscle element) was constructed by its origin and insertion where the force is transferred along its path. We used additional via points for muscle elements with a curved lines-of-action to represent their path better. At via points, muscle force is transferred in the direction of a line that bisects the angle formed by the muscle on either side of the point. The isometric strength of each muscle element was calculated by multiplying its physiological-cross sectional area (PCSA) by the specific muscle tension (100 N/cm², (Hansen et al. 2006)). For simplicity, muscle force-length and force-velocity relationships were not considered. The model relied on inverse dynamics and employs a static optimization criterion (cubic polynomial) for computing muscle forces (Rasmussen et al. 2001).

Figure 1. (a) The musculoskeletal model used in this study (Bayoglu et al. 2019), bones are displayed in beige and muscle-tendon elements in red. The reference frame depicted in the sagittal view is the global reference frame (${\psi }_{\text{global}}$). In the global frame, +X^g, +Y^g, and +Z^g-axes point anteriorly, cranially, and laterally (to the right side of the cadaver), respectively. The vertebral (${\psi }_{\text{vertebra}}$) and intervertebral disc joint reference frames (${\psi }_{\text{joint}}$) are shown for one motion unit. The orientations of the joint reference frames were defined as the average of their corresponding superior and inferior vertebrae reference frames (Senteler et al. 2016). The compressive force at the intervertebral disc is calculated along the y^j-axis of the joint reference frame while the shear is the force component on the xz plane. (b) Here the perturbations of the muscle attachment locations are illustrated. In this example, the measured location (nominal) was perturbed 2.5 mm along the negative Z^g-axis. The modified muscle path due to this perturbation is indicated by the red dashed line.

Table 1. Muscles included in the Twente Spine Model.

Region	Muscle
Cervical spine	Iliocostalis cervicis, Levator scapulae, Longissimus capitis, Longissimus cervicis Longus capitis, Obliquus capitis inferior, Obliquus capitis superior, Omohyoid Rectus capitis anterior, Rectus capitis lateralis, Rectus capitis posterior major Rectus capitis posterior minor, Scalenus anterior, Scalenus medius, Scalenus posterior Semispinalis capitis, Semispinalis cervicis, Splenius capitis, Splenius cervicis Sternocleidomastoid, Sternohyoid, Sternothyroid, Thyrohyoid
Perturbed muscles
Thoracic spine	Iliocostalis thoracis (IT), Rotatores (ROT), Semispinalis thoracis (SST) Serratus posterior superior (SPS), Spinalis thoracis (ST), Trapezius (TP) Rhomboideus major (RMMA), Rhomboideus minor (RMMI), Serratus anterior (SA) Subclavius (SCL)
Rib cage	Intercostales externi (IE), Intercostales interni (II), Levatores costarum (LEC) Subcostales (SCO), Transversus thoracis (TT)
Lumbar spine	Iliocostalis lumborum (IL), Longissimus thoracis (LT), Multifidus (MF) Obliquus externus abdominis (OEA), Obliquus internus abdominis (OIA), Psoas major (PM) Quadratus lumborum (QL), Rectus abdominis (RA), Serratus posterior inferior (SPI) Latissimus dorsi (LD)

In this study, only the muscles located in the thoracic spine, rib cage, and lumbar spine were perturbed. Abbreviations used for the muscle names are indicated inside the parenthesizes.

The Twente Spine Model’s estimations of intradiscal pressures have previously been validated against in vivo disc pressures and loads from instrumented vertebral implants both quantitatively and qualitatively (Bayoglu et al. 2019). Briefly, the same tasks–consisting of the flexion, extension, lateral bending, and axial rotation of the trunk–were simulated, and intradiscal pressures were calculated from predicted disc compressive forces. Good correlations between calculated and measured intradiscal pressures were found for the lumbar (L3/L4 and L4/L5) and thoracic (T6 – T8 and T9 – T11) discs. Furthermore, normalized resultant loads (as the percentage of the load at the upright position) at the L1 vertebra were also in good agreement with the experimental data (Rohlmann et al. 2008). However, the model overpredicted the normalized loads during trunk extension.

Several quasi-static trunk movements were simulated with the model: standing upright, flexion, lateral bending, and axial rotation (see Table 2). Each movement started from the upright posture, lasted ten seconds, and finished when the prescribed movement was achieved. We only specified the motion between the first thoracic vertebra and pelvis segments (hence, no motion at the cervical joints) and distributed the total motion across the thoracolumbar joints in accordance with previous reports (White and Panjabi 1978; Wong et al. 2006; Fujii et al. 2007; Rozumalski et al. 2008; Fujimori et al. 2012, 2014; Tafazzol et al. 2014).

Table 2. Changes in the compressive and shear forces at the T6/T7, T12/L1, and L4/L5 discs due to perturbations in muscle groups.

These values were calculated at the following trunk postures: upright position (0°), 30° and 60° flexion (F 30° and F 60°), 15° and 30° lateral bending (B 15° and B 30°), and 15° and 30° axial rotation (R 15° and R 30°). For each muscle, the maximal increases and decreases in force magnitudes were calculated for its attachment sites along six directions (positive and negative X^g, Y^g, and Z^g-axes). The force magnitudes are given in corresponding disc joint reference frames. Larger changes are emphasized by coloring, where red indicates relatively higher increases and blue indicates higher decreases.

In this study, we only investigated the trunk muscles since our focus was to assess the sensitivity of internal force predictions within this region. In total, 22,032 (918 attachment locations–total number of origins, insertions, and via points–per side × 6 directions × 4 movements) perturbations were simulated. In every simulation, a single attachment point was perturbed from its nominal location along one direction, and all the other attachments were left intact (Figure 1b). However, the attachment points shared by different muscle elements were perturbed simultaneously. Since this model was built right to left symmetrical on the sagittal plane, perturbations of the muscles residing at the right side were also reflected at the left side. Each location was perturbed ± 2.5 mm along the global (${\psi }_{\text{global}}$) +X^g, +Y^g, and +Z^g-axes.

We computed the relative changes (as the percent of the nominal value) in predicted compressive and shear forces at the middle-thoracic (T6/T7), thoracolumbar (T12/L1), and lumbar (L4/L5) discs due to the perturbations. The compression force was calculated along the y^j-axis of the intervertebral joint reference frame (${\psi }_{\text{joint}}$) while the anterior-posterior (AP) and medial-lateral (ML) shear forces were the components along the x^j and z^j-axes, respectively (Figure 1). The percent changes in the disc forces were computed for standing upright (0°), 30° and 60° flexion, 15° and 30° lateral bending, and 15° and 30° axial rotation.

Similar to earlier studies in the lower extremity (Redl et al. 2007; Carbone et al. 2012), a metric–overall sensitivity index–was calculated for each perturbation. This index is essentially a measure of the combined reaction of all muscle elements. Overall sensitivity index (OSI) is calculated as the total percent change in the muscle momentum contributions by using Eq. 1(1) $$\mathit{OSI}=\frac{{\Sigma }_i{\int }_0^T|F_{\mathit{perturbed},i}^{\mathit{mt}}\left(t\right)-F_{\mathit{nominal},i}^{\mathit{mt}}\left(t\right)|\mathit{dt}}{{\Sigma }_i{\int }_0^T{F}_{\mathit{nominal},i}^{\mathit{mt}}\left(t\right)\mathit{dt}}\times 100\%$$(1) : (1) $$\mathit{OSI}=\frac{{\Sigma }_i{\int }_0^T|F_{\mathit{perturbed},i}^{\mathit{mt}}\left(t\right)-F_{\mathit{nominal},i}^{\mathit{mt}}\left(t\right)|\mathit{dt}}{{\Sigma }_i{\int }_0^T{F}_{\mathit{nominal},i}^{\mathit{mt}}\left(t\right)\mathit{dt}}\times 100\%$$(1)

In this equation, $F_{\mathit{perturbed},i}^{\mathit{mt}}(t)$ and $F_{\mathit{nominal},i}^{\mathit{mt}}(t)$ stand for the perturbed and nominal muscle-tendon actuator forces (averages of the left and right sides, i: muscle element) and t is the duration between two consecutive solution steps (t = 1 second and T = 10 seconds). In Eq. 1(1) $$\mathit{OSI}=\frac{{\Sigma }_i{\int }_0^T|F_{\mathit{perturbed},i}^{\mathit{mt}}\left(t\right)-F_{\mathit{nominal},i}^{\mathit{mt}}\left(t\right)|\mathit{dt}}{{\Sigma }_i{\int }_0^T{F}_{\mathit{nominal},i}^{\mathit{mt}}\left(t\right)\mathit{dt}}\times 100\%$$(1) , predicted forces due to perturbations were subtracted from the corresponding nominal values, and the differences were integrated over each time step. To calculate the percent change, the cumulative change in momentum was divided by the nominal cumulative momentum. To avoid zero denominators, nominal muscle forces were set to 0.1 N when smaller than 0.1 N). Finally, grand maximal values of OSI were obtained for each muscle from their corresponding elements along six directions. OSIs were calculated for standing upright, flexion (60°), lateral bending (30°), and axial rotation (30°).

All analyses were run on a desktop computer with 32 cores of an Intel Xeon Gold 6130 (Intel Core, Intel, Santa Clara, CA) and 32 gigabytes of RAM. The simulations were batch-processed and post-processed in Python 3.6 environment using AnyPyTools library (Lund et al. 2019).

3. Results

Changes in the compressive and shear disc forces due to perturbations were calculated for each trunk posture and are presented in Table 2. For each muscle, the maximal increases and decreases in force magnitudes were obtained for its attachment sites along six directions (positive and negative X^g, Y^g, and Z^g-axes). In this table, changes are emphasized by coloring, where red indicates relatively higher increases and blue does higher decreases. Furthermore, nominal disc reaction forces and the corresponding minimum and maximum values are given in Table 3.

Table 3. Nominal intervertebral disc reaction forces [N].

		Upright	Flexion (30°)	Flexion (60°)	Lateral bending (15°)	Lateral bending (30°)	Axial rotation (15°)	Axial rotation (30°)
T6/T7	AP shear	8 (−3, 37)	−49 (−57, −41)	−65 (−84, −58)	17 (9, 36)	25 (14, 41)	42 (33, 57)	66 (40, 111)
	Compression	384 (357, 447)	344 (330, 354)	365 (344, 384)	408 (398, 459)	453 (443, 512)	420 (411, 477)	529 (519, 574)
	ML shear	0 (0,0)	0 (0,0)	0 (0,0)	−12 (−20, −4)	−30 (−42, −26)	−7 (−16, 1)	−15 (−52, 11)
T12/L1	AP shear	24 (21, 33)	−66 (−70, −62)	−248 (−256, −243)	26 (23, 33)	33 (31, 40)	27 (24, 33)	37 (29, 42)
	Compression	425 (409, 444)	584 (579, 590)	738 (706, 806)	439 (430, 457)	460 (456, 482)	459 (443, 486)	534 (507, 569)
	ML shear	0 (0,0)	0 (0,0)	0 (0,0)	18 (11, 20)	29 (20, 32)	6 (2, 9)	4 (2, 13)
L4/L5	AP shear	−294 (−350, −263)	−456 (−463, −450)	−592 (−620, −576)	−335 (−410, −325)	−440 (−483,−426)	−322 (−378, −305)	−352 (−426, −337)
	Compression	503 (477, 546)	821 (815, 828)	1245 (1213, 1311)	545 (533, 613)	652 (631, 693)	535 (522, 587)	586 (567, 671)
	ML shear	0 (0,0)	0 (0,0)	0 (0,0)	−24 (−27, −17)	−43 (−51, −33)	22 (18, 25)	41 (37, 45)

Minimum and maximum values due to perturbations are given inside the parentheses, respectively. The force magnitudes are given in corresponding disc joint reference frames.

The perturbations in most of the muscle groups had only negligible effects (typically below 5%) on predicted disc forces. However, the increases and decreases as high as 353% and 172% respectively were found for shear forces. These were 17% and 7% for the compression forces. Disc forces were most sensitive to perturbations in certain muscle groups located in the rib cage and lumbar regions. The perturbations in the psoas major attachments caused the largest increase in the compressive forces (84 N) at the L4/L5 joint during 30° axial rotation, and longissimus thoracis and rotatores resulted in the largest decrease (32 N) at the T12/L1 joint during 60° flexion. Intercostales externi caused the highest increase in the shear forces (45 N) at the T6/T7 joint, and quadratus lumborum led to the highest decrease (74 N) at the T12/L1 joint, both during 30° axial rotation. When all the joints and postures were considered, on average, quadratus lumborum had the most influence on predicted disc forces, followed by intercostales interni, longissimus thoracis, and intercostales externi. Furthermore, disc forces at the L4/L5 joint were the most sensitive against muscle attachment sites, followed by T6/T7 and T12/L1 joints. Although higher force magnitude changes were found for the compressive forces, percent changes were generally higher for the shear forces (anterior-posterior and medial-lateral). The computed disc forces were also influenced by trunk angle and posture but the effect differed between different muscle groups.

Grand OSI corresponding to each muscle is shown as stacked columns in Figure 2. Similar to the changes in disc forces, relatively higher values were found for certain muscles (quadratus lumborum, intercostales interni, intercostales externi, longissimus thoracis, and iliocostalis lumborum) located in the rib cage and lumbar regions. For all other muscle groups, cumulative OSI was below 10%.

Figure 2. Grand maximal values of OSI are shown as stacked column charts. In this depiction, the height of each bar corresponds to the sum of the maximal values for each muscle in each movement.

4. Discussion

It was previously asserted that insecurities in data and assumptions embedded in musculoskeletal models of the spine, such as vertebral geometry, ligament material properties, and spinal movement, would affect predictions of internal forces. Zander et al. (2016) found that variations in the center of rotation of the intervertebral discs considerably affected muscle activities and disc forces in the lumbar spine. Arshad et al. (2016) investigated the effects of using different spinal rhythms and intra-abdominal pressures. They concluded that there were significant effects on disc and muscle forces and that these increased with trunk flexion. Putzer, Ehrlich, et al. (2016) found that anatomical variations of vertebral body height, disc height, transverse process width, and the curvature of the lumbar spine greatly influenced the joint force predictions at the L4/L5 disc. Recently, Ignasiak, Ferguson, et al. (2016) found that modeling the thorax as a single lumped rigid body instead of a flexible structure resulted in moderately lower compressive disc forces but significantly changed muscle forces. Furthermore, Nussbaum et al. (1995) studied the effects of inter-individual variations in muscle lines-of-action (in a simple biomechanical model) on muscle and joint forces at an erect posture under externally applied moments on the frontal and sagittal planes. They found that predicted L3/L4 compression forces could change by more than 100 N while L3/L4 shear forces could change by more than 50 N.

In the current study, we investigated the sensitivity of muscle and intervertebral disc force predictions to potential errors in modeling muscle attachment sites (muscle origin, insertion, and via points). For this, we used a previously validated model of the entire human spine (Bayoglu et al. 2019), perturbed each attachment location ± 2.5 mm along six anatomical directions, and quantified the changes in muscle and disc forces during upright position, flexion, bending, and rotation of the trunk. Similar to previous studies (Redl et al. 2007; Carbone et al. 2012), we computed OSI for each perturbation. This index is a measure of the collective reaction of all muscles to maintain the joint moments. The error magnitude studied in this work can be introduced during the measurement of the attachment locations in cadaveric studies and their implementation in musculoskeletal models or due to inter-individual anatomical variations. Based on our experience in measuring and defining attachment locations of spinal muscles, this is a large but realistically possible error.

Our simulations indicated that small errors in modeling the muscle attachment sites caused slight to considerable changes in the computed disc forces (Tables 2 and 3). Although for the majority of muscles the changes were minor (less than 5%), for some muscles, for example, quadratus lumborum, the change in disc shear and compressive force was as high as 353% and 17%, respectively. In general, we found that the percent changes in disc shear forces were higher for the changes in the muscle lines-of-action than the compression forces. This was attributed to the fact that disc compressive forces were much larger than disc shear forces. Thus, for the same absolute force change, we would observe larger percentile changes in shear forces. Furthermore, intervertebral disc forces at the L4/L5 joint were the most sensitive against muscle attachment locations, and the T6/T7 joint showed higher sensitivity than the T12/L1 joint.

Our analyses also suggest that muscle recruitment patterns were most sensitive for the attachments of certain muscle groups positioned in the rib cage and lumbar spine (Figure 2). These muscles were quadratus lumborum, intercostales interni, longissimus thoracis, and intercostales externi. These muscle groups have multiple elements that attach to several bones and contribute to the moments at the costo-sternal, compound costo-vertebral and costo-transverse, and lumbar intervertebral joints during trunk movements. These findings indicate that the modified moment arms made these muscles less suitable to contribute to these joint moments and favored other muscles for maintaining the needed joint moments (Redl et al. 2007; Carbone et al. 2012). Similar to the changes in the disc forces, relatively higher sensitivities were found for these muscles. This indicates the dependence of spinal loads to muscle recruitment patterns and hence to muscle attachment locations. Other muscle groups had much lower sensitivity values implying that errors involved in modeling their locations would have a much smaller impact on muscle function and internal loads. Moreover, Figure 2 depicts that the changes in the muscle forces were also dependent on the simulated tasks. The computation of muscle forces was most sensitive in trunk lateral bending followed by upright posture, trunk rotation, and trunk flexion. In a similar study, Nussbaum et al. (1995) elaborated that spinal forces were substantially altered for moderate changes in the lines-of-action of the oblique muscles but were minimally changed for the vertically oriented muscles. We also found that disc forces were sensitive to external and internal oblique muscles as well as some vertically oriented muscles such as iliocostalis lumborum, longissimus thoracis, and quadratus lumborum. However, our analyses did not identify the oblique muscles as being most influential nor vertically oriented muscles as least influential. The discrepancies in the findings can perhaps be due to using different spinal morphology and modeling practices or to simplifications of the muscular anatomy.

Several (limiting) assumptions affect the findings of this study. Firstly, we used a simple muscle model which only took into account isometric strengths for modulating force exertion. We choose to not to represent force-length and force-velocity relationships as this would require further calibrating muscles’ optimal fiber and tendon slack lengths for which there is no data. Secondly, the movement at the intervertebral disc articulations was simplified by using spherical joints which were constrained for translational DoFs. Ghezelbash et al. (2015) found low to moderate impact on model predictions and suggested that intervertebral translational DoFs can be ignored for small trunk angles in flexion. Furthermore, Senteler et al. (2018) investigated the sensitivity of disc force predictions to center of rotation (COR) locations. They asserted that a posterior COR in an upright position and an anterior COR in flexion would optimize the lumbar joint loads. Thirdly, we did not incorporate spinal ligaments but implemented their collective rotational stiffnesses linearly in the intervertebral and rib cage joints. Thus, excluding the ligaments or simplifying their non-linear mechanical behaviors could alter the load sharing mechanism between the muscles and affect our findings (Putzer, Auer, et al. 2016; Ghezelbash et al. 2018). For example, Putzer, Auer, et al. (2016) found that implementing stiffer ligaments could cause a shift of the loads to the lower lumbar disc levels. Fourthly, the facet joints between the adjacent vertebrae were not modeled. It was previously discussed that these joints will decrease the forces carried through the anterior column during trunk extension (Dreischarf et al. 2013; Bruno et al. 2015). Extension was, however, not investigated in this study. Fifthly, the data presented in this work was obtained by using a subject-specific model of the spine where anatomical variations in spinal geometry and muscle morphology were not accounted for. Previous studies found that the spinal curvature had a major influence on calculated internal loads (Briggs et al. 2007; Putzer, Ehrlich, et al. 2016; Bruno et al. 2017). Briggs et al. (2007) reported strong correlations between thoracic spine curvature and the vertebral loads and thoracic curvature and muscle forces. They discussed that a higher kyphosis angle resulted in larger normalized flexion moments and the disc forces. Thus, generalizing our findings to other musculoskeletal representations of the spinal morphologies requires further study. Lastly, we used a frequently used optimization criterion which minimized muscle fatigue for muscle recruitment. Implementing a different performance criterion might affect the findings. Nevertheless, we do not feel these assumptions will alter the qualitative behavior of the model much.

The majority of the previous generic musculoskeletal models of the spine relied on the anatomical illustrations provided in the previous dissection studies for defining muscle lines-of-action. Our findings suggest that for most muscles this only moderately affects the results, but for some muscles this approach might lead to erroneous evaluation of the muscle function and spinal loads. We would thus recommend using anatomical datasets in which muscle attachment locations were measured precisely or using imaging methods to acquire muscle paths. In this line, certain care should also be taken when morphing an anatomical dataset for developing a personalized model of the spine, particularly for the muscles which affect spine biomechanics to a large extent as identified in this study.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This study was supported by a research grant by fonds NutsOhra (1302-032).