Muscle-driven and torque-driven centrodes during modeled flexion of individual lumbar spines are disparate

Lumbar spine biomechanics during the forward-bending of the upper body (flexion) are well investigated by both in vivo and in vitro experiments. In both cases, the experimentally observed relative motion of vertebral bodies can be used to calculate the instantaneous center of rotation (ICR). The timely evolution of the ICR, the centrode, is widely utilized for validating computer models and is thought to serve as a criterion for distinguishing healthy and degenerative motion patterns. While in vivo motion can be induced by physiological active structures (muscles), in vitro spinal segments have to be driven by external torque-applying equipment such as spine testers. It is implicitly assumed that muscle-driven and torque-driven centrodes are similar. Here, however, we show that centrodes qualitatively depend on the impetus. Distinction is achieved by introducing confidence regions (ellipses) that comprise centrodes of seven individual multi-body simulation models, performing flexion with and without preload. Muscle-driven centrodes were generally directed superior–anterior and tail-shaped, while torque-driven centrodes were located in a comparably narrow region close to the center of mass of the caudal vertebrae. We thus argue that centrodes resulting from different experimental conditions ought to be compared with caution. Finally, the applicability of our method regarding the analysis of clinical syndromes and the assessment of surgical methods is discussed.


Introduction
Relative movement between (lumbar) vertebrae occurs during most daily motions. In three dimensions, such a relative movement between two time instances can be represented by a finite helical (or screw) axis (FHA), i.e., an (instantaneous) axis of rotation that points toward the direction of possible translations. By intersecting the FHA with an anatomical plane, an instantaneous center of rotation (ICR) is obtained. The time evolution (path) of the ICR, the socalled centrode, had been investigated by multiple researchers with highly diverse and ambitious objectives, e.g. (1) recognizing general patterns in order to "give base-line references for potential diagnostic applications" , "predicting...injurious vectors" (Qiu et al. 2003), and finding an "indicator for mechanical disorders" (Schmidt et al. 2008) or "motion characteristics of the normal lumbar spine" (Yoshioka et al. 1990), (2) evaluating "the quality, rather than the quantity, of cervical spine movement" (Baillargeon and Anderst 2013), (3) hoping for the ICR to be "interpreted in terms of...anatomical and pathological Electronic supplementary material The online version of this article (https ://doi.org/10.1007/s1023 7-020-01382 -9) contains supplementary material, which is available to authorized users.

3
factors" (Bogduk et al. 1995), (4) describing the change in ICR location as a consequence of disk degeneration (Cossette et al. 1971;Ellingson and Nuckley 2015;Gertzbein et al. 1985), (5) attempting to relate the ICR location to the "choice of anterior and posterior instrumentation" (Haher et al. 1991) or certain implant parameters (Niosi et al. 2006), (6) demonstrating that "analysis for sagittal plane motion of the lumbar spine is possible" (Ogston et al. 1986), and (7) correlating ICR paths to facet forces (Rousseau et al. 2006). However, a recent review (Widmer et al. 2019) had revealed that, up to now, the ICR provides only faint criteria for the description of spinal kinematics under healthy and degenerative conditions.
In this work, we used elementary, individualized multibody simulation (MBS) models of the lumbar spine (Damm et al. 2019), performing flexion movements, to introduce a statistical criterion that may serve as a first step toward describing, detecting, and eventually understanding the cause and effect of the centrode's location: the (weighted) confidence ellipse as introduced in Sect. 2.4. The underlying assumption was that similar (individual) spinal structures, under the same loading conditions, yield a similar relative motion of segments. Of course, changes in the individual geometries cause different relative motion patterns and thus different centrodes. Yet, we show that these differences were small under the same loading condition, but distinguishable under varying conditions. Our method was utilized to identify centrode locations depending on (1) individual geometries, (2) force application modes, and (3) material properties that are influenced by clinical syndromes. First, the location of the centrode across individual spine geometries was calculated, both with and without preload representing the upper body weight. This individualized modeling is particularly worthwhile, since, on the one side, modelers from MBS and finite element (FE) communities usually employ generic models (Abouhossein et al. 2013;Qiu et al. 2003;Senteler et al. 2017;Schmidt et al. 2008), which cannot account for structure-based deviations. Experimenters, on the other side, conduct helpful individual measurements, but have no elaborated model on hand (Cossette et al. 1971;Gertzbein et al. 1985;Haher et al. 1991;Niosi et al. 2006;Ogston et al. 1986). Second, the differences in conducting physiologically-based (muscle-driven) and artificial (torquedriven) movement on the centrode were investigated, likewise with and without preload. Third, our method could aid to assess the influence of clinical syndromes or treatments on the centrode's location, as is discussed on the example of modeling the surgical fixation of vertebrae.

Model
In total, seven individual lumbar spine models (L1-SA, cf. Fig. 1) were investigated. The term 'individual' here  and muscle (red lines) insertion points were set at typical anatomic landmarks. Ligaments, intervertebral disks (gray ellipsoids) and facet joints (blue planes) transmit forces. Passive muscles are depicted in a paler red tone, whereas active muscles were depicted in vibrant red. The line of action of the preload force, representing the upper body weight, was located between the femoral heads (gravity line, silver). a, b The neutral position in lateral view, together with the acting forces ( F preload , F muscle ) and an external torque ( T external ), respectively. c Depicts the corresponding ventral view. Slight asymmetries, resulting from individual geometries can be detected, see also 1 3 refers to the vertebral surfaces, including inter-vertebral distances, and ligament insertion points that were extracted from digital image data (DICOM). Geometries of healthy patients (six male and one female with an age of 32.6 ± 7.04 years) were provided by the university clinic in Mainz. After semi-automatic segmentation, vertebral surfaces were loaded into a MBS tool (SIMPACK: Dassault Systèmes, Vélizy-Villacoublay, France) under preservation of in vivo distances for intervertebral disks and facet joints. Ligament insertion points were set manually (Schünke et al. 2015) and counter-checked by neurosurgeons from the university clinic in Mainz. All passive force-transmitting structures, i.e., intervertebral disks (IVD), ligaments and facet joints, were modeled as nonlinear spring-damper elements (Damm et al. 2019). Here, particularly force-length and force-angle characteristics of ligaments and IVD, respectively, were extracted from step-wise reduction experiments (Heuer et al. 2007) and validated using data on range of motion (Heuer et al. 2007) as well as intradiscal pressure (Wilke et al. 1996). Ligaments were found to be significantly less stiff than suggested by classical data sets, e.g. (Chazal et al. 1985;Shirazi-Adl et al. 1986;White and Panjabi 1990). Indeed, recent work, combining in vivo experiments and computer simulation, support this finding (Mörl et al. 2020). Ligament pre-strain was set according to literature data, cf. (Damm et al. 2019, sect. 2.2.4).
Active forces were transmitted by Hill-type muscle models (Guenther et al. 2007;Haeufle et al. 2014;Guenther 2016, 2018) of M. psoas major and M. multifidus, each strand modeled as a one-dimensional point-topoint element with no deflection. Since neither kinematic nor kinetic data were available to conduct an elaborated parameter estimation, tendon and fiber lengths parameters were taken from literature, particularly (Christophy et al. 2012). For every muscle strand, values for optimal fiber length and tendon slack length were adapted from their Table 1, columns seven and nine, respectively. In order to maintain the fiber-to-tendon length ratio as an important functional measure (Mörl et al. 2015), both quantities were equally scaled to match individual geometries, e.g. each increased by 5% if the distance from origin to insertion was 5% higher than the literature reference. The fiber-to-tendon length ratios for each muscle strand were 2.85 for M. psoas and 2.47 for M. multifidus. For one model, exemplary optimal fiber lengths and tendon slack lengths for both muscle groups are displayed in Table 1. M. psoas maximum force was set to 80 N for each strand, the median value from Christophyet al. (Christophy et al. 2012, table 1, fourth column) for non-IVD parts. M. multifidus force was set to 21 N for each strand, likewise the median value for non-laminar parts. Remaining parameters of the muscle model were assumed to be generic and taken from (Guenther et al. 2007, table 2) Table 1 Optimal fiber lengths (OFL) and tendon slack lengths (TSL) of each strand within the M. psoas and M. multifidus group, exemplary for one lumbar spine model Strands from M. psoas (Ps) are labeled with the corresponding insertion point at the lumbar spine, whereas strands of the M. multifidus (Mf) are labeled with both origin and insertion, cf. Fig. 1. The listing serves the purpose of estimating the magnitudes of each muscle strand as well as demonstrating the slight asymmetry of sinister and dexter geometries. The fiber-to-tendon length ratios (OFL/TSL) for each group, which were applied for all seven models, are also given

Calculating the instantaneous FHA and ICR during flexion
For comparability, we hereinafter investigated only flexion of the modeled lumbar spines, i.e., the movement that occurs during forward bending of the upper body, see Fig. 1d and the supplementary video file. This flexion was driven either by muscle forces or an external torque applied at the COM of vertebra L1 (Abouhossein et al. 2013), see also next Sect. 2.3. We determined the finite helical axis (FHA) between two vertebrae using a least-squares method on spatial coordinates (Kwon 2008;Spoor and Veldpaus 1980). Therefore, four markers (ligament insertion points of ligamentum supraspinale, flavum sinister, intertransversale sinister and anterior longitudinal dexter) on each vertebral body were tracked at each time instance (every millisecond), relative to the subjacent vertebra. In a subsequent step, the intersection point of the FHA with the corresponding anatomical plane-in our case the sagittal plane-was calculated and defined to be the instantaneous center of rotation (ICR). Figure 2 depicts the situation for an exemplary L4-L5 segment. The reduction to a two-dimensional quantity is reasonable here, because the observed motion (flexion) only takes place in the sagittal plane itself. This means the vertebrae can be assumed to undergo no substantial relative translation along the helical axis, which ought to be approximately perpendicular to the sagittal plane, and thus perform a pure rotation. We further note: First, the term 'instantaneous' is rather to be read as 'approximately instantaneous' or 'finite', because no rigorous differential-geometric approach has been applied. However, observed time intervals were small and thus this term was chosen in order to retain consistency with the literature. Second, most in vivo/vitro studies (Bogduk et al. 1995;Cossette et al. 1971;Haher et al. 1991;Yoshioka et al. 1990) and even elaborated FE models (Qiu et al. 2003;Shirazi-Adl et al. 1986;Schmidt et al. 2008) calculate only elementary, two-dimensional rotation with Reuleaux' method (Reuleaux 1875), which is neither applicable in arbitrary three-dimensional motion nor captures translational movement of vertebrae.

Simulation protocol: muscle-driven and torque-driven flexion
In order to investigate the differences in the paths of the ICR, when either driven by muscle (psoas) forces or 10 Nm torque, we defined four different simulation scenarios of lumbar spine flexion. Two scenarios with and two without 500 N preload, representing the upper body weight: (1) no preload and muscle-driven (NPM), (2) preload and muscle-driven (PM), (3) no preload and torque-driven (NPT), (4) preload and torquedriven (PT). Pre-load forces were supposed to act on L1 and along the gravity line (Le Huec et al. 2011), i.e., the line that vertically bisects the two femoral heads and is thought to go through the center of gravity of the upper body, cf.  Since CT data were obtained from decumbent patients, those scenarios were supposed to represent lying (unloaded) and standing (preloaded) positions. Torque was applied on the COM of the upmost vertebral body (L1). Total simulation time was 2.5 s in each scenario for each spine. M. psoas major started in a fully activated equilibrium ( q ≈ 1 ) and was fully stimulated ( = 1 ) for the whole simulation period, for the notation see (Rockenfeller and Guenther 2018). M. multifidus was neither activated (q = q 0 ≈ 0.01) nor stimulated ( = 0) and therefore had only passive, antagonistic contributions on the ICR paths in all cases, cf. also (Zwambag,Brown 2020). Although the high stimulation value for M. psoas major, is far from physiological reality, it was chosen (a) for comparability of all simulation outputs and (b) to obtain a significant motion at all with this reduced model, see also the discussion.

Confidence ellipses for the instantaneous center of rotation
To quantify the effects of individual geometry on ICR paths, we calculated a 95% confidence ellipse for each pair of adjacent vertebra in each of the afore-defined scenarios. The center point of the ellipses was calculated as the mean of the twodimensional path data, and the semi-axes were represented by the two eigenvectors of the respective covariance matrix (Draper 1998;Galton 1886;Spruyt 2014). The lengths of the semi-axes thus indicate the corresponding standard deviations. We noted that variations in movement of the vertebrae were larger at the beginning than at the end of the 2.5 s simulation, because the flexed equilibrium had been reached. To account for that fact, we additionally calculated exponentially weighted confidence ellipses of the n (in our case n = 2500 ) data points by the discrete mapping w ,n ∶ {1, … , n} → [0, 1] , with Weighted mean ( ̄ w ) and covariance matrix ( Here, = 10∕n was found to be a suitable compromise between the location and the narrowness of the final confidence ellipses. The larger the more weight is attributed to the early data points. In the non-weighted case, w ,n (k) = 1∕n for all k ∈ {1, … , n}.   Fig. 3a shows the ICR-time courses-and corresponding confidence ellipses for each level of the lumbar spines-that resulted from simulating the muscle-driven scenario without any preload (NPM). In each level of the individual, muscle-driven lumbar spine models, the ICR-time courses showed similar behavior: starting inferior-posterior to the center of mass of the caudal vertebra and moving in a superior-anterior direction. The first semi-axis lengths of the corresponding (non-weighted) ellipses increased with caudal position and were significantly longer (up to factor 5) than the second semi-axis lengths, cf. Table 2. Noticeably, the orientation angle of these ellipses did not vary substantially across spine levels ( 32 • -37 • ). The lower the spine level (L1-L2 to L5-SA), the more ellipse centers were found to move in superior-anterior direction, with the center of the L3-L4 ellipse almost congruent to the L4 COM. Introducing the weight function (Eq. (1)) for the purpose of including early ICR data, resulted in an increase in orientation angle ( 33 • -50 • ) as well as an increase particularly in first semi-axis lengths. ROM values from 5.1 • to 14.5 • were found comparable to the NPT scenario.

Results
NPT Fig. 3b shows the ICR-time courses that resulted from simulating the torque-driven scenario without any preload (NPT). In any level of the torque-driven spine models, with exception of L5-SA, the ICR-time courses exhibited a similar behavior: for the whole time horizon of simulation, the ICR was located in a narrow region, superior-anterior and close to the COM of the caudal vertebra. Consequently, the semi-axis lengths were distinctly shorter than in the muscle-driven scenarios and the quotients between the first and second semi-axis lengths were substantially smaller, cf. Table 2. Eventually, this quotients became approximately 1 on the L4-L5 level, transforming the ellipses almost into circles. Contrary to the muscle-driven scenarios, in the NPT scenario the weighted ellipses were even  Table 2 ICR confidence ellipse parameters for muscle-and torque-driven flexion scenarios with and without preload (column one: NPM, NPT, PM, PT) as displayed in Fig. 3 Column two: range of motion (ROM) for the whole lumbar spine models. Columns 4-9: center points, first and second semi-axes lengths, as well as orientation angles for the non-weighted (columns 4-6) and weighted (columns 7-9) ellipses at each spinal level (L1-L2 to L5-SA) narrower than the non-weighted ellipses, indicating that the variation in the centrode's location took place well after the beginning of simulation. Due to the narrowness and the resulting indistinguishability of the two semi-axes, the orientation angles became a random number on the interval −90 • … + 90 • . ROM values of around 10.7 • were already mentioned to be comparable to the NPM scenario, but also matched the literature data of cadaver experiments (Heuer et al. 2007, fig. 4). PM Fig. 3c shows the ICR-time courses that resulted from simulating the muscle-driven scenario with a preload of 500 N (PM) acting on L1 and along the gravity line, as described in Sect. 2.3. The different levels of the individual loaded lumbar spine models showed a comparable ICR-time courses as in the NPM mode, i.e., the ICR-time courses started inferior-posterior to the COM of the caudal vertebra and moved in a superior-anterior direction. Likewise, the lower the spine level, the more ellipse centers moved in superior-anterior direction, here with the center of the L2-L3 ellipse almost congruent to the L3 COM. However, in levels L1-L2 and L2-L3, few early ICR data points lay markedly posterior to the COM, thus forming a tail-shaped path toward the ellipse center. Consequently, the respective first semi-axis lengths of the weighted ellipses were longer than for the non-weighted ellipses, and for both cases longer in the PM than in the NPM scenario. Yet, the width of the ellipses (the second semi-axis lengths) showed no such differences and likewise increased in caudal direction, cf. Table 2. Ellipse centers for the PM scenario were found to be near to those from the NPM scenario. The orientation angle of the non-weighted and weighted ellipses increased in caudal direction from 7 • and 9 • to 33 • and 40 • , respectively. The latter angles (L5-SA) were thus similar to the NPM scenario. The ROM ranged from 25.5 • to 39.1 • and was comparable to the ROM of the PT scenario. Furthermore, the ROM was about four times larger in the scenarios with than without preload.
PT Fig. 3d shows the ICR-time courses that resulted from simulating the torque-driven scenario with a preload of 500 N (PT). Comparable to the NPT scenario, each level of the PT models yielded similar centrodes, which were located in a narrow region, superior-anterior and close to the COM of the caudal vertebra. Here, the lengths of the first and second semi-axes were on average 1.5-2 times longer than in the NPT scenario, but still small compared to the muscle-driven scenarios. As in the NPT scenario, the range of quotients between the first and second semi-axis lengths was smaller in NPM ( 1.4 … 2.7 ) than in the muscle-driven scenarios ( 1.5 … 14 ). In addition, there was no noticeable difference between the semi-axis lengths of the weighted and non-weighted ellipses. Contrary to the muscle-driven scenarios, the orientation angles of the non-weighted ellipses were negative, except for the level L4-L5. The orientation angles of the weighted ellipses increased, as in the PM scenario, in a caudal direction from −78 • to 39 • . ROM values between 32.6 • and 41.6 • were comparable to the PM scenario. In both NPM and NPT, the size of the ellipses increases in caudal direction, indicating an amplified translation instead of pure rotation (the narrower the ellipse, the purer the rotation).

On torque-driven experiments and physiological insights
Bending forward, reaching for a crate of beer and lifting it up incorrectly may result in high loading peaks within the lumbar spine (Nachemson 1965;Wilke et al. 1999). Due to unfavorable lever arms with regard to the joints, the gravitational forces of the body parts can cause high torques within the human body, which must be compensated for by the muscles. As a consequence, the internal muscle forces, transmitted according to Newton's law (actio=reactio), eventually generate high compressive and shear forces in the various spinal structures. When experimentally isolating certain structures (e.g. vertebrae, IVDs and ligaments) while leaving others out (e.g. muscles), experimental observations may greatly differ from physiological reality. For example, isolated (lumbar) spines had been observed to "buckle" under compressive loads not even close to in vivo magnitudes (Crisco 1989). To compensate for the missing supporting structures, when applying high loads on cadaver specimen, the concept of follower load was established (Patwardhan et al. 1999;Rohlmann et al. 2009a). At this, a guiding rail ought to ensure a purely compressive force transmission and prevent the occurrence of shear forces, which would lead to buckling. Consequently, the forces were applied "tangent to the spinal curve, passing through the center of rotation of each segment" (Patwardhan et al. 1999, fig. 1), which was supposed to be located perfectly in between the vertebrae. It was found that the path of the follower load influences the model output and should be optimized in a sense that it passes through the centers of rotation between vertebrae (Dreischarf et al. 2010).
In both in vitro experiment and follower-load model, torque had not been physiologically induced by muscle forces, but artificially "applied" by spine testers to induce a flexion movement. We herein revealed that this method differs substantially from muscle-driven movements by means of the corresponding centrodes: torque-driven centrodes, regardless the individual spine geometry, can be found in a narrow region superior-anterior to the caudal vertebra's COM. Muscle-driven centrodes show a more individual behavior and stretch over a wider range. These observations hold true for non-preload and preload scenarios alike, under our model assumptions. Of course, our presented model is far from capturing every physiological aspect of in vivo force transmission, since we omitted most of the whole body's structures. However, our approach may serve as a starting point for pursuing centrode-based investigations. For example, the herein introduced concept of confidence ellipses may be utilized to assess the influence of model parameter changes (sensitivities) on the centrode location. These investigations could include the influence of ligament stiffness or failure (Abouhossein et al. 2013;Alapan et al. 2013;Putzer et al. 2016), joint forces (Senteler et al. 2017), implant positioning (Dreischarf et al. 2015;Rohlmann et al. 2010) or variable load application (Rohlmann et al. 2009b).
On the one hand, our findings support the application of the follower load concept for recreating in vitro experiments: torque-induced centrodes (and thus the path of the follower load) are virtually inert to individual geometries, instants of time, or loading scenarios. Hence, once established, the follower load can remain unchanged during the whole simulation. On the other hand, our findings speak against the application of the follower load concept for recreating in vivo experiments: the ICR location is known to change during physiological motion (Aiyangar et al. 2017, fig. 2), as is likewise visible in our Fig. 3a, c. This change in ICR location indicates the existence of translational rather than purely rotational movements of vertebrae relative to each other, which cannot be captured while utilizing a follower load.
Summarizing, our findings suggest that modeled spinal motion have to be compared with caution regarding their impetus. When aiming for physiological insights, muscle-driven models ought to be utilized. Here, it might be worth investigating whether the method of muscle control, e.g. inverse dynamic (Happee 1994), EMG-driven (Lloyd and Besier 2003) or forward dynamic (Rupp et al. 2015), has significant influence on the respective centrodes. Likewise, incorporating muscle deflection on larger muscles (Hammer et al. 2019) might lead to altered centrodes. When aiming at recreating in vitro experiments that require stabilizing follower load, torque-driven models constitute a more natural choice, but their results cannot directly be transferred to physiological reality, as statements on possible medical implications (see the Introduction and the next Sect. 4.2) might not satisfy the expectation.

Centrodes from a medical point of view
The study of spinal motion is of utmost importance when aiming to understand the formation of disorders and the effects of surgical interventions. In clinical practice, hypermobile segments or degenerative structures commonly undergo fixing procedures. Stabilizing only one or few segments is hereby generally considered successful, although the individual spinal motion pattern is not examined in detail. If long constructs are required, or when aiming to preserve or restore "healthy" motion by application of dynamic implants (motion preserving implants), an exact balance of the resulting forces and thus profound knowledge of the motion pattern is, however, mandatory. Numerous research on motion patterns has been conducted on the cervical spine (Amevo et al. 1991;Anderst et al. 2015;Wachowski et al. 2017). This area is not only less complex than the lumbar area (due to less soft tissue involved), but also the region where most dynamic implants-especially disk prostheses-are used. The highest loads and consequently the location where degenerative changes occur first, is however be found the lumbar spine (Auerbach et al. 2008).
Changes in (lumbar) spinal kinematics have been observed following surgical procedures (spondylodesis using different techniques (Nomoto et al. 2019), facetectomy (Zeng et al. 2017), implantation of disk prothesis (Yue et al. 2019) or pedicle screw-based dynamic implants (Prud'homme et al. 2015)), but do also occur naturally due to degeneration or trauma (Amevo et al. 1992) as well as in obese patients (Rodriguez-Martinez et al. 2016). In addition, several studies, in vitro and in vivo, have been conducted to analyze lumbar spinal kinematics and to determine the centrode under healthy and degenerative conditions, see (Widmer et al. 2019) for a review. So far, however, there exists neither mechanistic nor statistical criteria linking the mere observation to a quantitative kinematic assessment, let alone to predicting the effects of surgical interventions. Given sufficient experimental data, the herein presented concept of confidence ellipses could help correlating centrodes to their corresponding clinical syndromes or medical treatments.
CT scans and bending fluoroscopy are generally available for most spinal patients. These image data allow for example to assess the grade of disk degeneration (Quint and Wilke 2008), but other structural properties-as the individual stiffness of certain ligaments or the strength of supporting muscles-can only be approximated. Further evaluation of medical image data, e.g. water-fat MRI (Schlaeger et al. 2018), could provide estimates for individualized muscle parameters, such as the maximum force. Further, stereo X-ray films ) could allow for a preciser tracing of the ICR location. Overall, the utilization of individual data constitutes an important step toward an individual spine model that could eventually be used to predict the effect of surgical interventions and to optimize operative plans before surgery (Kantelhardt et al. 2015). Implants and their positions could thus be selected, not solely but among others, on basis of individual centrode simulations, see also the next section.

Consecutive fixation systematically alters the centrode: an exemplary scenario
Modeling load changes in lumbar spines, as a result of surgical implant placement, are utilized to assess the required medical procedure a priori (Xu et al. 2019). Additional criteria could be derived by classifying the respective centrode paths. In the literature, however, we identified two major issues in existing studies. First, lumbar spine models (Kiapour et al. 2012), which are validated using in vivo data (Pearcy and Bogduk 1988), where the center of rotation had been calculated only between the start and the endpoint of the movement, omit all mid-range dynamical information (Dombrowski et al. 2018). Second, models (Abouhossein et al. 2013), which are validated using torque-driven in vitro data (Rousseau et al. 2006), where one of the vertebrae had been fixed, obtain centrodes in between the vertebrae, which is contrary to in vivo findings , fig. 2). These two issues can also be found combined (Naserkhaki et al. 2018;Schmidt et al. 2008).
In this section, we cursory glance at possible consequences of consecutively fixating spinal segments, i.e., inserting rigid implants starting from the sacrum and moving cranially up to L2. Figure 4 visualize the scenario and the results from a single NPT model that was utilized for a purely exemplary purpose, thus no error ellipses were calculated. Although we stated above that muscle-driven models ought to be consulted when aiming at confidable physiological findings, more sophisticated structures are required in our model to allow for significant ROM in the ultimately fixated scenario. Let SA-X denote a fixation from the SA upward to segment X, i.e., SA-L2 refers to a spinal unit where SA to L2 are rigidly connected by implants. The ROM of the flexion induced by 10 Nm expectedly decreased during consecutive fixation: from ROM SA-SA = 10.5 • to ROM SA-L2 = 0.3 • . The corresponding centrodes showed a systematic behavior: the further away segments were from the fixation, the narrower and more inferior-posterior the centrodes were located, ultimately approaching the COM of the caudal body. Within these centrodes, the migration path of the ICR was found to evolve superior-anterior. The segments, which were located directly cranial to the fixation showed a hookshaped centrode right in the middle of the respective IVD.
These findings, although conducted with a rather simplified model, suggest a critical view on the aforementioned issues: mid-range dynamics as well as multi-(> 2 )-body analyses could play a crucial role in assessing operational methods such as optimizing implant positioning (Haher et al. 1991;Niosi et al. 2006). Fig. 4 a Two-dimensional centrode locations for five consecutive degrees of fixation (colored dots), obtained from running NPT scenarios on an individual lumbar spine model. The lateral view on the lumbar spine is shown in the background and crucial regions are highlighted by a zoom. b Semitransparent three-dimensional depiction of the corresponding implant (or pedicle screw) placements within the vertebral bodies, without further active or passive structures visible. Same colors indicate the same degree of fixation: (1) SA-SA (green, standard NPT as in Fig. 3b), (2) SA-L5 (violet), (3) SA-L4 (yellow), (4) SA-L3 (red), (5) SA-L2 (blue). ICR coordinates were computed as described in Sect. 2.2

Supplementary material
Three video files showing an exemplary flexion movement, from neutral position to full bending (cf. Fig. 1a, d), are provided. Each video depicts a different view on the lumbar spine (frontal, lateral and lateral-frontal). At every time instance, the corresponding FHA are shown.
Author's Contribution RR developed the concept of the manuscript, wrote the core text and discussed physiological implications of the method. AM conducted the model calculations and discussed the effects of consecutive fixation. ND set up the model as well as the basic idea of this work. MK and SRK contributed equally on collecting image data, identifying ligament insertion points, and discussing the findings from a medical point of view. RF helped developing the calculation of the FHA and contributed geometrical interpretations of the ICR path locations. KG helped developing the model, provided the MBS graphics and conducted a final review of the manuscript.
Funding Open Access funding provided by Projekt DEAL. A.M. is supported by funding from Mechanical Systems Engineering Laboratory, EMPA-Swiss Federal Laboratories for Materials Science and Technology, Duebendorf, Switzerland.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflicts of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creat iveco mmons .org/licen ses/by/4.0/.