Accurate estimation of joint loading is of high significance to study the biomechanics of the spine. Several musculoskeletal (MSK) spine models have been introduced in literature (de Zee et al., 2007; Christophy et al., 2012; Bruno et al., 2015; Rupp et al., 2015; Cazzola et al., 2017; Mörl et al., 2020; Guo et al., 2021; Silvestros et al., 2022; Meszaros-Beller et al., 2023) with the goal to predict muscle forces, muscle activation patterns and internal loading conditions during human movement.

There are two multibody (MB) approaches typically used to determine these quantities: the forward-dynamics (FD) and the inverse-dynamics (ID) approach. Both approaches aim to provide the optimal solution to the redundancy problem, a mathematical overdeterminacy of the MSK system arising from a greater number of muscles crossing a joint than number of degree of freedom (DOF) specifying the joint movement (Pandy, 2001). Which approach is appropriate to be used depends on the intended purpose of the modelling study.

The main goal of the FD approach is to gain an intrinsic understanding of movement control, i.e., how the nervous system and the muscles use sensory information to produce a coordinated movement. FD works in a forward sense in which muscle forces governed by the underlying muscle-tendon dynamics act on the skeletal geometry in response to muscle stimulation initiated by the central nervous system. The predicted movement is the result of a dynamic interplay of all modelled structural components. The drawback of the FD approach is that it requires accurate modelling of all load-bearing structures, i.e., muscles, ligaments, the intervertebral disc (IVD), and the control scheme predicting individual muscle stimulation. Applied on the spine, this is particularly difficult regarding the high number of articulations and muscles and ligaments involved in stabilising the mechanically instable multi-joint structure of the upright human spine (Smit, 2020).

Due to the structural complexity of the spine, the FD approach has been rarely used to simulate spine biomechanics but has risen to new heights in recent years. A detailed lumbar spine model with a lumped upper body was introduced by Rupp et al. (2015), validated with respect to its passive stiffness properties in a lying posture (Mörl et al., 2020) and recently extended to include articulated thoracic vertebrae in order to quantify the load sharing between individual biological structures of the functional spinal unit under gravity (Meszaros-Beller et al., 2023). Another FD lumbar spine model was presented in Damm et al. (2020). In this study, ligament material properties were reverse engineered from experimental stiffness measurements conducted on cadaveric functional spinal units (Heuer et al., 2007). Both, Damm et al. (2020) and Mörl et al. (2020) highlighted the importance of having physiologically accurate material properties of passive spinal structures, i.e., of ligaments in particular. Further, Müller et al. (2021) used the FD approach to undertake a quantitative investigation of the effect of varying lumbar lordosis angle and muscle activation on the load distribution. Finally, Guo et al. (2021) recently presented a full spine model with an articulated cervical, thoracic and lumbar region to demonstrate the muscle activity-dependent change in intra-abdominal pressure and and its unloading effect on the spine.

On the other hand, the ID approach serves a more descriptive, analytic purpose: the computation of generalised net joint forces and torques for a specific motion task of interest (Pandy, 2001) and the corresponding net contribution of muscle forces to the joint loads which are assumed to compensate the remaining net joint forces and torques from the ID analysis. Thus, ID results are often interpreted as net muscle forces acting in the functional spinal unit. Specifically, ID works in an inverse sense in which the motion is a priori known from experimentally acquired motion capture data to derive segment body velocities and accelerations. In more advanced scenarios, normalised electromyography recordings monitoring the muscle activity for the desired movement are taken into account to assist in solving the optimisation problem in ID (Pizzolato et al., 2015; Molinaro et al., 2020; Silvestros et al., 2022).

The majority of MSK spine models today make use of the ID approach. The detailed characterisation of back muscles in Christophy et al. (2012) set the baseline in modelling of the spine upon which many current models are based (Senteler et al., 2014; Bruno et al., 2015; Dao et al., 2015; Rupp et al., 2015; Raabe and Chaudhari, 2016; Mörl et al., 2020; Overbergh et al., 2020). Bruno et al. (2015) combined, revised and completed model geometry and muscle architecture based on previous models (Vasavada et al., 1998; Christophy et al., 2012) to develop a fully articulated thoracolumbar spine model. Despite the impressive anatomical detail of the spine models presented by Christophy et al. (2012) and Bruno et al. (2015), the studies lacked physiological kinematic input. Generally, the advantage of the ID approach is that it allows analysis of advanced motion tasks such as high-impact (Cazzola et al., 2017; Silvestros et al., 2022), running (Raabe and Chaudhari, 2016), lifting (Beaucage-Gauvreau et al., 2019) or throwing activities (Molinaro et al., 2020) as well as the study of pathologic movement patterns (Overbergh et al., 2020). While current ID models of the spine are becoming increasingly complex, they face the challenge of relying on accurate and comprehensive kinematic data, which most of the current motion capture technologies cannot provide for the spine. Assessment of spinal motion from skin-mounted reflective markers is particularly difficult due to the inevitable movement artefacts prevalent on the spine during dynamic movements (Zemp et al., 2014; Mahallati et al., 2016; Papi et al., 2017). Significant efforts were made to estimate the magnitude of skin movement along the spine, however, have failed to give a clear relationship between a particular movement and the displacement of respective spine markers (Zemp et al., 2014). Moreover, the number of identifiable landmarks through palpation of the back is limited to the spinous process that on its own is not sufficient to provide reliable information on the position and orientation of the corresponding vertebra (Galbusera and Wilke, 2018; Millar et al., 2019).

As a result, for the analysis of motion tasks, the complexity of ID models is commonly drastically reduced by i) assuming interconnecting three DOF spherical joints and ii) lumping the motion of various adjacent vertebrae. The latter is typically realised by implementing linear kinematic coupling constraints (Christophy et al., 2012; Cazzola et al., 2017) further limiting the total DOF of the system, recently reviewed in Alemi et al. (2021). Furthermore, most of the current ID spine models neglect the contribution of passive joint stiffness produced by ligaments and IVDs (Christophy et al., 2012; Bruno et al., 2015; Cazzola et al., 2017; Beaucage-Gauvreau et al., 2019; Molinaro et al., 2020; Silvestros et al., 2022), thereby, highly overestimating the net muscle forces necessary to hold the spine in place.

The major aim of this study was to quantify the effect of modelling individual passive structures (i.e., ligaments and IVDs) on the remaining joint forces and torques carried by muscles in the human spine using the ID approach. For this purpose, the recently developed generic spine model (Hammer et al., 2022; Meszaros-Beller et al., 2023), implemented in the demoa software environment (Schmitt, 2022), was used. The thoracolumbar spine model including six DOF intervertebral joints, a detailed musculature, intersegmental ligaments and IVDs, previously used in FD simulations of a forward flexion-extension movement (Meszaros-Beller et al., 2023), was transferred into the OpenSim MSK modelling platform. Using the full kinematic description obtained from the FD simulations, systematic ID analysis was performed in a step-wise approach increasing the model complexity by adding individual elements and compare the ID analysis results to a standard ID “plain” model neglecting the role of spinal ligaments and IVD stiffness.

A secondary aim of this study was to compare the model performance in two different MSK modelling environments (i.e., demoa and OpenSim). For this reason, solutions for the equivalent modelling of individual structures (i.e., muscles, ligaments and IVDs) in OpenSim were found. Under consideration of identical geometry and soft tissue properties, the ID model kinematics and kinetics were cross-validated against FD simulation results using the in silico derived motion data.

The novelty of this work comprises the capability to use a sophisticated generic spine model (Meszaros-Beller et al., 2023) across two different modelling environments exploiting the strength of each environment and MB approach, e.g., the possibility to remove biological structures under the conservation of movement in the ID approach. Moreover, in a quantitative investigation, this study has shown that neglecting passive spinal structures leads to a significant overestimation of remaining joint forces and torques that muscles must balance.

IK and ID analysis was performed using the generic OpenSim spine model with and without the contribution of individual structures according to Table 2 using in silico motion data for a flexion-extension movement. For the interpretation of ID analysis results, the reader is referred to Supplementary Figure 1.

Current ID models of the spine (Christophy et al., 2012; Bruno et al., 2015; Cazzola et al., 2017; Beaucage-Gauvreau et al., 2019; Molinaro et al., 2020; Silvestros et al., 2022) commonly face two major limitations that is i) the lack of accurate kinematic data and ii) the absence of passive elements, i.e., ligaments and IVDs. While the latter can be solved by modelling the individual structures or by the incorporation of a “lumped” joint stiffness (Wang et al., 2020), the accurate and reliable measurement of spinal motion is made challenging due to skin movement artefacts, marker misplacement (Papi et al., 2017) and the lack in identifiable bony landmarks from the back that is limited to the spinous process (Galbusera and Wilke, 2018). As a result, commonly gross spinal motion is recorded and applied on models with reduced model complexity, e.g., using simplified three DOF spherical joints and linear kinematic coupling constraints for the thoracic and lumbar spinal region (Alemi et al., 2021), respectively. These simplifications, however, are in conflict with the increasing complexity of spine models and introduce intrinsic errors in the computation of kinematics that may propagate to even larger errors in the computation of joint loads and muscle forces (Byrne et al., 2020; Alemi et al., 2021). Currently there is no solution to this problem.

In the present study, in silico motion data obtained from a FD simulation allowed to perform ID analysis of a detailed fully articulated thoracolumbar spine model with six DOF joints avoiding error-prone experimental data. For this purpose, the recently developed generic baseline model (Meszaros-Beller et al., 2023) was transferred with respect to its geometric and soft tissue properties into an established MSK modelling platform, OpenSim 4.3. Note that the FD simulation kinematics might differ from human spine flexion. However, due to the aforementioned limitations, experimental datasets are not suited for this comparative approach yet.

The minimal cross-platform losses in accuracy due to the transformation of rotations from degrees to radians of ± 0.008° were considered negligible. The advantage of using in silico motion data over experimental marker-based motion data typically acquired in a motion capture laboratory is that one has access to precise intervertebral kinematics that is not biased by common limitations of marker-based motion capture techniques or any post processing, i.e., filtering and gap-filling. Three markers per vertebra were successfully tracked over the course of the flexion-extension motion in the FD simulation in demoa and used to run IK analysis on the same model implemented in OpenSim. Given the motion data were obtained from the same model, no marker error was present. In Figure 2, the identical model kinematics were verified via comparison of individual lumbar joint angles between the FD and ID approach.

In addition to the difficulty in obtaining accurate kinematic data described above, the majority of ID spine models neglect passive elements such as ligaments and IVDs that contribute with a posterior torque component to the system dynamics counteracting the predominant anterior loading of the spine (Supplementary Figure 1). As a consequence, neglecting passive elements results in an overestimation of remaining joint loads and, subsequently, predicted net contribution of muscle forces to the joint force and torque. At the same time, commonly employed optimisation processes are known to underestimate muscle co-contraction (Cazzola et al., 2017; Beaucage-Gauvreau et al., 2019).

For the first time, this study assessed the effect of individual passive structures (i.e., ligaments and IVDs) on the ID analysis results using FD-derived in silico kinematics in a step-wise approach increasing the model complexity (Table 2). As shown in Figures 3A–D, the incorporation of ligaments and/or IVDs changed the magnitude and the pattern of the ID results. Thus, it can be concluded that passive elements likely affect predicted muscle recruitment patterns. Typically, this is obtained through static optimisation (SO), an optimisation process following the ID analysis in which the ID-derived generalised net joint forces and torques, or remaining loads after consideration of passive soft tissue elements are distributed to individual muscle fascicles. Even though, SO analysis lied beyond the scope of this study and individual muscle contribution were not examined here, it will be considered in future work. However, the muscle forces used to produce the spine kinematics are presented in Meszaros-Beller et al. (2023) and raw data is available from Hammer et al. (2022).

It is worth noting that the simplest model analysed in this study, i.e., the “plain” model neglecting all passive structures was representative for the standard ID spine models currently available in literature. If only the VB and segment masses were incorporated in the spine model, the corresponding ID results in Figure 3A imply the patently wrong statement that with forward flexion muscles would have to push and pull at the same time in order to withstand gravitational load and the anterior torque generated through the upper body. Without additional decompressive elements, the joints would need to be simplified to three DOF spherical joints inhibiting translational joint movement to avoid this problem, see Supplementary Section 1.

Regarding the IVD and ligament implementation in OpenSim and demoa, ID-derived structural contribution was cross-validated at every lumbar level against FD simulation results to ensure their identical implementation. Both ligaments (Figure 4A) and IVDs (Figure 4B) contributed equally to the net joint axial force F _y and the net joint torque in the sagittal plane T _z in the FD and ID approach over the course of the flexion-extension movement.

With respect to the ID analysis results, ligaments and IVDs both significantly affected the estimated remaining joint loading (Figures 3B–D): Between L1/2 and L4/5 the modelling of passive structures reduced both net contribution of muscle forces, compressive loading and anterior torque, on average by −200% and −75%, respectively. The relatively high joint decompression through passive elements can be partly attributed to the IVD offset force F IVD,0 that was responsible for a −101% reduction in compressive loading entirely compensating for the gravitational load as intended. These results demonstrate that passive structures contribute to a redistribution of spinal loading among the different structures and should not be omitted in MSK spine models.

As anticipated, the implementation of ligaments and IVDs, respectively, had an opposite effect on the axial loading which can be explained through their line of action: while IVDs are ‘pushing’ vertebrae away from each other under compression, ligaments act in the opposite direction and restrain their ‘sliding apart’ through tensile forces (see Supplementary Figure 1). Thereby ligaments act in line with the muscles. This is in accordance with the results presented.

Furthermore, the presented net ligament contribution to the axial force reveals the strong compressive force onto the IVDs which would be underestimated in a spine model without ligaments included. Thus, considering ligaments does not only affect calculated remaining joint loads and corresponding net muscle forces in the ID analysis but already the predicted bone-on-bone force (Winter 2009) in FD simulations. In our case, this force acting onto the cartilage between bony segments comprises only the IVD loads. Its estimation has the highest clinical relevance of all modelled soft tissue elements in this study, e.g., for the design and development of biomaterials for IVD replacement.

The validity of the “all elements” spine model presented in this study was previously discussed (Meszaros-Beller et al., 2023, Section 4.1). Even though the structural complexity of the “all elements” spine model is one rarely seen among ID spine models, limitations include the absence of the intra-abdominal pressure, the limited modelling of the ribcage, and the rudimentary modelling of the facet joints (Meszaros-Beller et al., 2023, Section 4.1). In the future, the presented model can be used for various comparative studies, e.g., exploring different neural control solutions and muscle recruitment patterns employed in the FD and ID approach. Further, by having access to precise kinematics the inaccuracies introduced through kinematic coupling constraints (Alemi et al., 2021) could be evaluated. The work presented in this study cross-validated the kinematics and kinetics of the developed generic spine model and its implementation in OpenSim providing the necessary baseline for further investigation.

Summarising, the presented approach has demonstrated the capability of cross-platform analysis of an identical detailed generic spine model by using the FD and ID approach avoiding common limitations in motion capture. Moreover, the used step-wise approach allowed us to quantitatively investigate the effect of passive structures on the ID results. IVDs and ligaments together have shown to significantly reduce remaining compressive joint loading implicating reduced muscle forces needed to balance the net joint forces and torques. This work provides the necessary baseline towards exploring different neural control solutions employed in the FD and ID approach.