The anatomy and mechanical properties of the underlying passive hybrid FE-MB LSS model were the same as reported previously (Remus et al., 2021). It consisted of five rigid lumbar vertebrae L1–L5, the rigid sacrum S1, the fiber-reinforced FE intervertebral discs, the facet joints, and the pre-tensioned ligaments. The discs were composed of the quasi-incompressible nucleus pulposus and the surrounding annulus fibrosus with five crisscrossed collagen fiber rings. Hyperelastic material models (Yeoh and Mooney-Rivlin) were used to describe the complex nonlinear stress-strain behavior. The entire passive model is available on GitHub (https://github.com/RemusR9/artisynth_lumbosacralSpineModel) and Zenodo (https://doi.org/10.5281/zenodo.4453702). The anatomical basis was the male Visible Human Project (VHP) (Spitzer et al., 1996) which was adapted in the posture of the LSS (Roussouly et al., 2005). The LSS was extended to include the thorax, humeri, and pelvis. Because the focus was on the LSS, the thoracic region was represented as a single lumped rigid body (Arjmand and Shirazi-Adl, 2006a; Ignasiak et al., 2016), consisting of the vertebrae C7 to T12, the ribs, and the superior segments of the humerus. Thorax and L1 were rigidly connected. With respect to lordosis and a sacral angle of 36.3°, the pelvic alignment was adjusted according to Le Huec et al. (2019) to 45.36° angle of pelvic incidence and 9.06° pelvic tilt to allow sagittal balance in upright standing. The resulting actual tilt of the iliac crest was thus about 2.2° and classifies as an active tight posture balanced between back and abdominal muscles (Schünke et al., 2018). Pelvis and sacrum were rigidly connected because influences of sacroiliac joint on LSS biomechanics for adults are negligible (Le Huec et al., 2019). Upper body segment parameters (Vette et al., 2011) of the male VHP, including center of mass (COM) and moments of inertia, were transformed into the present global coordinate system. The joint centers of the lumbar and thoracic spine serve as points of reference for this. The lumped segment parameters were applied to each rigid model component (Supplementary Material S1). Parameters cranial to L1 were concentrated in COM_UB and applied to the thorax. The entire model was symmetrical about the mid-sagittal plane.

The active musculoskeletal model extended the passive hybrid FE-MB LSS model to include the muscles of the lower back and abdomen. Our workflow was based on a semi-automated procedure for generating subject-specific musculoskeletal MB models of lower extremities (Ascani et al., 2015; Modenese et al., 2018). We adapted the procedure to efficiently and reproducibly transfer reference model muscle data to a target model. The input data used in the following (Figure 2) are provided in Supplementary Materials S2–S4.

For this, we initially created three skeletal models in NMSBuilder: 1) the modified skeleton of the VHP as the target anatomy (according to the passive anatomical model described in 2.1), plus 2) the OpenSim skeleton “Lumbar_C_238” from Christophy et al. (2012) and 3) the “Twente Spine Model” skeleton by Bayoglu et al. (2017) both as reference models. To describe bone geometries used to register a cloud of muscle points, supervised virtual palpations of anatomical landmarks (LMs) were performed by three operators experienced in back anatomy (Figure 2A). To identify and mark the anatomical characteristics with LMs (Van Sint Jan, 2007), the palpation dictionary (Supplementary Material S2) with anatomical descriptions of the LM positions (created in advance) was available to each operator after a technical introduction. There was the option to use an illustrated documentation of the palpated target anatomy during the initial palpation process (Supplementary Material S3). The registration atlas (Supplementary Material S4) facilitates the workflow in NMSBuilder by providing all LM names per bone in sorted order. All LMs of a bone were aggregated as a LM cloud. To verify repeatability, LM positions were examined for their mean squared distances from the respective mean for each model separately (Figure 2B). A consistent set of anatomical LM clouds was selected for each model for the registration of the reference muscle points (Figure 2C).

The implemented main muscle groups (Figure 3) of the lower back (Christophy et al., 2012) included latissimus dorsi (LD), quadratus lumborum (QL), multifidus (MF), obliquus internus abdominis (IO), obliquus externus abdominis (EO), psoas major (PM), rectus abdominis (RA), iliocostalis thoracis (IT), iliocostalis lumborum (IL), longissimus thoracius (LT), and longissimus lumborum (LL). To complete the abdominal muscles compressing an abdominal cavity (Cresswell et al., 1992; Stokes and Gardner-Morse, 1999; Daggfeldt and Thorstensson, 2003; Essendrop and Schibye, 2004), we included the transversus abdominis (TA) from Bayoglu et al. (2017). Muscles can be divided into functional, one-dimensional fascicles and were modeled in three different forms to best represent muscle paths (de Zee et al., 2007): as a straight line directly connecting insertion and origin point, redirected by means of via-points, and non-linear wrapped around auxiliary bodies. Reference marker clouds consisted of the respective muscle points with reference LM clouds for each bone. Once combined and imported to the target model, an affine transformation was used to register a reference marker cloud on the target model. Raw specific target muscle point clouds were created for each bone. Because anthropometric differences of the bones remained despite the registration, not all muscle points touched the surfaces of the target bones (Figure 2D). Excluding via points, the muscle points were shifted to the bone surface with the least distance using a customized Matlab script (Modenese et al., 2018). All muscles of the right side of the body were generated from the automatically corrected muscle points in NMSBuilder.

Wrapped muscles were redirected by three geometric auxiliary bodies (Table 1; Figure 1). Left and right wrapping body were rigidly attached to the thorax and redirected the long posterior muscles. The lumbar wrapping body was connected to the sacrum via a linked hinge and slider joint. Both joint centers were located in the L5 superior articular process (x = −0.02 m, z = 0.0345 m). The degrees of freedom of the kinematic joint chain 1) rotation about the joint y-axis and 2) translation along the slider joint rotating with the hinge joint, were controlled in extension by the position changes of vertebrae L1 and L2: 1) The lumbar wrapping body was rotated such that the angular sum of the vectors S1-L1 and S1-lumbar wrapping body to the perpendicular on S1 remained constant. 2) Anterior translation with a ratio of 0.54–1, when the z-distance between L2 and S1 reduced. In flexion, only the positional change of L4 was used for the sole translational displacement anteriorly along 2) with a ratio of 1.3 to 1.

For musculoskeletal finalization, the muscle paths and attachment points were manually verified with the VHP image data (Figure 2E) and anatomical descriptions (Hansen et al., 2006; Bogduk, 2012; Adams et al., 2013; Schünke et al., 2018). The generated muscle paths resulted directly from the muscle points transferred to the target anatomy. Minor adjustments were made by relocating muscle points on the bone surfaces, considering muscle wrapping. The final musculoskeletal target model created in NMSBuilder included all bones, the right-sided muscles (Table 2, Supplementary Material S5), and the four auxiliary bodies (Table 1). The right-side musculature comprised 129 (119 without TA) muscle fascicles. For the transfer to ArtiSynth, an OpenSim (Delp et al., 2007) model file (.osim) was exported. The ArtiSynth “OpenSimParser” enabled the import and subsequent integration of the musculoskeletal data. The left-side musculature was mirrored on the sagittal plane.

Each muscle fascicle of the 12 implemented muscle groups was described as a tension-only force-generating spring-damper system and modeled via a Hill-type muscle model (Zajac, 1989). Muscle parameters were specified using the “Millard2012AxialMuscle” material, with tendons assumed to have no compliance and non-linear normalized force-length as well as force-velocity curves by Millard et al. (2013). The reference parameters for physiological cross-sectional area (PCSA), fiber to tendon length, sarcomere length, optimal fiber length, pennation angle, and tendon slack length were taken from the respective sources (Christophy et al., 2012; Bayoglu et al., 2017). Multiplying the reference PCSA with the model specific muscle tension K (cf. section 2.7) gave the maximum isometric muscle force F 0 r e f M . Because the muscle volume is proportional to the body mass b and that the muscle-fiber length is proportional to the musculotendon length l M T , the maximum isometric reference muscle force F 0 r e f M was scaled according to Eq. 1 (Correa and Pandy, 2011). F 0 M = F 0 r e f M · b b r e f · l r e f M T l M T (1)

Musculotendon lengths were extracted from the three fully built musculoskeletal models (cf. section 2.2). A sarcomere length of 2.8 µm and a tendon slack length when the muscle is in neutral position was assumed (Christophy et al., 2012). The muscle modeling parameters and muscle attachments used are provided in Supplementary Material S5.

To define the movements of the abdominal muscles attached to the elliptical abdominal plate (AP), we implemented an open kinematic chain of three linked joints (Figure 4A). The kinematic chain connected AP to sacrum and was composed of three joints: 1) universal joint (UJ), 2) hinge joint (HJ), and 3) cylindrical joint (CJ). Two massless auxiliary bodies in the centers of HJ and CJ (Table 3) were used to connect the joints. The origin of UJ was located in the center of vertebral body L5. The complete kinematic chain was: Sacrum–UJ–Auxiliary body 1–HJ–Auxiliary body 2–CJ–AP. For full kinematic control of the five degrees of freedom, the joint coordinates θ ₁ , θ ₂ , θ ₃ , φ, and Z, visualized in Figure 4B, were determined by the angular changes α _j about the three principal axes j of the thorax with respect to the pelvis and the constantly held volume of the abdominal cavity. This kinematic correlation was implemented via a polynomial function according to Eq. 2. C C   α j = x , y , z = p 1 α j + p 2 α j 2 (2)

The polynomial coefficients p _i (Table 3) were determined using the procedure visualized in Figure 5 plus downstream optimization by combining in vivo measurements and simulations. For this purpose, motion capture recordings (Vicon Vero v2.2 system with 12 cameras) of four healthy, trim men (30.75 ± 1.5 years, BMI 22.67 ± 1.62) were made. Ethic approval was obtained from the Ethics Committee of the Medical Faculty of the Ruhr-University Bochum (23-7801 04/20/23) and written informed consent was obtained from all participants prior to inclusion in the study. Ten reflective markers with a diameter of 14 mm were securely attached to the volunteers’ skin with fixing tape directly or interconnected with a plastic triangle, as shown in Figure 5A. At the same locations virtual marker pendants were added to the simulation model. Starting in an upright position, the participants performed three upper body movements in flexion-extension, axial rotation, and lateral bending. Each to their maximum voluntary range of motion. Participants were instructed to move slowly, at their own pace, and to move primarily out of their lower back. Before starting the measurements, participants undertook at least one practice movement in each direction to ensure smooth recording. During the recordings, attention was paid to ensure that the pelvis tilted little, there was no breathing during the movements, and the trunk muscles were slightly tensed. To process the raw measurement data, individual 3D marker motion curves were extracted starting from upright standing (Figures 5B, C), trimmed based on the maximum thorax range of motion (ROM) (flexion = 30°, extension = −20°, lateral bending = 15° and axial rotation = 10°), and projected to the respective principal motion plane. In order to obtain one valid reference curve for each direction of motion, the individual curves of the marker LAL41 were combined for each coordinate with the modified dynamic time-warping (DTW) method (Wang and Gasser, 1997; Bender and Bergmann, 2012). Lines of best fit with the matching polynomial interpolation degree for each coordinate are listed in Supplementary Material S1 and visualized for flexion in Figure 5D.

Kinematic similarity between AP and in vivo data was evaluated by comparing the components of the motion curves of the marker LAL41. To find the ten polynomial coefficients p_i (Table 3), that provide the highest kinematic compliance, five multicriteria optimizations (patternsearch from the Matlab 2020b optimization toolbox) with global settings were performed. One in each direction of motion independently to determine initial values for p_i for the final optimization in all directions combined. Target values to be minimized were the unnormalized distances from modified DTW (Bender and Bergmann, 2012) between both LAL41 coordinate components of polynomial curves and simulated virtual marker trajectories. In all cases the unnormalized distances were summed up unweighted via a cost function.

To model a resulting IAP from the abdominal muscle contraction, a second dynamically active and quasi-massless, abdominal plate (AP_IAP) was integrated at the same location as AP. This was required because only dynamically active bodies calculate forces beyond their motion. As visualized in Figure 6A, EO, IO, and TA were applied to AP_IAP. RA was redirected via AP without generating a direct influence on the IAP (Essendrop and Schibye, 2004). AP and AP_IAP were connected by a “FrameSpring” (Figure 6B), which is a six dimensional spring that generates restoring forces and moments between rigid bodies. Their linear stiffness in F_AP direction was set to 90 kN/m to ensure a stable calculation and to limit the relative AP displacement to 1.5 cm at maximum tension of all abdominal muscles. For small relative movements in the other five degrees of freedom, 1 MN/m and 500 Nm/rad were used. Diaphragm and pelvic floor were assumed to be rigid (Stokes et al., 2011). Considering a constant volume (V _AC = const.), the abdominal cavity was calculated as a cylinder with height h _AC and diameter d _AC (Figure 6C). The force F _IAP resulting from the IAP was applied to the thorax. Acting cranially, the application point of F _IAP was the center of the diaphragmatic surface A _D , 5.1 cm anterior T12. In order to always act perpendicularly on the diaphragm, the force vector changed its direction with the thorax. The abdominal muscles only partially enclosed the simplifying cylinder barrel and were redirected laterally in such a way that the resulting forces were not considered. Under the geometric assumption of an effective muscle entwinement of c _w = 60%, the relation between F _AP and F _IAP is given in Eq. 3. F I A P = 1 c w · d A C 4   h A C F A P = 1 4   π d A C 2 p I A P (3)

An IAP of 4 mmHg (Andersson et al., 1977; Schultz et al., 1982; Nachemson et al., 1986; Mueller et al., 1998) was applied as offset in the unloaded state and the maximum IAP was limited to 200 mmHg (Essendrop, 2003).

To solve the muscle redundancy problem and mimic a physiological muscle recruitment, the tracking controller (TC) integrated in ArtiSynth (Stavness et al., 2010) was used throughout the simulation time. The tracking-based inverse controller provided an optimized solution to the forward dynamics simulation by finding a set of muscle activations at each time step that drives the hybrid FE-MB model along with other constraints through a target movement trajectory (Stavness, 2010; Stavness et al., 2012). Erdemir et al. (2007) classified this approach as “forward dynamics assisted data tracking”. Thus, all postures were generated purely muscle-actuated, without prescribing complete kinematics to the dynamic bones. For the TC motion target term, we specified five target frames for the thorax and the vertebral bodies L2 to L5 (Figure 7B). Its rotational components were used as motion target values. To neglect spinal compression, but constrain displacements in sagittal plane, the translation in x-direction was additionally considered for the thorax target frame (Center of rotation: x = −0.015 m, z = 0.1 m). While the thorax target frame determined the posture of the model, the vertebral body target frames had the subordinate function of reducing oscillations and stabilizing the LSS. Thus, the thorax target values were weighted 15 times higher in the TC cost function. Two additional terms were added to the cost function: The sum of the muscle excitation group activations squared (Stavness et al., 2010) and the rate of muscle excitation changes. The rate of change was added to avoid muscle activation spikes and associated force peaks on the FE discs, which may lead to instabilities. With respect to the thorax target values, the squared excitation term was weighted with 0.667 and the damping term with 0.3×10⁻⁴. To reduce the number of muscles that were controlled individually by the TC, all muscle exciters were allocated to 12 activation groups for the left and right side, respectively (Table 2).

Aim of our first calibration was to determine a maximum specific muscle tension K by evaluating the maximum isometric back extension that the model could exert. Because in vivo measurements demonstrated that the trunk musculature of healthy men generate maximum isometric back-extension torques of 210 Nm in −10° extension and 260 Nm in 10° flexion about the L5/S1 level (Daggfeldt and Thorstensson, 2003) we applied an increasing force equivalent to the thorax at T7: 700 N in extension ( α y * = −10°), 870 N in flexion ( α y * = 10°). The force rotated with the thorax and always acted perpendicularly on it. Gravity was deactivated to mimic the upper body lying on its side. According to Eq. 4 we evaluated equilibrium via the tracking error ∆ α y as the rotational deviation between the thorax target frame and the thorax. ∆ α y = | α y * − α y | (4)

Regarding the upright posture after settling (see section 2.9 Validation), α y * defines the rotation of the thorax target frame and α y the predicted rotation of the dynamic thorax in sagittal plane. Reported values for K vary between 10 and 100 N/cm² (Daggfeldt and Thorstensson, 2003; Hansen et al., 2006). In models of the lower back, most commonly values of 46 (Bogduk et al., 1992; Stokes and Gardner-Morse, 1995; Christophy et al., 2012), 60 (Arjmand and Shirazi-Adl, 2006a), 90 (Arshad et al., 2016), or 100 N/cm² (Bruno et al., 2015; Bayoglu et al., 2019; Favier et al., 2021a; Lerchl et al., 2022; Malakoutian et al., 2022) were assumed. Thus, we tested values of 46, 73, and 100 N/cm² for K, expected ∆ α y ≤ 1.0°, and that no muscle was fully activated by the TC.

For 46 N/cm², equilibrium could not be established in any posture. At 73 N/cm² equilibrium was achieved only in flexion, with a tracking error of 0.85°. With K = 100 N/cm², ∆ α y was 0.8° in extension and 0.55° in flexion. Furthermore, the model could generate a maximum extension moment of 265 Nm and a maximum flexion moment of 430 Nm for ∆ α y = 1.0°. This showed a much higher trunk stability in flexion, only considering the tracking error. Because differences in the L4/5 IDP between 73 and 100 N/cm² in flexion were less than 3.0%, we assumed K = 100 N/cm² for all muscles. In upright posture, when the abdominal muscles attached to AP_IAP without TA were activated to 80% or with TA to 65%, the maximum IAP was generated.

In the second model calibration, segmental rotation contributions for the vertebral target frames (Figure 7B) were determined with respect to their influence on the biomechanical model responses. Consequently, the rotational constraints of the vertebral target frames were part of the optimization problem solved by the TC and influenced the model responses. The deviations between the rotation specifications and the resulting intervertebral rotations (IVRs) of the vertebrae L2 to L5 as well as the calculated IDP L4/5 were used for evaluation. Often referred to as spinal rhythm, in vivo studies have shown that the segmental rotation contributions to the overall lumbar motion are level specific and non-linear over the entire motion (Breen et al., 2021). Based on literature sources we tested seven different fixed segmental target frame rotation contributions (Figure 8A):

R1) Proportionally equal contributions of 20% each

In cases where no segmental contribution was measured for L1-L2, we assumed 80%–100% of the contribution of L2-L3, following complete in vivo measurements, and adjusted the sum of all five contributions to 100%. Therefore, and because of the rigid connection between T12 and L1 limiting the mobility of L1, we evaluated only the results of L2-S1.

Segmental target frame rotation contributions R1 to R7 and the resulting IVRs are visualized in Figure 8A. The rotational deviations between vertebral target frames and vertebrae were smallest for R4 and R5. Considering all three postures, the root-mean-square errors for R4 and R5 were 9.89% and 9.34%, respectively. Consistent with Arshad et al. (2016), the biomechanical influences of spinal rhythms increased with greater flexion. Directly correlated with this were the IDPs at all levels. Comparison of the percent IDP changes on L4/5 with in vivo measurements from Wilke et al. (2001) showed good agreement for all cases except R2, R3, and R6 (Figure 8B). In upright posture, IDP spread was less than 2.8% and IVR deviations were below 0.3°. The choice of segmental contribution was almost irrelevant in phases of low thoracic tilt (−5° < α y * < 5°). Similar to the experimental setup of Breen et al. (2021), the pelvis in our model was stationary. Furthermore, because segmental target frame rotation contributions R4 and R5 provided the most suitable boundary conditions, we used R5 with the lower proportion of L1-L2 for model validation: 23.4% at L1-L2, 29.3% at L2-L3, 25.9% at L3-L4, 16.0% at L4-L5, and 5.4% at L5-S1.

For validation, the calculated biomechanical model responses of 13 LCs (Table 4) were compared with in vivo studies (Nachemson, 1965; Schultz et al., 1982; Sato et al., 1999; Wilke et al., 2001; Arjmand and Shirazi-Adl, 2006a; Takahashi et al., 2006; Wong et al., 2006). These LCs were selected due to the availability of the L3/4 or L4/5 IDP measurements. The LCs comprised sagittal postures of the thorax target frame graded in 10° in the range −10° ≤ α y * ≤ 30° (Figures 7C, D), each with or without the external load F _Load . For example, when compared to the measurements from Wilke et al. (2001), the external load approximated holding a 20 kg crate in both hands with arms bent or extended in front of the upper body. In Takahashi et al. (2006) 5 kg were attached to both wrists of the vertically hanging outstretched arms. The force F _Load , always acting vertically, moved with respect to the shoulders and their distance from the origin was l _Load (Figure 7A). Dead weight displacements of the arms were neglected. Origin and sternum had a distance of 0.133 m in neutral posture. Biomechanical model outputs used for validation were the IDPs of the three caudal lumbosacral levels, the predicted muscle forces, and the IAP. Each LC was simulated with muscle group TA (w/TA) and without TA (w/oTA). This is due to the combination of two independent data sets for the abdominal muscles and the verification of the sensitivity of the biomechanical model responses.

All simulations started in a neutral state without gravity and without external forces with the LSS and thorax in upright position (Figure 7A). To increase computational stability, the resting muscle tone was set to 0.1%, which affected IDPs by less than 1% compared to no muscle tone. The sacrum was stationary, maximum step size was set to 0.01 s, and a first order backward integrator (“ConstrainedBackwardEuler”) was used. Within the first simulation steps, the gravity was ramped up to 9.81 m/s² in the negative z-direction. This was followed by a settling of the model, done by a thorax target frame rotation of α y * = 2° and the return to the upright posture. This small reversible movement of all dynamic bones established an energetically favorable condition with minimal IDPs between passive system stiffness, dead weight, and muscle forces. We defined the resulting orientation of the vertebrae as the stable, upright reference posture (α _y = 0°). Due to the rigid connection of thorax and L1, α _y also represented the absolute rotation of L1. Subsequently, the poses were set by moving the target frames steadily, optionally followed by an increase of F _Load in form of a smoothstep function. In Supplementary Material S6, the simulation procedure is visualized as an example for LC F30₃. Results were evaluated when the entire model reached a static state (equilibrium).

The passive LSS without surrounding active musculature is an unstable system that tends to buckle even under low vertical loading (Patwardhan et al., 1999). The co-contraction of abdominal and lumbar muscle groups combined with the activation of the diaphragm and the pressure that builds up provides an mechanism for stabilization (Hodges et al., 2005; Kuo et al., 2021). Spinal instability, clinically defined as loss of the spine’s ability to maintain its displacement patterns during physiological loads, is an important cause of low back pain (Panjabi, 2003). Destabilizing factors may include force reducing and response degrading dysfunction of MF and TA (Kuo et al., 2021). Persistently increased abdominal muscle tension, in turn, leads to abdominal hypertension (Tayebi et al., 2021) that can lead to health problems such as pain or organ deterioration (van Ramshorst et al., 2011). The physiological interplay of all components is therefore of particular relevance for spinal stability. Within the scope of our musculoskeletal modeling, however, respective interactions in the lower trunk could often only be represented in a simplified way, which is discussed in the following.

The active model built consisted of a trunk musculature in addition to the passive FE-MB LSS. Due to the unavoidable combination of data for anatomy and materials as well as calibration and validation, our model represented an average healthy and middle-aged man without individual characteristics. Consequently, patient-specific conclusions were limited. However, to achieve better comparability of results and to benefit from the wide availability of data, dimensions and masses were based on the Male VHP. For the musculature, we combined morphology and material reference data from two independent sources which thus included the relevant parts of the autochthonous back muscles as well as the abdominal muscles (Table 2). Reason for the combination was our consideration of the IAP (Essendrop and Schibye, 2004) and the omitted TA in the primarily used OpenSim model of Christophy et al. (2012). Limitations of the second data set by Bayoglu et al. (2017) were the fusion of L5 with the sacrum, whereby the authors could not exclude alterations of the measured muscle groups. Because intertransversarii and interspinalis muscles are mainly considered as proprioceptive sensors rather than force actuators (de Zee et al., 2007), we did not supplement them with additional reference data. In addition, compared to other recent models (Lerchl et al., 2022), fascicles not included in the model of Christophy et al., such as MF attachments to the thoracic spine, were not manually added. To integrate the muscle morphologies into our model, a published and established semi-automated codified registration procedure (Ascani et al., 2015; Modenese et al., 2018) was modified and utilized. A challenge in our case was the considerably higher number of muscle groups and fascicles as well as their often non-linear course. Likewise, the degree of morphological detail of the reference bones was different. The sacrum of the OpenSim model, for example, had a very low level of detail, which is why the repeated palpations showed high local deviations. Absolute highest LM variations were found for the largest element to be palpated, the thorax at angulus costae of the 7th and 9th rib and the 7th and 11th corpus costae (see palpation dictionary in Supplementary Material S2). LMs for vertebrae and pelvis could be determined for all data with the highest repeatability. In general, depending on the 3D bone geometry quality and morphology, LMs could be determined with varying accuracy. Compared to the VHP image data imported into NMSBuilder, the paths of the QL fascicles required the most manual correction. After registration, these passed almost completely through the transverse cross-section of the E.S. at level L4. Little to no correction was required for the remaining muscle groups. Overall, the use of the adapted semi-automated registration procedure (Figure 2) had proven to be reliable and efficient. Furthermore, the workflow may become relevant for easy modification of the musculature or for transferring the muscle morphology used here to deviating bone geometries. Also, patient-specific musculoskeletal models can be generated in this way based on one reference model (Modenese et al., 2018).

Four auxiliary rigid bodies were integrated into the model for redirecting or attaching muscle fascicles (Table 1). Of these, the lumbar wrapping body and the AP were kinematically controlled via open kinematic chains. Advantage of using linked joints were the inherent kinematic determinacy with easy manipulability of the chain as well as using scaling factors during investigations and optimizations. The kinematics of the lumbar wrapping body were implemented according to the smaller cylinder integrated in the OpenSim model so that the lumbosacral muscles were kept near the lordotic spine in extension. We did not integrate an equivalent of the larger of the two lumbar wrapping cylinders into the presented model. The smaller cylinder redirected all muscle fascicles for the LCs considered.

Consistent with the state of the art in MB models, the AP replaced all ligamentous structures of the central abdominal wall (Schünke et al., 2018) and defined the anterior muscle attachment points and abdominal muscle movements (de Zee et al., 2007; Christophy et al., 2012; Malakoutian et al., 2018; Liu et al., 2019b). Kinematics optimization of the AP was based on motion capture measurements performed only on a small group of young men. It must therefore be assumed that no universally valid abdominal kinematics was implemented. However, the interindividual high concordance of the measured marker trajectories suggested to us a suitable approximation for this target group. Not all recorded marker movements of our in vivo study were used for the presented optimization. Investigations by us, not described further here, showed that consideration of further measured values did not add value to the derivation of abdominal kinematics. Other measured values included the pose of the abdominally applied triangle and the movements of marker SPP4 (Figure 5). Their consideration can be useful, for example, in more complex movements beyond the three main planes or in pathologies. A limitation might be, however, that we based abdominal kinematics solely on the position and orientation of the pelvis and the thorax. As already realized via SPP4, further vertebral body displacements should be considered in the future and thus the absolute distance between lumbar spine and abdominal wall should be taken into account. In addition, it should be noted that abdominal and upper body movements can be measured easily and noninvasively using markers, but both palpation of bony LMs and displacement of the skin can lead to erroneous results.

To consider IAP in simulation models of the lower back, various approaches have been implemented (Han et al., 2012; Arshad et al., 2017; Malakoutian et al., 2018; Liu et al., 2019b). The use of a kinematic chain with two abdominal plates to directly utilize compressively acting abdominal muscle forces from IO, EO, and TA represented a new approach. Since AP and AP_IAP had the same geometry and were in the same pose in the model, we counted only AP for the auxiliary bodies. The mechanical effect of the IAP was modeled as a cranially aligned force acting perpendicularly on the diaphragm that rotated with the thorax. To convert between the total force acting on the AP_IAP and the IAP, we assumed a constant volume of the abdominal cavity idealized as a cylinder. Both the upper body posture and the correlating AP position therefore had an influence on the IAP. The computed volume change via the IAP thus affected the stability of the system and was part of the overall muscle activation problem. The entwinement factor implemented for scaling was based on the geometrical assumption and the redirection of the TA fascicles by via points. Resulting compression forces acting on the via points were not considered. This resulted in considering only the part of the muscle forces acting perpendicularly on the AP_IAP. Especially because the factor c_w is variable, its variation should be systematically investigated in further studies. Instead of the direct muscle force measurement, the implementation of a fully analytical approach is alternatively conceivable, as in a fundamental work (McGill and Norman, 1987). For the estimation of the IAP, the dynamic AP_IAP could be omitted, and redirections of TA would become irrelevant. Only the activity of the muscles of the abdominal wall and the present geometric conditions would be input variables. For a discussion of the difficulties involved, we refer the reader to the work of McGill and Norman (1987). Another limitation of our approach to be considered were the equidistant distances of the abdominal muscles in all postures as a result of the anterior attachments to the rigid AP_IAP. Potential effects on pressure as well as activation resulted in flexion could not be assessed in more detail at this point. Thus, contraction of the trunk muscles leading to deformations of the abdominal cavity and the associated elevation of the anterior abdominal wall (Todros et al., 2020) could not be studied either. In addition, we neglected the pressure-increasing influences on the abdominal cavity resulting from diaphragm contraction, which may contribute to spinal stabilization (Hodges et al., 2001) as well as the influence of breathing to the control of abdominal muscles (Hodges et al., 2015).

The hybrid modeling approach included five fibre-reinforced FE discs and ten FE inferior articular facets in addition to the rigid bones, muscles, and ligaments. Compared to MB models, this led to several downsides. The modeling effort was considerably increased and is close to FE models (Affolter et al., 2020; Pickering et al., 2021). Manual segmentation and meshing were required, which, along with the level of intervertebral disc detail, is currently a bottleneck for anatomical modification. With the initial aim of directly calculating the IDP and predicting the complex non-linear kinematic responses of the LSS, we conducted a sensitivity study (Remus et al., 2021) to reduce the complexity of the disc. More detailed analysis of the discs is possible, but would require a finer mesh and more parameters, increasing the complexity of the model. The use of FE meshes for discs in active hybrid models could therefore allow the influence of motion and loading on disc degeneration to be studied and, for example, fatigue damage to be predicted (Subramani et al., 2020). Partial or complete replacement of rigid vertebrae with FE bodies is also possible. In the future, this may allow the investigation of structural-mechanical issues under in vivo similar loads, for example, to study bone adaptations due to mechanical stimuli (Smit et al., 1997; Lipphaus and Witzel, 2019; Favier et al., 2021b), and the osseous integration of spinal cages (Knapik et al., 2012; Malakoutian et al., 2015; Abbasi-Ghiri et al., 2022). With regard to the already integrated FE bodies and contacts, the computing time has to be mentioned. Using a desktop PC with Intel i7-10700K @ 3.80 GHz, 32 GB Ram and 1 TB SSD running Windows 11 Pro 64-bit, LC F30₃ took an average of 30 min to calculate. Comparable load cases are calculated in ArtiSynth in several seconds by a purely forward dynamics, muscle-driven MB model (Malakoutian et al., 2018).

We calibrated the maximum specific muscle tension and segmental target frame rotation contributions because corresponding literature data from experiments and simulations showed wide variations. The iterative procedure within the calibration, required to ensure consistency of results, was not discussed in this paper. Each simulation started in the same neutral state without external forces, but with ligaments pre-tensioned and the TC activated. This was the state in which the model was geometrically built. During the first simulation steps, gravity was ramped, followed by the settling phase in which the bones were able to move into an energetically favorable state. This was necessary because due to the forward dynamics method (Pandy, 2001) bone poses were not kinematically controlled and the upright stable is unequal to the initial neutral state. In contrast to inverse dynamic models (Christophy et al., 2012; Ignasiak et al., 2016; Khoddam-Khorasani et al., 2018; Liu et al., 2018; Favier et al., 2021a; Rajaee et al., 2021; Lerchl et al., 2022), the kinematic inputs were not complete and only used to move the five target frames (Figure 7B). The model dynamics described how the components advanced in time from one state to another. The exact motions of the active hybrid model, which include the IVRs, were not known in advance. By controlling the muscles using the TC, deviations between all target frames and respective bones were minimized for the selected coordinates. For the four vertebral target frames, these were only their rotations in sagittal plane. Resulting muscular forces directly moved the entire spine into various postures and provided spinal stability (muscle-driven approach).

All LCs studied could be generated by the postural specification of the thorax target frame alone. However, as a sole boundary condition, we found that for some LCs with α _y ≥ 20, no stable state for the lumbar vertebrae L2-L5 could be obtained by the TC. As a result, vertebral oscillations were not reduced by the muscles, which in the case of LC F30₂ resulted in inverted finite elements due to excessive distortions of the discs. Consequently, the specifications of the vertebral target frame rotations primarily served to stabilize the spine with the premise of keeping the influence on the model dynamics minimal (low weighting and calibration of segmental target frame rotation contribution). In upright position, sensitivity to IDPs were negligible (ΔIDP < 0.03 MPa) because IVRs and thus muscle activation patterns varied little. In flexion and extension, the influence on the calculated IDPs and IVRs were evident for the studied segmental rotation contributions of the target frames. As measured in vivo (Breen et al., 2021) non-linear segmental rotation contributions should be considered in the future given the relevance of our findings. Ligament strains, IDPs, or facet joint contact forces may represent alternative stability criteria besides the thorax target frame.

The validation of realistic simulation models requires experimental confirmation. In vivo, however, such data are difficult to acquire or may never be obtained. Given the large number of components and the inevitable combination of data from multiple sources, an absolute overall validation of our active hybrid model was not possible at this stage. Sole approach was the comparison of particular biomechanical model responses with experimental results under defined boundary conditions. Based on available publications, we therefore considered absolute IDP and IAP values and relative changes in muscle forces in 13 LCs. Moreover, to increase model validity, we extensively validated the passive LSS with experimental in vitro studies in all spatial directions in advance (Remus et al., 2021). This included ROMs, IVRs, IDPs, facet joint contact forces, instantaneous centers of rotation, functional spinal unit stiffnesses, and intervertebral disc bulges. This was done under the assumption of component validation (Lewandowski, 1982), based on the fact that models consisting of well-validated submodels are likely to be valid. Advantages compared to a direct validation of the entire active model included the higher number of useable in vitro studies, simplified loading conditions, and the influence analysis of LSS structures in case of a stepwise reduction of the anatomy (Heuer et al., 2007). Even though in vivo loading conditions cannot be correctly mimicked using a moment and follower load (Khoddam-Khorasani et al., 2018), muscle activation patterns represent boundary conditions too complex for fundamental, and detailed validation of the basic biomechanical responses of passive LSSs. The validation in this study was limited to symmetrical LCs in the sagittal plane. As other simulation models of the LSS (Liu et al., 2018; Malakoutian et al., 2018) have mostly been validated with the same in vivo studies as we have used, there are no significant differences in the biomechanical model responses. For more comprehensive model validation, additional movements such as axial rotations and lateral flexions as well as non-symmetric loads like holding a weight laterally and combined movements outside the anatomical planes need to be investigated.