Muscle-driven forward dynamic active hybrid model of the lumbosacral spine: combined FEM and multibody simulation

Most spine models belong to either the musculoskeletal multibody (MB) or finite element (FE) method. Recently, coupling of MB and FE models has increasingly been used to combine advantages of both methods. Active hybrid FE-MB models, still rarely used in spine research, avoid the interface and convergence problems associated with model coupling. They provide the inherent ability to account for the full interplay of passive and active mechanisms for spinal stability. In this paper, we developed and validated a novel muscle-driven forward dynamic active hybrid FE-MB model of the lumbosacral spine (LSS) in ArtiSynth to simultaneously calculate muscle activation patterns, vertebral movements, and internal mechanical loads. The model consisted of the rigid vertebrae L1-S1 interconnected with hyperelastic fiber-reinforced FE intervertebral discs, ligaments, facet joints, and force actuators representing the muscles. Morphological muscle data were implemented via a semi-automated registration procedure. Four auxiliary bodies were utilized to describe non-linear muscle paths by wrapping and attaching the anterior abdominal muscles. This included an abdominal plate whose kinematics was optimized using motion capture data from upper body movements. Intra-abdominal pressure was calculated from the forces of the abdominal muscles compressing the abdominal cavity. For the muscle-driven approach, forward dynamics assisted data tracking was used to predict muscle activation patterns that generate spinal postures and balance the spine without prescribing accurate spinal kinematics. During calibration, the maximum specific muscle tension and spinal rhythms resulting from the model dynamics were evaluated. To validate the model, load cases were simulated from −10° extension to +30° flexion with weights up to 20 kg in both hands. The biomechanical model responses were compared with in vivo literature data of intradiscal pressures, intra-abdominal pressures, and muscle activities. The results demonstrated high agreement with this data and highlight the advantages of active hybrid modeling for the LSS. Overall, this new self-contained tool provides a robust and efficient estimation of LSS biomechanical responses under in vivo similar loads, for example, to improve pain treatment by spinal stabilization therapies.

1 Introduction

Low back pain is a leading cause of severe activity limitations in the personal and professional lifes of sufferers (Duthey, 2013) and affects most adults during their lifetime (Andersson, 1998; Setchell et al., 2019). Many of the underlying causes and mechanisms that lead to back pain are not yet fully understood (Panjabi, 2006; Reeves et al., 2019). With the complexity of the interplay between a variety of anatomical structures, the lumbar spine is an important topic of biomechanical research (Martin et al., 2018). To increase our understanding, a wide range of high-quality in vivo and in vitro spine studies have been performed over the past decades (Sato et al., 1999; Wilke et al., 2001; Heuer et al., 2007; Oxland, 2016). Even latest experimental methods, however, are reaching their limits in studying internal mechanics like muscle forces or stress states in soft tissues, the overall context of the stabilizing functions of individual muscles (Santaguida and McGill, 1995; Freeman et al., 2010; Arbanas et al., 2013), and thus the possible causes of pain (Panjabi, 2003; Dreischarf et al., 2015; Hodges et al., 2015; Reeves et al., 2019). Also, in terms of cost reduction, ethical justifiability, and broader patient populations, numerical simulation methods are increasingly being used with varying focus, level of detail, and solution approach (Galbusera and Wilke, 2018; Knapik et al., 2022).

Most validated and extensively used biomechanical spine models belong to either the musculoskeletal multibody (MB) (de Zee et al., 2007; Christophy et al., 2012; Bruno et al., 2015; Senteler et al., 2016; Malakoutian et al., 2018; Bayoglu et al., 2019; Lerchl et al., 2022; Meszaros-Beller et al., 2023) or implicit finite element (FE) method (Schmidt et al., 2006; Ayturk and Puttlitz, 2011; Dreischarf et al., 2014; Naserkhaki et al., 2018; Turbucz et al., 2022). Musculoskeletal MB models are used to analyze mechanisms consisting of rigid bodies subject to mechanically simplifying joints and constraints. Studies thus examined interindividual muscle activation patterns, joint reaction forces, or vertebral movements in varying postures during different activities. Resulting muscle redundancy problems are commonly solved inverse dynamically (de Zee et al., 2007; Erdemir et al., 2007; Bayoglu et al., 2019). New approaches omit the specification of accurate kinematic data for the spine by forward dynamics simulation of a fully articulated spine (Rupp et al., 2015; Malakoutian et al., 2018; Meszaros-Beller et al., 2023) and including stiffening effects caused by compressive loads (Wang et al., 2020). With strength in dynamic full body approaches, MB models are limited in the study of structural behavior including the load distribution between discs, facet joints, and ligaments as well as the direct prediction of intradiscal pressures (IDP) (Azari et al., 2018). FE analyses, on the other hand, enable a detailed investigation of the structural behavior of the passive ligamentous spine and thus an estimation of internal strains and stresses by discretizing deformable components and contacts (Ayturk and Puttlitz, 2011). Effects of muscle forces and body weight, however, need to be modeled by forces and moments as simplified boundary conditions that do not represent realistic in vivo posture or loading conditions, essential to provide a realistic mechanical environment (Favier et al., 2021b; Rajaee et al., 2021; Abbasi-Ghiri et al., 2022).

In particular, clinical problems often exceed the sole capabilities of FE or MB models (Lipphaus et al., 2021), which is why model coupling (Khoddam-Khorasani et al., 2018; Liu et al., 2018; Khoddam-Khorasani et al., 2020; Favier et al., 2021b; Honegger et al., 2021; Kumaran et al., 2021; Panico et al., 2021) is increasingly used. At least two separate spine models using different methods and solvers, are coupled to use results of the other model as complementary input data. Mostly, loading modes are first calculated using an inverse dynamic musculoskeletal MB model and then applied to the passive elements of an implicit FE model for a subsequent detailed structural mechanical analysis (Nispel et al., 2023). That allows the prediction of non-linear internal strains and stresses under combined loading modes that better mimic in vivo loads. However, a coupling that leads to valid results exceeds the purely technical challenge of a manual or automated data transfer between two models built in different programs. For results from one model to be used as a valid input to another, both models must be similar, or at best the same, in their mechanical response and morphology. Even if these conditions seem obvious, it is essentially the reason why FE and MB models are coupled: Their approaches to building passive motion segments are inherently different. Liu et al. (2018) described the challenges of adjusting the biomechanical responses of both models under similar loads, arising, for example, from the fact that the intervertebral discs of the MB model were simulated by spherical joints and thus do not allow for deformations, while they were deformable in the FE model. Large movements or high loads reduced the model synchronicity, which required individual corrective measures to be taken (Liu et al., 2018). When using an upstream FE model to initially estimate the nonlinear spine stiffness of the MB model, only simplified load cases are available and the resulting FE motion is the input kinematics for the MB model (Panico et al., 2021). To circumvent this problem, Khoddam-Khorasani et al. (2018) solved a coupled model staggered-iteratively until the FE and MB solutions were convergent. Common to current coupled simulations is that the MB models used use an inverse dynamics approach. The results are thus mainly dependent on accurate kinematic data, which most of the current technologies can neither provide (Meszaros-Beller et al., 2023) nor are available for an in vivo equivalent of a MB model.

A hitherto less established way to overcome drawbacks of coupled model simulations are active hybrid FE-MB models. They represent an approach to combining the complementary strengths and limitations of both types of modeling for more predictive models (Lloyd et al., 2019; Knapik et al., 2022). Passive hybrid models provide the inherent combination and interaction of multiple elastic FE and rigid bodies in one model (Stavness et al., 2011). Earlier motivation for passive hybrid spine modeling was the increase of computational efficiency, with advantages for a simplified model structure and increased usability in clinical routine (Deacy et al., 2010; Moramarco et al., 2010; Coombs et al., 2011; Dicko et al., 2015). Likewise, this enabled dynamic solving of complex biomechanical systems with large deformations using explicit FE environments (Shirazi-Adl et al., 2002; Knapik et al., 2012; Rao, 2012). A hybrid model that comprises additional force actuators such as muscles can be referred to as an active hybrid model (Alizadeh et al., 2019). In 2012 Knapik et al. introduced an active hybrid model of the lumbosacral spine (LSS) built in a multi-body dynamic simulation environment to investigate the effects of a total disc replacement at L5/S1 level. Therefore, muscle forces were represented as force vectors driven with electromyographic (EMG) data from a healthy subject. This allowed comparison of stresses in the intervertebral discs and facet joints as well as the range of motions of the individual segments before and after disc replacement. However, their model did not focus on the computational prediction of muscle forces to drive and stabilize the spine. Another active hybrid approach was presented by Rajaee et al. (2021). In their model, the active musculature of a MS model (Arjmand and Shirazi-Adl, 2006a) was integrated into a FE model of the thoracolumbar spine (Shirazi-Adl, 1994) with the aim of fully incorporating the nonlinear stiffness properties of the passive spine. Muscle activities were determined using static optimization algorithm.

In summary, current active hybrid and coupled LSS models have in common that mostly by using established simulation environments, the potentials of the different modeling approaches can only be utilized to a limited extent. Understanding of spinal stability highlights this, as the spinal stabilization system is based on the interaction of three subsystems (Panjabi, 2003; Reeves et al., 2019): the intrinsic stability of the passive spine, the dynamic stability of the spinal muscles surrounding the spine, and the neural control unit that coordinates muscle responses. In more detail, this also includes, for example, the stabilizing influence of the intra-abdominal pressure (IAP) (Hodges et al., 2005), the exact morphology of the facet joints (Khoddam-Khorasani et al., 2018; Knapik et al., 2022), the load sharing (Liu et al., 2018), and the non-linear stiffening of the intervertebral discs under compression (Edwards et al., 1987; Wang et al., 2020). Consequently, an inherent and sufficiently detailed combination of the relevant passive and active spinal structures is necessary, without prescribing accurate kinematics.

Such models, capable of biomechanically investigating the concept of spinal stability under both physiological and pathological conditions (e.g., fusions, injuries, and degenerations), are much needed. Towards this goal, we aim to:

1) Develop a novel forward dynamic active hybrid FE-MB model of the LSS that incorporates the three subsystems involved in spinal stability and moves and balances through a muscle-driven approach.

2) Calibrate and validate the biomechanical model responses of this active hybrid model with literature data measured in vivo during various activities.

2 Material and methods

In the following, we cover our approach to building the active hybrid model of the LSS (Figure 1), followed by calibrating the maximum specific muscle tension, calibrating the lumbar segmental rotation contributions of the vertebral target frames, and the systematic testing procedure of the models’ mechanical responses under different load cases (LCs) for validation. The Java-based open-source framework ArtiSynth (www.artisynth.org), a physics simulator that supports the combined simulation of MB and FE models (Lloyd et al., 2012), was used to implement and run the active hybrid model. To perform musculoskeletal geometry modeling and visualize biomedical data in the preprocessing phase we used the software application NMSBuilder (v2.1) (Valente et al., 2017). Downstream evaluations of simulation results were carried out with Matlab (R2020b, MathWorks Inc., US). To encourage open science in spine biomechanics, the model data and procedures developed are described in detail or made freely available.

FIGURE 1

FIGURE 1. Muscle-driven forward dynamic active hybrid lumbosacral spine model with muscles in red and cyan colored auxiliary rigid bodies (Table 1): 1) Abdominal plate, 2) lumbar wrapping body, 3) left and 4) right thoracic wrapping body. Visualization from three views (A) right lateral, (B) right front, and (C) left rear.

2.1 Passive anatomical model

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.

2.2 Musculoskeletal geometry

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.

FIGURE 2

FIGURE 2. Block diagram of the semi-automated muscle registration process with the relevant in- and output data. The process outlined in the gray box is divided into five general steps (A–E). For each step, the most relevant software, NMSBuilder or Matlab, is indicated. Vertebra L1 serves as an illustration in steps (A–D). The manual adjustment of the muscle paths in e is shown for psoas major.

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.

FIGURE 3

FIGURE 3. Visualization of the muscles of the active hybrid lumbosacral spine model: Latissimus dorsi (LD), quadratus lumborum (QL), multifidus (MF), obliquus internus abdominis (IO), obliquus externus abdominis (EO), transversus abdominis (TA), psoas major (PM), rectus abdominis (RA), iliocostalis thoracis (IT), iliocostalis lumborum (IL), longissimus thoracius (LT), and longissimus lumborum (LL). (A) All muscles differentiated by color in anterior and posterior view. (B) One side of each muscle plus the respective muscle points of both muscle sides (black dots) from three views each: anterior, left lateral, and posterior.

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.

TABLE 1

TABLE 1. Geometric description of auxiliary rigid bodies for muscle interactions in reference to the global coordinate system. Dimensions and transformations are given in x, y, z coordinates following the NMSBuilder modeling convention. Visualization of the numbered auxiliary bodies in Figure 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.

TABLE 2

TABLE 2. Overview of the muscle data for the right side of the model.

2.3 Musculotendon models

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_{ref}}^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^{MT}$, the maximum isometric reference muscle force $F_{0_{ref}}^M$ was scaled according to Eq. 1 (Correa and Pandy, 2011).

$F_0^M=F_{0_{ref}}^M·\frac{b}{b_{ref}}·\frac{l_{ref}^{MT}}{l^{MT}}(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.

2.4 Abdominal kinematic

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.

$CC\left({\alpha }_{j=x,y,z}\right)=p_1{\alpha }_j+p_2{\alpha }_j^2(2)$

FIGURE 4

FIGURE 4. Implementation of the open kinematic chain to control the abdominal plate. (A) Visualization of the three joints (HJ, UJ, CJ) with their local coordinate systems connecting the abdominal plate with the sacrum (cf. Table 3). Virtual representatives of the motion capture markers are shown as red dots. (B) Lateral view of the three joints with their five controlled coordinates, visualized as arrows with black frames: Rotation angles θ₁, θ₂, θ₃, and φ as well as translation distance Z.

TABLE 3

TABLE 3. Abdominal kinematic chain joint settings with polynomial coefficients p_i. Given are the global joint centers (x, y, z) and the rotational transformations in respective axis order. Axes after an elementary rotation are denoted with apostrophes. Case differentiations were made for flexion and extension.

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.

FIGURE 5

FIGURE 5. Procedure for determining the polynomial coefficients p_i of the abdominal kinematic chain. (A) Marker setup applied to the participants skin. All marker named with a trailing index of 1-3 were fixed to rigid equilateral triangles with respective distances of 60 mm. In lateral view line 1 is crossing the spinous process of L4 and line 2 runs at the level of the upper edge of the iliac crest. (B) Sample representation of processed motion capture measurement data for one participant in flexion and extension. Color coded motion curves from neutral posture to flexion for markers LAL41, CSI1, and SPP4. (C) Summarization of three motion curves of a participant for the marker LAL41. (D) All consulted motion capture curves in sagittal plane in flexion for marker LAL41 with the calculated DTW reference curves as well as the graph of the polynomial functions. Plotted for this purpose are the trajectory components of the same marker at the abdominal plate in the simulation model with optimized p_i (Table 3).

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.

2.5 Intra-abdominal pressure

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_{IAP}=\frac{1}{c_w}·\frac{d_{AC}}{4h_{AC}}{F}_{AP}=\frac{1}{4}\pi d_{AC}^2{p}_{IAP}(3)$

FIGURE 6

FIGURE 6. Implementation of the intra-abdominal pressure (IAP). (A) Differentiation between the kinematically controlled abdominal plate (AP) and the dynamic active AP (AP_IAP). The muscle via points of RA were attached to AP. The muscle attachment points of IO, EO, and TA were located on AP_IAP. Obligatory muscle groups shown in red were IO, EO, and RA. TA was optional and could be deactivated (w/oTA) in the simulation. (B) Abdominal cavity with anatomical dimensions, acting forces, and boundary conditions. (C) Idealized cylinder with geometric relations between surfaces and forces.

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).

2.6 Simulation procedure

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).

FIGURE 7

FIGURE 7. Simulation procedure details. (A) Lumped segment parameters, loads applied to the model, and definition of dimensions in upright posture with respect to the origin and the global coordinate system. (B) To better represent start and end situations of a simulation, thorax and vertebral body target frames, colored in turquoise, are shown in flexion with ${\alpha }_y^{*}$ = 30° while the model is in the unloaded neutral state. During simulation, all target frames were successively moved so that the motion target term could be continuously minimized by the TC. For the final static condition, the model and the thorax target frame are shown in (C) +30° flexion (F30₁) and (D) −10° extension (E10₁).

2.7 Calibration of maximum specific muscle tension

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 (${\alpha }_y^{*}$ = −10°), 870 N in flexion (${\alpha }_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 $∆{\alpha }_y$ as the rotational deviation between the thorax target frame and the thorax.

$∆{\alpha }_y=|{\alpha }_y^{*}-{\alpha }_y|(4)$

Regarding the upright posture after settling (see section 2.9 Validation), ${\alpha }_y^{*}$ defines the rotation of the thorax target frame and ${\alpha }_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 $∆{\alpha }_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², $∆{\alpha }_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 $∆{\alpha }_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.

2.8 Calibration of segmental target frame rotation contributions

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

R2) Relative contributions averaged over 10°–30° ROM from Wong et al. (2006)

R3) Percentages published by Arjmand and Shirazi-Adl (2006a)

R4) Mean relative segment contributions over 30%–100% L2-S1 ROM in flexion from Breen et al. (2021) with L1-L2 equal to L2-L3

R5) Mean relative segment contributions over 30%–85% L2-S1 ROM in flexion from Breen et al. (2021) with L1-L2 equal to 80% L2-L3

R6) Mean relative segment contributions over 15%–70% L2-S1 ROM in return from Breen et al. (2021) with L1-L2 equal to 80% L2-L3

R7) Relative segmental contributions averaged over 30°–10° flexion to upright from Aiyangar et al. (2015) with L1-L2 equal to 80% L2-L3

FIGURE 8

FIGURE 8. Resulting spinal rhythms for the segmental target frame rotation contributions R1 to R7 used in calibration. All values are given in relation to the respective upright posture (α_y = 0°). (A) Comparison of the target frame rotation contributions for L2-L3 to L5-S1 with the resulting IVRs for the three simulated thoracic tilts ${\alpha }_y^{*}$ = −10, +15, and +30°. (B) IDP changes for level L4/5 compared with in vivo measurements for the same thoracic tilts.

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° < ${\alpha }_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.

2.9 Validation

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° ≤ ${\alpha }_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.

TABLE 4

TABLE 4. Experimental plan of the model validation. LCs were taken from in vivo studies. l_Load is specified for the given posture in the respective static state.

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 ${\alpha }_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).

3 Results

The resultant lumbar lordosis of the musculoskeletal biomechanical model in upright position (N0₁) measured from the superior endplate of L1 to the superior endplate of S1 was 47.5°. The IVRs between neutral vertebral orientation and the upright reference posture were from L1 to L5: 4.82, 1.88, −2.71, −3.37, and 3.12°. For all LCs, $∆{\alpha }_y$ was less than 0.22°. The estimated IDPs of the three caudal levels in different postures without additional load are shown in Figure 9A and with different F_Load in Figure 9B. The load on the lumbar spine was lowest in N0₁ and a minimum pressure was predicted at level L3/4 with 0.59 MPa. From this, the resulting IDPs for L4/5 and L5/S1 increased approximately linearly in flexion up to 20°. Above that, the pressure response was no longer linear. The IDP increased caudally per level. A maximum IDP of 2.39 MPa occurred at level L5/S1 in 30° flexion holding 8 kg with outstretched arms (F30₂). Comparing the absolute IDP values w/TA as well as their relative changes, there is a high agreement with the in vivo measurements (Schultz et al., 1982; Sato et al., 1999; Wilke et al., 2001; Takahashi et al., 2006), shown in Figure 9. However, in vivo values for upright standing without external load (Takahashi et al., 2006) and with 8 kg in outstretched arms (N0₂) (Schultz et al., 1982) are about half lower than simulated. This was not the case for the data of Wilke et al. (2001) and Sato et al. (1999) as well as for all pressure measurements in sagittal deflection. Except for 10° flexion, predicted IDP deviations were on average lower than ±4.4% for the three caudal levels in LCs w/oTA compared to LCs w/TA. Least changes with maximum −1.5% resulted for all LCs in 30° flexion. For N0₁, the IDPs w/oTA were increased in caudal direction by 3.7, 3.3, and 2.6%. Mean deviations for F10₁ and F10₂ were 19.1% and 6%, respectively.

FIGURE 9

FIGURE 9. Comparison between IDP results for intervertebral discs L3/4, L4/5, and L5/S1 w/TA with in vivo literature data. (A) LCs without additional loads in different sagittal postures. Given are the thorax target angles α_y and the predicted L4-L5 IVRs to compare with measurements by Sato et al. (1999). (B) LCs with external load (F_Load ≠ 0 N) in different postures.

The muscle forces determined by the TC to drive and stabilize the model are visualized in Figures 10A–C for the validation LCs w/TA from Table 4. For better comparability with in vivo literature data (Schultz et al., 1982; Arjmand and Shirazi-Adl, 2006a; Takahashi et al., 2006), the paraspinal muscle forces excluding MF were summarized as erector spinae (E.S.) and abdominal muscles (A.M.) excluding TA. MF and TA are visualized separately. The sum of all muscle forces in N0₁ including TA was 351 N. In 10, 20, and 30° flexion (F10₁, F20₁, and F30₁) the sums were 468, 636, and 1084 N, respectively, and 374 N in E10₁. The resulting muscle activation patterns were individual for each LC. That the highest back muscle forces occurred when the thorax was flexed 30° and 8 kg was held in both hands with outstretched arms (F30₂) is consistent with the measurements of Schultz et al. (1982). Similarly, the estimated force of the abdominal muscles decreased with increasing flexion and F_Load. Arjmand and Shirazi-Adl (2006a) also found in vivo an increased activity of E.S. and MF in flexion. Correlating with higher weight in the hands, both their activity increased. The increase in E.S. activity in flexion with and without a load of 10 kg measured in vivo (Takahashi et al., 2006) correlates closely with our calculated force changes (Figures 10A, B). Consistent with this in vivo study, additional loads did not induce significant changes in the predicted abdominal muscle activity. High agreement was also seen for MF inactivity in upright position: Holding a light to heavy weight (N0₂ - N0₅) barely led to force changes in the model as well as changes in muscle activity based on EMG measurements by Arjmand and Shirazi-Adl (2006a). To keep the upper body upright, only the muscle groups of the E.S. were activated to a greater extent. In general, in agreement with EMG measurements (Arjmand and Shirazi-Adl, 2006a), Figures 10A, B show a decrease in abdominal muscle activity in flexion. In contrast, Takahashi et al. (2006) measured in vivo an up to 2.7-fold increase of RA activity in 20° flexion without additional load. In our model, RA was only activated to a greater extent in maximum flexion (F30₁) and generated a force that was almost 3.5 times higher than in upright posture. In extension (E10₁), RA generated twice the force of the upright posture. The force variations of PM and QL between LCs were small, but in relation to N0₁ maximal for F30₂.

FIGURE 10

FIGURE 10. (A–C) Predicted muscle forces of the left side of the model w/TA. Forces of IT, IL, LT, and LL are combined into erector spinae (E.S.), and IO, EO, and RA are combined into abdominal muscles (A.M.). (A) Sagittal postures without additional loads. (B) LCs in which 10 kg were carried with vertically hanging arms. (C) Comparison of LCs with different lifting activities. (D) Muscle force changes w/oTA compared to values from (A–C) w/TA. Values with * exceed a change of 200%.

For comparison of the calculated muscle forces of all LCs w/oTA, the percentage changes to w/TA are visualized in Figure 10D. Generally, the least changes were seen for E.S. and the largest for A.M. When holding heavy loads in upright position (N0₂, N0₄, N0₅), MF and QL exerted forces reduced by up to 24%. Other muscle forces as well as the IDPs were approximately constant. For the upright position without heavy loads (N0₁, N0₃), forces of A.M. were increased by up to 50% and QL as well as MF by up to 10%. This had no effect on the IDPs but led to IVR deviations of up to 0.15°. In case of 10° flexion, A.M. were increased primarily due to EO by 400% and 500% for F10₁ and F10₂, respectively. In F10₁, the forces for MF and QL w/oTA were also increased by 28% and 49%. This resulted in a 16% higher IDP for L4/5 (Figure 9A). In flexion ≥20° and extension, the predictions w/oTA resulted in no relevant force changes for MF, E.S., and QL (−2.6% ± 2.8%). Only A.M and PM showed correlative changes in their forces. Effects on IVRs were smaller than 0.2°. Changes for the predicted IDPs were <1.8% for flexions greater than 10° and <3.8% in the upright position.

The IAPs calculated by our model are shown in comparison with in vivo measured values (Schultz et al., 1982; Nachemson et al., 1986; Mueller et al., 1998) in Figure 11. Highest IAP with 16.30 mmHg and lowest IAP with 7.15 mmHg occurred in F30₂ and F20₁. The calculated pressure values w/TA for the upright posture demonstrated a high agreement with measurements of Schultz et al. (1982) and Mueller et al. (1998). w/oTA the IAP was reduced by up to 65%. For F30₂, the calculated IAP was half as high as measured by Schultz et al. (1982).

FIGURE 11

FIGURE 11. Comparison of predicted IAP w/TA and w/oTA with literature data.

4 Discussion

This paper presents the modeling procedure, calibration, and validation of a novel muscle-driven active hybrid FE-MB model of the LSS. The modeling procedure included the extension of our previously validated passive hybrid LSS model (Remus et al., 2021) with trunk muscles and auxiliary bodies for muscle attachments and wrapping. Morphological muscle data were integrated via virtual palpation as part of a semi-automated registration procedure. Redundant muscle activations were predicted via a forward dynamics assisted data tracking algorithm to move and balance the LSS in different postures. The auxiliary bodies comprised an AP whose kinematics were optimized using motion capture data of the torso performed in this study. Assuming an incompressible abdominal cavity, the IAP was directly calculated from the posteriorly acting forces of the abdominal muscles on the AP. The aim of the two calibrations in this study was to obtain spinal stability. For this purpose, a specific muscle tension was first determined, followed by the determination of segmental rotation contributions for the vertebral target frames as an optimization criterion. To validate the forward dynamic model, 13 sagittally symmetric LCs from −10° extension to +30° flexion with and without different loads up to 20 kg held with both hands were simulated. Biomechanical model responses used for validation with in vivo literature data comprised IDP at L3/4 to L5/S1, IAP, and muscle activation patterns.

4.1 Active hybrid lumbosacral spine model

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).

4.2 Calibration

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.

4.3 Validation

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.

4.3.1 Muscle activation pattern

The most common approaches to distinguish for biomechanical muscle redundancy problems are optimization-driven and EMG-driven models (Mohammadi et al., 2015; Dreischarf et al., 2016). In EMG-driven models, muscle activities recorded via electrodes from in vivo measurements are used to solve the present redundancy problem (Knapik et al., 2012). If the solution of muscular redundancy is purely mathematical, it is an optimization-driven approach. The basic assumption is that at least one cost function is optimized by the central nervous system, while equilibrium constraints as well as upper and lower limits for muscle forces are fulfilled (Dreischarf et al., 2016). EMG-based optimizations combine both approaches to account for vital equilibrium constraints in addition to physiological EMG signals from human subjects to improve model predictions (Mohammadi et al., 2015). Solution methods for computing muscle forces in musculoskeletal models include forward and inverse optimization and optimal control (Erdemir et al., 2007; Stavness, 2010; Malakoutian et al., 2018). We used an optimization-driven approach to determine the muscle forces through which the spine is moved and stabilized. Considering equilibrium constraints via five target frames, the muscle activation patterns were calculated by the TC (Stavness, 2010; Stavness et al., 2010) integrated in ArtiSynth. Due to the redundancy of the muscle activities, the optimization problem was underdetermined (Stavness et al., 2010). The additional consideration of orientation constraints via target frames and muscle activity changes resulted in a multicriteria optimization. Despite the quadratic consideration of muscle cross-sections to select an efficient activation pattern, it is possible that only local optima for muscle activation patterns were found.

The muscle parameters of the two reference models (Christophy et al., 2012; Bayoglu et al., 2017) were based on in vitro studies and estimations. The parameters for pennation angle, passive stiffness, sarcomere length, and specific muscle tension were constant for most of the muscle groups. However, a recent simulation study (Malakoutian et al., 2022) showed the relevance of these paraspinal muscle parameters on the predicted spinal loads. To use reference material parameters, we scaled these linearly based on muscle cross-sections and body weight. Due to this, as well as the assumption of same parameters for the most muscle groups, a systematic error in our simulation results could not be excluded. We did not investigate possible correlations in more detail. Sensitivity studies with muscle group-specific parameters should therefore be pursued. Muscle activities determined were sensitive to varying weighting parameters. These were chosen iteratively in advance in such a way that all the LCs examined were solved in a stable manner and considered relative to one another, the weighting factors for the following terms were minimal with descending relevance: vertebral targets, muscle damping, thorax target, and muscle excitation. For a qualitative validation of the predicted muscle forces we used published EMG data, because in vivo measured muscle forces are generally not available (Erdemir et al., 2007). High agreement was shown in the relative changes of the muscle groups for all simulated LCs. Compared to another simulation study (Arshad et al., 2016), the sum of our local muscle forces (IL, LL, PM, MF, QL) was only slightly lower at about 272 N in the upright posture (N0₁). With 30° flexion of the upper body, the estimated local muscle force tripled comparably. In flexion, moreover, the passive muscle force component became greater because the sarcomere lengths were thus extended.

4.3.2 Intra-abdominal pressure

Compression of the abdominal cavity between the diaphragm and the pelvic floor is primarily regulated by activation of the enclosing abdominal musculature TA, EO, and IO (Goel and Weinstein, 1990). For validation, we computed all LCs without and with the predominantly transverse TA. Reason for this was the combination of the geometric and anatomical reference muscle data from two independent sources. It should be noted that the conversion dimensioning between the forces F_IAP and F_AP was performed only w/TA. The influence of TA turned out to be dependent on the LC considered and the biomechanical model response. It could be seen that IAP w/oTA was considerably underestimated compared to w/TA. Difficulty for validation was caused by limited in vivo data with a wide range of measured values. These span from 1.5 to 7.5 mmHg for upright standing without load held (Andersson et al., 1977; Schultz et al., 1982; Nachemson et al., 1986), which was probably measured most frequently. Measured values also differed by BMI (Cobb et al., 2005) and sex (Essendrop and Schibye, 2004). The percent increase in IDPs calculated in the model w/oTA was greater in upright posture than in flexion and extension. Comparable to another simulation study (Liu et al., 2019b), for w/TA and the resulting elevated IAP at greater flexion angles, E.S. activation as a global muscle group was reduced and IDP was reduced by up to 1.8%. In general, our results on the IDPs, which were in agreement with previous findings (McGill and Norman, 1987; Kumar, 1997; Arjmand and Shirazi-Adl, 2006b), suggest that the unloading function of the IAP is not to be overestimated. Much of the pressure was generated by the abdominal muscles, but they also had a compressive effect on the lower back. There is no way to contract the abdominal muscles without increasing the IAP (Cholewicki et al., 2002). In view of the effects on the other muscle groups as well as the current implementation, we considered TA to be a relevant muscle group in the simulation of sagittal-symmetric LCs. Furthermore, due to its role as a moment generator around the longitudinal axis (Cresswell et al., 1992), the relevance of TA in axial rotation and asymmetric LCs, is of particular importance. Consideration of the TA is also important for the examination of individuals with low back pain, for which this muscle group is often much more active (Knapik et al., 2022).

We attributed the underestimation of IAP in flexion and lifting activities (Figure 11) primarily to the fact that no explicit consideration of co-activation of the abdominal musculature was integrated as part of the solution of the muscle redundancy problem. Therefore, abdominal muscle activity also tended to be underestimated in other models with an optimization approach (Arjmand and Shirazi-Adl, 2006c; Liu et al., 2018). To overcome this, constant activations (Khoddam-Khorasani et al., 2020) or situation-dependent lower activity limits (El-Rich et al., 2004) for the abdominal muscles can be implemented in the future. Alternatively, the force F_IAP, which results mathematically from the compressive forces of the abdominal muscles, can be considered as an additional variable target when solving the muscle redundancy problem. This can result in the abdominal muscles being more activated by the TC to reach a targeted IAP, depending on the upper body posture and the weighting in the cost function.

4.3.3 Intradiscal pressure

In vivo IDP measurements (Schultz et al., 1982; Sato et al., 1999; Wilke et al., 2001; Takahashi et al., 2006) are useful and as mentioned before one of the few ways to gain quantitative insight into the biomechanics of the lumbosacral spine. The pressure directly reflects the load acting on the spinal level in question (Nachemson, 1981). In addition to the body’s dead weight and external masses, stabilizing muscular forces are also an integral part of the load. Therefore, the use of in vivo IDP measurements for the validation of simulation models was all the more relevant. In our model, the intervertebral discs were modeled as healthy and non-degenerated fiber-reinforced FE bodies with almost incompressible hyperelastic material behaviors. Assuming that the nucleus pulposus within the annulus fibrosus behaves hydrostatically throughout its volume (Nachemson, 1960), we calculated the IDP from the negative mean of the normal stresses of all FE nodes inside a nucleus pulposus (Remus et al., 2021). This assumption implies that in vivo only non-degenerated intervertebral discs can be reliably measured by means of a pressure transducer, because the proteoglycan-water gel must be present and mobile to a sufficient extent (McNally and Adams, 1992; Wilke et al., 2001). Degenerative alterations as well as the presence of collagen fibers in the gel that can cause anisotropic pressure distributions (McNally and Adams, 1992) are possible reasons for the partially considerably varying IDP in the literature. In general, IDP decreases when the intervertebral disc degenerates (Sato et al., 1999). The number of studies that measured IDP in vivo in different situations, however, was limited and primarily restricted to the levels L3/4 and L4/5. Because no suitable data were present for cranial lumbar spine levels and L1 and T12 were rigidly connected in the model, we considered only the three caudal spine levels in the validation. As in many biomechanical studies (Tafazzol et al., 2014), the thoracic region was represented as a single lumped rigid body in the current mode. Despite resulting variations in muscle forces, the influence on lumbosacral model responses was shown to be small in a simulation study (Ignasiak et al., 2016).

In comparison with the in vivo IDPs visualized in Figure 9, we found predominantly high agreement with the values calculated by our model. Only in upright position with 5 kg in each hand and without load were there relevant differences. In both cases, the model overestimated the IDP at L3/4 and L4/5. This did not apply to other LCs in the upright posture. As mentioned above, however, it was possible that the IDP measured by pressure transducers were reduced by degenerative changes or anisotropies. Also, not further traceable relieving postures of the examined subjects in upright standing could be plausible. However, especially the comparison with the measured values of Wilke et al. (2001), which were measured at L4/5 in a volunteer similar in body weight and height to the model, showed a very high agreement in all LCs. In general, the predicted load on the lumbar spine increased as the magnitude of segmental rotations increased, which was the case in both extension and flexion. Among others, also Andersson et al. (1977) and Nachemson (1965) observed an increase in pressure of comparable magnitude for flexions without, as well as with additional weights in the hands. The comparatively sharper increase in IDPs at 30° flexion we attributed to the fixed pelvis in the model. This is a limitation compared to other models (Liu et al., 2019a; Favier et al., 2021a; Honegger et al., 2021) that take into account a lumbopelvic ratio for larger movements and more activities. In vivo studies had shown (Tafazzol et al., 2014) that the upper body flexion is characterized by a simultaneous rotation of vertebrae and pelvis. The contribution of pelvis rotations increased with larger flexions. Small flexions were predominantly accomplished by vertebral rotations. A lumbopelvic rhythm should therefore be implemented in the future, especially for larger flexion angles.

5 Conclusion

This work is a new approach to computational spine biomechanics that combines optimization-driven trunk musculature, FE intervertebral discs, ligaments, and facet joints in one model. The resulting capability to use the interplay of passive and neurally coordinated active mechanisms to simulate spinal stability represents an advance on the state-of-the-art. Despite the discussed simplifications, the realized muscle-driven forward dynamic active hybrid model enables valid, robust, and efficient estimation of biomechanical responses of the lumbosacral spine under in vivo similar loads. Furthermore, the findings motivate a future application of the presented methods to develop patient-specific and pathological active hybrid spine models to study more complex load cases. This can be of therapeutic interest to identify conditions with great deflections but little stress on the intervertebral discs. More individualized therapies for muscle disorders through targeted strengthening or unloading are also conceivable. However, an important field of research will be the understanding and correlation to clinically relevant pain conditions.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Ethics statement

The studies involving humans were approved by the Ethics Committee of the Medical Faculty of the Ruhr-University Bochum (23-7801 04/20/23). The studies were conducted in accordance with the local legislation and institutional requirements. The participants provided their written informed consent to participate in this study.

Author contributions

The study was conceptualized and designed by RR. RR developed, calibrated, and validated the active hybrid model, interpreted and visualized the results, and wrote the original draft. SS and RR performed the experimental motion capture after RR requested the ethics approval. RR processed and evaluated the measurement data for model building. AL provided a clinical perspective and assisted RR in developing, testing, and implementing the semi-automated muscle registration procedure. MN and BB supervised the project, and MN managed the funding acquisition. All authors contributed to the article and approved the submitted version.

Funding

We acknowledge support by the DFG Open Access Publication Funds of the Ruhr-University Bochum.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbioe.2023.1223007/full#supplementary-material

Supplementary Material S1 | Tables with the locations of the COMs applied to the models’ components and polynomial curve fit coefficients s_i for the marker LAL41 DTW reference curves.

Supplementary Material S2 | Virtual palpation dictionary for identifying anatomical characteristics of the human LSS, pelvis, abdomen, and thorax.

Supplementary Material S3 | Illustrated description of palpated target anatomies including general remarks and descriptions for the workflow in NMSBuilder.

Supplementary Material S4 | Registration atlas to facilitate the workflow by providing all LM names per bone in sorted order.

Supplementary Material S5 | Muscle modeling parameters used in the model with the muscle attachments in the global coordinate system for the right side of the body.

Supplementary Material S6 | Video showing LC F30₃ for exemplary visualization of the forward dynamic simulation procedure with target frames and without musculature in right lateral view.

References

Affolter, C., Kedzierska, J., Vielma, T., Weisse, B., and Aiyangar, A. (2020). Estimating lumbar passive stiffness behaviour from subject-specific finite element models and in vivo 6DOF kinematics. J. Biomechanics 102, 109681. doi:10.1016/j.jbiomech.2020.109681

PubMed Abstract | CrossRef Full Text | Google Scholar

Abbasi-Ghiri, A., Ebrahimkhani, M., and Arjmand, N. (2022). Novel force–displacement control passive finite element models of the spine to simulate intact and pathological conditions; comparisons with traditional passive and detailed musculoskeletal models. J. Biomechanics 141, 111173. doi:10.1016/j.jbiomech.2022.111173

CrossRef Full Text | Google Scholar

Adams, M. A., Bogduk, N., Burton, K., and Dolan, P. (2013). The biomechanics of back pain. Edinburgh: Elsevier.

Google Scholar

Aiyangar, A., Zheng, L., Anderst, W., and Zhang, X. (2015). Apportionment of lumbar L2–S1 rotation across individual motion segments during a dynamic lifting task. J. Biomechanics 48, 3709–3715. doi:10.1016/j.jbiomech.2015.08.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Alizadeh, M., Knapik, G. G., Mageswaran, P., Mendel, E., Bourekas, E., and Marras, W. S. (2019). Biomechanical musculoskeletal models of the cervical spine: A systematic literature review. Clin. Biomech. (Bristol, Avon) 71, 115–124. doi:10.1016/j.clinbiomech.2019.10.027

PubMed Abstract | CrossRef Full Text | Google Scholar

Andersson, G. B. J. (1998). Epidemiology of low back pain. Acta Orthop. Scand. 69, 28–31. doi:10.1080/17453674.1998.11744790

PubMed Abstract | CrossRef Full Text | Google Scholar

Andersson, G. B. J., Ortengren, R., and Nachemson, A. L. (1977). Intradiskal pressure, intra-abdominal pressure and myoelectric back muscle activity related to posture and loading. Clin. Orthop. Relat. Res. 129, 156–164. doi:10.1097/00003086-197711000-00018

PubMed Abstract | CrossRef Full Text | Google Scholar

Arbanas, J., Pavlovic, I., Marijancic, V., Vlahovic, H., Starcevic-Klasan, G., Peharec, S., et al. (2013). MRI features of the psoas major muscle in patients with low back pain. Eur. Spine J. 22, 1965–1971. doi:10.1007/s00586-013-2749-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Arjmand, N., and Shirazi-Adl, A. (2006a). Model and in vivo studies on human trunk load partitioning and stability in isometric forward flexions. J. Biomechanics 39, 510–521. doi:10.1016/j.jbiomech.2004.11.030

PubMed Abstract | CrossRef Full Text | Google Scholar

Arjmand, N., and Shirazi-Adl, A. (2006b). Role of intra-abdominal pressure in the unloading and stabilization of the human spine during static lifting tasks. Eur. Spine J. 15, 1265–1275. doi:10.1007/s00586-005-0012-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Arjmand, N., and Shirazi-Adl, A. (2006c). Sensitivity of kinematics-based model predictions to optimization criteria in static lifting tasks. Med. Eng. Phys. 28, 504–514. doi:10.1016/j.medengphy.2005.10.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Arshad, R., Zander, T., Bashkuev, M., and Schmidt, H. (2017). Influence of spinal disc translational stiffness on the lumbar spinal loads, ligament forces and trunk muscle forces during upper body inclination. Med. Eng. Phys. 46, 54–62. doi:10.1016/j.medengphy.2017.05.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Arshad, R., Zander, T., Dreischarf, M., and Schmidt, H. (2016). Influence of lumbar spine rhythms and intra-abdominal pressure on spinal loads and trunk muscle forces during upper body inclination. Med. Eng. Phys. 38, 333–338. doi:10.1016/j.medengphy.2016.01.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Ascani, D., Mazzà, C., Lollis, A. de, Bernardoni, M., and Viceconti, M. (2015). A procedure to estimate the origins and the insertions of the knee ligaments from computed tomography images. J. Biomechanics 48, 233–237. doi:10.1016/j.jbiomech.2014.11.041

PubMed Abstract | CrossRef Full Text | Google Scholar

Ayturk, U. M., and Puttlitz, C. M. (2011). Parametric convergence sensitivity and validation of a finite element model of the human lumbar spine. Comput. Methods Biomech. Biomed. Engin 14, 695–705. doi:10.1080/10255842.2010.493517

PubMed Abstract | CrossRef Full Text | Google Scholar

Azari, F., Arjmand, N., Shirazi-Adl, A., and Rahimi-Moghaddam, T. (2018). A combined passive and active musculoskeletal model study to estimate L4-L5 load sharing. J. Biomechanics 70, 157–165. doi:10.1016/j.jbiomech.2017.04.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Bruno, A. G., Bouxsein, M. L., and Anderson, D. E. (2015). Development and validation of a musculoskeletal model of the fully articulated thoracolumbar spine and rib cage. J. Biomech. Eng. 137, 081003. doi:10.1115/1.4030408

PubMed Abstract | CrossRef Full Text | Google Scholar

Bayoglu, R., Geeraedts, L., Groenen, K. H. J., Verdonschot, N. J. J., Koopman, B. F., and Homminga, J. (2017). Twente spine model: A complete and coherent dataset for musculo-skeletal modeling of the lumbar region of the human spine. J. Biomechanics 53, 111–119. doi:10.1016/j.jbiomech.2017.01.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Bayoglu, R., Guldeniz, O., Verdonschot, N., Koopman, B., and Homminga, J. (2019). Sensitivity of muscle and intervertebral disc force computations to variations in muscle attachment sites. Comput. Methods Biomech. Biomed. Engin 1, 1135–1143. –9. doi:10.1080/10255842.2019.1644502

PubMed Abstract | CrossRef Full Text | Google Scholar

Bender, A., and Bergmann, G. (2012). Determination of typical patterns from strongly varying signals. Comput. Methods Biomech. Biomed. Engin 15, 761–769. doi:10.1080/10255842.2011.560841

PubMed Abstract | CrossRef Full Text | Google Scholar

Bogduk, N. (2012). Clinical and radiological anatomy of the lumbar spine. Edinburgh, New York: Churchill Livingstone.

Google Scholar

Bogduk, N., Macintosh, J. E., and Pearcy, M. J. (1992). A universal model of the lumbar back muscles in the upright position. Spine 17, 897–913. doi:10.1097/00007632-199208000-00007

PubMed Abstract | CrossRef Full Text | Google Scholar

Breen, A., Carvalho, D. de, Funabashi, M., Kawchuk, G., Pagé, I., Wong, A. Y. L., et al. (2021). A reference database of standardised continuous lumbar intervertebral motion analysis for conducting patient-specific comparisons. Front. Bioeng. Biotechnol. 9, 863. doi:10.3389/fbioe.2021.745837

CrossRef Full Text | Google Scholar

Cholewicki, J., Ivancic, P. C., and Radebold, A. (2002). Can increased intra-abdominal pressure in humans be decoupled from trunk muscle co-contraction during steady state isometric exertions? Eur. J. Appl. Physiol. 87, 127–133. doi:10.1007/s00421-002-0598-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Christophy, M., Faruk Senan, N. A., Lotz, J. C., and O'Reilly, O. M. (2012). A musculoskeletal model for the lumbar spine. Biomech. Model Mechanobiol. 11, 19–34. doi:10.1007/s10237-011-0290-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Cobb, W. S., Burns, J. M., Kercher, K. W., Matthews, B. D., James Norton, H., and Todd Heniford, B. (2005). Normal intraabdominal pressure in healthy adults. J. Surg. Res. 129, 231–235. doi:10.1016/j.jss.2005.06.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Coombs, D. J., Rao, M., Bushelow, M., James, D., Peter, L., and Paul, R. (2011). Simulation of lumbar spine biomechanics using abaqus. SIMULIA Customer Conference.

Google Scholar

Correa, T. A., and Pandy, M. G. (2011). A mass-length scaling law for modeling muscle strength in the lower limb. J. Biomechanics 44, 2782–2789. doi:10.1016/j.jbiomech.2011.08.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Cresswell, A. G., Grundström, H., and Thorstensson, A. (1992). Observations on intra-abdominal pressure and patterns of abdominal intra-muscular activity in man. Acta Physiol. Scand. 144, 409–418. doi:10.1111/j.1748-1716.1992.tb09314.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Daggfeldt, K., and Thorstensson, A. (2003). The mechanics of back-extensor torque production about the lumbar spine. J. Biomechanics 36, 815–825. doi:10.1016/S0021-9290(03)00015-0

PubMed Abstract | CrossRef Full Text | Google Scholar

de Zee, M., Hansen, L., Wong, C., Rasmussen, J., and Simonsen, E. B. (2007). A generic detailed rigid-body lumbar spine model. J. Biomechanics 40, 1219–1227. doi:10.1016/j.jbiomech.2006.05.030

CrossRef Full Text | Google Scholar

Deacy, J. S., Rao, M., Smith, S., Petrella, A. J., Laz, P. J., and Rullkoetter, P. J. (2010). “Combined rigid-deformable modeling of lumbar spine mechanics,” in Proceedings of the ASME Summer Bioengineering Conference 2010: Presented at 2010 ASME Summer Bioengineering Conference, Naples, Florida, USA (New York, N.Y: ASME), June 16-19, 2010, 701–702. doi:10.1115/SBC2010-19672

CrossRef Full Text | Google Scholar

Delp, S. L., Anderson, F. C., Arnold, A. S., Loan, P., Habib, A., John, C. T., et al. (2007). OpenSim: open-source software to create and analyze dynamic simulations of movement. IEEE Trans. Biomed. Eng. 54, 1940–1950. doi:10.1109/TBME.2007.901024

PubMed Abstract | CrossRef Full Text | Google Scholar

Dicko, A. H., Tong-Yette, N., Gilles, B., Faure, F., and Palombi, O. (2015). Construction and validation of a hybrid lumbar spine model for the fast evaluation of intradiscal pressure and mobility. Int. Sci. Index, Med. Health Sci. 9, 134–145. doi:10.5281/zenodo.1099356

CrossRef Full Text | Google Scholar

Dreischarf, M., Albiol, L., Zander, T., Arshad, R., Graichen, F., Bergmann, G., et al. (2015). In vivo implant forces acting on a vertebral body replacement during upper body flexion. J. Biomechanics 48, 560–565. doi:10.1016/j.jbiomech.2015.01.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Dreischarf, M., Shirazi-Adl, A., Arjmand, N., Rohlmann, A., and Schmidt, H. (2016). Estimation of loads on human lumbar spine: A review of in vivo and computational model studies. J. Biomechanics 49, 833–845. doi:10.1016/j.jbiomech.2015.12.038

PubMed Abstract | CrossRef Full Text | Google Scholar

Dreischarf, M., Zander, T., Shirazi-Adl, A., Puttlitz, C. M., Adam, C. J., Chen, C. S., et al. (2014). Comparison of eight published static finite element models of the intact lumbar spine: predictive power of models improves when combined together. J. Biomechanics 47, 1757–1766. doi:10.1016/j.jbiomech.2014.04.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Duthey, B. (2013). Background paper 6.24 Low back pain. WHO.

Google Scholar

Edwards, W. T., Hayes, W. C., Posner, I., White, A. A., and Mann, R. W. (1987). Variation of lumbar spine stiffness with load. J. Biomech. Eng. 109, 35–42. doi:10.1115/1.3138639

PubMed Abstract | CrossRef Full Text | Google Scholar

El-Rich, M., Shirazi-Adl, A., and Arjmand, N. (2004). Muscle activity, internal loads, and stability of the human spine in standing postures: combined model and in vivo studies. Spine 29, 2633–2642. doi:10.1097/01.brs.0000146463.05288.0e

PubMed Abstract | CrossRef Full Text | Google Scholar

Erdemir, A., McLean, S., Herzog, W., and van den Bogert, A. J. (2007). Model-based estimation of muscle forces exerted during movements. Clin. Biomech. (Bristol, Avon) 22, 131–154. doi:10.1016/j.clinbiomech.2006.09.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Essendrop, M., and Schibye, B. (2004). Intra-abdominal pressure and activation of abdominal muscles in highly trained participants during sudden heavy trunk loadings. Spine 29, 2445–2451. doi:10.1097/01.brs.0000143622.80004.bf

PubMed Abstract | CrossRef Full Text | Google Scholar

Essendrop, M. (2003). Significance of intra-abdominal pressure in work related trunk-loading. Dissertation. Denmark. National Institute of Occupational Health, Department of Physiology.

Google Scholar

Favier, C. D., Finnegan, M. E., Quest, R. A., Honeyfield, L., McGregor, A. H., and Phillips, A. T. M. (2021a). An open-source musculoskeletal model of the lumbar spine and lower limbs: A validation for movements of the lumbar spine. Comput. Methods Biomech. Biomed. Engin 24, 1310–1325. doi:10.1080/10255842.2021.1886284

PubMed Abstract | CrossRef Full Text | Google Scholar

Favier, C. D., McGregor, A. H., and Phillips, A. T. M. (2021b). Maintaining bone health in the lumbar spine: routine activities alone are not enough. Front. Bioeng. Biotechnol. 9, 661837. doi:10.3389/fbioe.2021.661837

PubMed Abstract | CrossRef Full Text | Google Scholar

Freeman, M. D., Woodham, M. A., and Woodham, A. W. (2010). The role of the lumbar multifidus in chronic low back pain: A review. PM R. 2, 142–146. doi:10.1016/j.pmrj.2009.11.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Galbusera, F., and Wilke, H.-J. (2018). Biomechanics of the spine: Basic concepts, spinal disorders and treatments. London: Academic Press.

Google Scholar

Goel, V. K., and Weinstein, J. N. (1990). Biomechanics of the spine: Clinical and surgical perspective. Boca Raton, Fla: CRC Press.

Google Scholar

Han, K.-S., Zander, T., Taylor, W. R., and Rohlmann, A. (2012). An enhanced and validated generic thoraco-lumbar spine model for prediction of muscle forces. Med. Eng. Phys. 34, 709–716. doi:10.1016/j.medengphy.2011.09.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Hansen, L., Zee, M. de, Rasmussen, J., Andersen, T. B., Wong, C., and Simonsen, E. B. (2006). Anatomy and biomechanics of the back muscles in the lumbar spine with reference to biomechanical modeling. Spine 31, 1888–1899. doi:10.1097/01.brs.0000229232.66090.58

PubMed Abstract | CrossRef Full Text | Google Scholar

Heuer, F., Schmidt, H., Klezl, Z., Claes, L. E., and Wilke, H.-J. (2007). Stepwise reduction of functional spinal structures increase range of motion and change lordosis angle. J. Biomechanics 40, 271–280. doi:10.1016/j.jbiomech.2006.01.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Hodges, P. W., Eriksson, A. E. M., Shirley, D., and Gandevia, S. C. (2005). Intra-abdominal pressure increases stiffness of the lumbar spine. J. Biomechanics 38, 1873–1880. doi:10.1016/j.jbiomech.2004.08.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Hodges, P. W., Cresswell, A. G., Daggfeldt, K., and Thorstensson, A. (2001). In vivo measurement of the effect of intra-abdominal pressure on the human spine. J. Biomechanics 34, 347–353. doi:10.1016/S0021-9290(00)00206-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Hodges, P. W., Ferreira, P. H., and Ferreira, M. L. (2015). “Chapter 14 - lumbar spine: treatment of motor control disorders,” in Pathology and intervention in musculoskeletal rehabilitation. Editors D. J. Magee, J. E. Zachazewski, W. S. Quillen, and R. C. Manske (s.l. Elsevier Health Sciences), 520–560.

Google Scholar

Honegger, J. D., Actis, J. A., Gates, D. H., Silverman, A. K., Munson, A. H., and Petrella, A. J. (2021). Development of a multiscale model of the human lumbar spine for investigation of tissue loads in people with and without a transtibial amputation during sit-to-stand. Biomech. Model Mechanobiol. 20, 339–358. doi:10.1007/s10237-020-01389-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Ignasiak, D., Ferguson, S. J., and Arjmand, N. (2016). A rigid thorax assumption affects model loading predictions at the upper but not lower lumbar levels. J. Biomechanics 49, 3074–3078. doi:10.1016/j.jbiomech.2016.07.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Khoddam-Khorasani, P., Arjmand, N., and Shirazi-Adl, A. (2020). Effect of changes in the lumbar posture in lifting on trunk muscle and spinal loads: A combined in vivo, musculoskeletal, and finite element model study. J. Biomechanics 104, 109728. doi:10.1016/j.jbiomech.2020.109728

PubMed Abstract | CrossRef Full Text | Google Scholar

Khoddam-Khorasani, P., Arjmand, N., and Shirazi-Adl, A. (2018). Trunk hybrid passive-active musculoskeletal modeling to determine the detailed T12-S1 response under in vivo loads. Ann. Biomed. Eng. 46, 1830–1843. doi:10.1007/s10439-018-2078-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Knapik, G. G., Mendel, E., and Marras, W. S. (2012). Use of a personalized hybrid biomechanical model to assess change in lumbar spine function with a TDR compared to an intact spine. Eur. Spine J. 21 (5), S641–S652. doi:10.1007/s00586-011-1743-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Knapik, G. G., Mendel, E., Bourekas, E., and Marras, W. S. (2022). Computational lumbar spine models: A literature review. Clin. Biomech. (Bristol, Avon) 100, 105816. doi:10.1016/j.clinbiomech.2022.105816

PubMed Abstract | CrossRef Full Text | Google Scholar

Kumar, S. (1997). The effect of sustained spinal load on intra-abdominal pressure and EMG characteristics of trunk muscles. Ergonomics 40, 1312–1334. doi:10.1080/001401397187397

PubMed Abstract | CrossRef Full Text | Google Scholar

Kumaran, Y., Shah, A., Katragadda, A., Padgaonkar, A., Zavatsky, J., McGuire, R., et al. (2021). Iatrogenic muscle damage in transforaminal lumbar interbody fusion and adjacent segment degeneration: A comparative finite element analysis of open and minimally invasive surgeries. Eur. Spine J. 30, 2622–2630. doi:10.1007/s00586-021-06909-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuo, Y.-L., Kao, C.-Y., and Tsai, Y.-J. (2021). Abdominal expansion versus abdominal drawing-in strategy on thickness and electromyography of lumbar stabilizers in people with nonspecific low back pain: A cross-sectional study. Int. J. Environ. Res. Public Health 18, 4487. doi:10.3390/ijerph18094487

PubMed Abstract | CrossRef Full Text | Google Scholar

Lipphaus, A., and Witzel, U. (2019). Biomechanical study of the development of long bones: finite element structure synthesis of the human second proximal phalanx under growth conditions. Anat. Rec. Hob. 302, 1389–1398. doi:10.1002/ar.24006

PubMed Abstract | CrossRef Full Text | Google Scholar

Le Huec, J. C., Thompson, W., Mohsinaly, Y., Barrey, C., and Faundez, A. (2019). Sagittal balance of the spine. Eur. Spine J. 28, 1889–1905. doi:10.1007/s00586-019-06083-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Lerchl, T., El Husseini, M., Bayat, A., Sekuboyina, A., Hermann, L., Nispel, K., et al. (2022). Validation of a patient-specific musculoskeletal model for lumbar load estimation generated by an automated pipeline from whole body CT. Front. Bioeng. Biotechnol. 10, 862804. doi:10.3389/fbioe.2022.862804

PubMed Abstract | CrossRef Full Text | Google Scholar

Lewandowski, A. (1982). Issues in model validation. Angew. Syst. 3.

Google Scholar

Lipphaus, A., Uttich, E., Remus, R., Hoffmann, A., Bender, B., and Witzel, U. (2021). “Applications and limitations of finite element method, multibody simulation, and hybrid modeling in bone biomechanics,” in Osteology. III. MuSkITYR symposium (Georg Thieme Verlag).

Google Scholar

Liu, T., Khalaf, K., Adeeb, S., and El-Rich, M. (2019a). Effects of lumbo-pelvic rhythm on trunk muscle forces and disc loads during forward flexion: A combined musculoskeletal and finite element simulation study. J. Biomechanics 82, 116–123. doi:10.1016/j.jbiomech.2018.10.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, T., Khalaf, K., Adeeb, S., and El-Rich, M. (2019b). Numerical investigation of intra-abdominal pressure effects on spinal loads and load-sharing in forward flexion. Front. Bioeng. Biotechnol. 7, 428. doi:10.3389/fbioe.2019.00428

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, T., Khalaf, K., Naserkhaki, S., and El-Rich, M. (2018). Load-sharing in the lumbosacral spine in neutral standing & flexed postures - a combined finite element and inverse static study. J. Biomechanics 70, 43–50. doi:10.1016/j.jbiomech.2017.10.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Lloyd, J. E., Sánchez, A., Widing, E., Stavness, I., Fels, S., Niroomandi, S., et al. (2019). “New techniques for combined FEM-multibody anatomical simulation,” in New developments on computational methods and imaging in biomechanics and biomedical engineering (Cham: Springer), 75–92.

CrossRef Full Text | Google Scholar

Lloyd, J. E., Stavness, I., and Fels, S. (2012). “ArtiSynth: A fast interactive biomechanical modeling toolkit combining multibody and finite element simulation,” in Soft tissue biomechanical modeling for computer assisted surgery. Editor Y. Payan (Springer), 355–394.

CrossRef Full Text | Google Scholar

Malakoutian, M., Sánchez, C. A., Brown, S. H. M., Street, J., Fels, S., and Oxland, T. R. (2022). Biomechanical properties of paraspinal muscles influence spinal loading-A musculoskeletal simulation study. Front. Bioeng. Biotechnol. 10, 852201. doi:10.3389/fbioe.2022.852201

PubMed Abstract | CrossRef Full Text | Google Scholar

Malakoutian, M., Street, J., Wilke, H.-J., Stavness, I., Fels, S., and Oxland, T. R. (2018). A musculoskeletal model of the lumbar spine using ArtiSynth – development and validation. Comput. Methods Biomechanics Biomed. Eng. Imaging & Vis. 6, 483–490. doi:10.1080/21681163.2016.1187087

CrossRef Full Text | Google Scholar

Malakoutian, M., Volkheimer, D., Street, J., Dvorak, M. F., Wilke, H.-J., and Oxland, T. R. (2015). Do in vivo kinematic studies provide insight into adjacent segment degeneration? A qualitative systematic literature review. Eur. Spine J. 24, 1865–1881. doi:10.1007/s00586-015-3992-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Martin, J. T., Gullbrand, S. E., Fields, A. J., Purmessur, D., Diwan, A. D., Oxland, T. R., et al. (2018). Publication trends in spine research from 2007 to 2016: comparison of the orthopaedic research society spine section and the international society for the study of the lumbar spine. JOR Spine 1, e1006. doi:10.1002/jsp2.1006

PubMed Abstract | CrossRef Full Text | Google Scholar

McGill, S. M., and Norman, R. W. (1987). Reassessment of the role of intra-abdominal pressure in spinal compression. Ergonomics 30, 1565–1588. doi:10.1080/00140138708966048

PubMed Abstract | CrossRef Full Text | Google Scholar

McNally, D. S., and Adams, M. A. (1992). Internal intervertebral disc mechanics as revealed by stress profilometry. Spine 17, 66–73. doi:10.1097/00007632-199201000-00011

PubMed Abstract | CrossRef Full Text | Google Scholar

Meszaros-Beller, L., Hammer, M., Riede, J. M., Pivonka, P., Little, J. P., and Schmitt, S. (2023). Effects of geometric individualisation of a human spine model on load sharing: neuro-musculoskeletal simulation reveals significant differences in ligament and muscle contribution. Biomech. Model Mechanobiol. 22, 669–694. doi:10.1007/s10237-022-01673-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Millard, M., Uchida, T., Seth, A., and Delp, S. L. (2013). Flexing computational muscle: modeling and simulation of musculotendon dynamics. J. Biomech. Eng. 135, 021005. doi:10.1115/1.4023390

PubMed Abstract | CrossRef Full Text | Google Scholar

Modenese, L., Montefiori, E., Wang, A., Wesarg, S., Viceconti, M., and Mazzà, C. (2018). Investigation of the dependence of joint contact forces on musculotendon parameters using a codified workflow for image-based modelling. J. Biomechanics 73, 108–118. doi:10.1016/j.jbiomech.2018.03.039

PubMed Abstract | CrossRef Full Text | Google Scholar

Mohammadi, Y., Arjmand, N., and Shirazi-Adl, A. (2015). Comparison of trunk muscle forces, spinal loads and stability estimated by one stability- and three EMG-assisted optimization approaches. Med. Eng. Phys. 37, 792–800. doi:10.1016/j.medengphy.2015.05.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Moramarco, V., Pérez del Palomar, A., Pappalettere, C., and Doblaré, M. (2010). An accurate validation of a computational model of a human lumbosacral segment. J. Biomechanics 43, 334–342. doi:10.1016/j.jbiomech.2009.07.042

PubMed Abstract | CrossRef Full Text | Google Scholar

Mueller, G., Morlock, M. M., Vollmer, M., Honl, M., Hille, E., and Schneider, E. (1998). Intramuscular pressure in the erector spinae and intra-abdominal pressure related to posture and load. Spine 23, 2580–2590. doi:10.1097/00007632-199812010-00013

PubMed Abstract | CrossRef Full Text | Google Scholar

Nachemson, A. L. (1960). Lumbar intradiscal pressure: experimental studies on post-mortem material. Acta Orthop. Scand. 31, 1–104. doi:10.3109/ort.1960.31.suppl-43.01

PubMed Abstract | CrossRef Full Text | Google Scholar

Nachemson, A. L., Andersson, G. B. J., and Schultz, A. B. (1986). Valsalva maneuver biomechanics. Effects on lumbar trunk loads of elevated intraabdominal pressures. Spine 11, 476–479. doi:10.1097/00007632-198606000-00015

PubMed Abstract | CrossRef Full Text | Google Scholar

Nachemson, A. L. (1981). Disc pressure measurements. Spine 6, 93–97. doi:10.1097/00007632-198101000-00020

PubMed Abstract | CrossRef Full Text | Google Scholar

Nachemson, A. L. (1965). The effect of forward leaning on lumbar intradiscal pressure. Acta Orthop. Scand. 35, 314–328. doi:10.3109/17453676508989362

PubMed Abstract | CrossRef Full Text | Google Scholar

Naserkhaki, S., Arjmand, N., Shirazi-Adl, A., Farahmand, F., and El-Rich, M. (2018). Effects of eight different ligament property datasets on biomechanics of a lumbar L4-L5 finite element model. J. Biomechanics 70, 33–42. doi:10.1016/j.jbiomech.2017.05.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Nispel, K., Lerchl, T., Senner, V., and Kirschke, J. S. (2023). Recent advances in coupled MBS and FEM models of the spine—a review. Bioengineering 10, 315. doi:10.3390/bioengineering10030315

PubMed Abstract | CrossRef Full Text | Google Scholar

Oxland, T. R. (2016). Fundamental biomechanics of the spine-What we have learned in the past 25 years and future directions. J. Biomechanics 49, 817–832. doi:10.1016/j.jbiomech.2015.10.035

PubMed Abstract | CrossRef Full Text | Google Scholar

Pandy, M. G. (2001). Computer modeling and simulation of human movement. Annu. Rev. Biomed. Eng. 3, 245–273. doi:10.1146/annurev.bioeng.3.1.245

PubMed Abstract | CrossRef Full Text | Google Scholar

Panico, M., Bassani, T., Villa, T. M. T., and Galbusera, F. (2021). The simulation of muscles forces increases the stresses in lumbar fixation implants with respect to pure moment loading. Front. Bioeng. Biotechnol. 9, 745703. doi:10.3389/fbioe.2021.745703

PubMed Abstract | CrossRef Full Text | Google Scholar

Panjabi, M. M. (2006). A hypothesis of chronic back pain: ligament subfailure injuries lead to muscle control dysfunction. Eur. Spine J. 15, 668–676. doi:10.1007/s00586-005-0925-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Panjabi, M. M. (2003). Clinical spinal instability and low back pain. J. Electromyogr. Kinesiol 13, 371–379. doi:10.1016/s1050-6411(03)00044-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Patwardhan, A. G., Havey, R. M., Meade, K. P., Lee, B., and Dunlap, B. (1999). A follower load increases the load-carrying capacity of the lumbar spine in compression. Spine 24, 1003–1009. doi:10.1097/00007632-199905150-00014

PubMed Abstract | CrossRef Full Text | Google Scholar

Pickering, E., Pivonka, P., and Little, J. P. (2021). Toward patient specific models of pediatric IVDs: A parametric study of ivd mechanical properties. Front. Bioeng. Biotechnol. 9, 632408. doi:10.3389/fbioe.2021.632408

PubMed Abstract | CrossRef Full Text | Google Scholar

Rajaee, M. A., Arjmand, N., and Shirazi-Adl, A. (2021). A novel coupled musculoskeletal finite element model of the spine - critical evaluation of trunk models in some tasks. J. Biomechanics 119, 110331. doi:10.1016/j.jbiomech.2021.110331

PubMed Abstract | CrossRef Full Text | Google Scholar

Rao, M. (2012). Explicit finite element Modeling of the human lumbar spine. Dissertation. Denver: University of Denver, Faculty of Engineering and Computer Science.

Google Scholar

Reeves, N. P., Cholewicki, J., van Dieën, J. H., Kawchuk, G., and Hodges, P. W. (2019). Are stability and instability relevant concepts for back pain? J. Orthop. Sports Phys. Ther. 49, 415–424. doi:10.2519/jospt.2019.8144

PubMed Abstract | CrossRef Full Text | Google Scholar

Remus, R., Lipphaus, A., Neumann, M., and Bender, B. (2021). Calibration and validation of a novel hybrid model of the lumbosacral spine in ArtiSynth-The passive structures. PLOS ONE 16, e0250456. doi:10.1371/journal.pone.0250456

PubMed Abstract | CrossRef Full Text | Google Scholar

Roussouly, P., Gollogly, S., Berthonnaud, E., and Dimnet, J. (2005). Classification of the normal variation in the sagittal alignment of the human lumbar spine and pelvis in the standing position. Spine 30, 346–353. doi:10.1097/01.brs.0000152379.54463.65

PubMed Abstract | CrossRef Full Text | Google Scholar

Rupp, T. K., Ehlers, W., Karajan, N., Günther, M., and Schmitt, S. (2015). A forward dynamics simulation of human lumbar spine flexion predicting the load sharing of intervertebral discs, ligaments, and muscles. Biomech. Model Mechanobiol. 14, 1081–1105. doi:10.1007/s10237-015-0656-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Santaguida, P. L., and McGill, S. M. (1995). The psoas major muscle: A three-dimensional geometric study. J. Biomechanics 28, 339–345. doi:10.1016/0021-9290(94)00064-B

PubMed Abstract | CrossRef Full Text | Google Scholar

Sato, K., Kikuchi, S., and Yonezawa, T. (1999). In vivo intradiscal pressure measurement in healthy individuals and in patients with ongoing back problems. Spine 24, 2468–2474. doi:10.1097/00007632-199912010-00008

PubMed Abstract | CrossRef Full Text | Google Scholar

Schmidt, H., Heuer, F., Simon, U., Kettler, A., Rohlmann, A., Claes, L. E., et al. (2006). Application of a new calibration method for a three-dimensional finite element model of a human lumbar annulus fibrosus. Clin. Biomech. (Bristol, Avon) 21, 337–344. doi:10.1016/j.clinbiomech.2005.12.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Schultz, A. B., Andersson, G. B. J., Ortengren, R., Haderspeck, K., and Nachemson, A. L. (1982). Loads on the lumbar spine. Validation of a biomechanical analysis by measurements of intradiscal pressures and myoelectric signals. J. Bone Jt. Surg. Am. 64, 713–720. doi:10.2106/00004623-198264050-00008

PubMed Abstract | CrossRef Full Text | Google Scholar

Schünke, M., Schulte, E., and Schumacher, U. (2018). PROMETHEUS Allgemeine Anatomie und Bewegungssystem: LernAtlas der Anatomie. Stuttgart: Georg Thieme Verlag.

Google Scholar

Senteler, M., Weisse, B., Rothenfluh, D. A., and Snedeker, J. G. (2016). Intervertebral reaction force prediction using an enhanced assembly of OpenSim models. Comput. Methods Biomech. Biomed. Engin 19, 538–548. doi:10.1080/10255842.2015.1043906

PubMed Abstract | CrossRef Full Text | Google Scholar

Setchell, J., Costa, N., Ferreira, M., and Hodges, P. W. (2019). What decreases low back pain? A qualitative study of patient perspectives. Scand. J. Pain 19, 597–603. doi:10.1515/sjpain-2019-0018

PubMed Abstract | CrossRef Full Text | Google Scholar

Shirazi-Adl, A. (1994). Biomechanics of the lumbar spine in sagittal/lateral moments. Spine 19, 2407–2414. doi:10.1097/00007632-199411000-00007

PubMed Abstract | CrossRef Full Text | Google Scholar

Shirazi-Adl, A., Sadouk, S., Parnianpour, M., Pop, D., and El-Rich, M. (2002). Muscle force evaluation and the role of posture in human lumbar spine under compression. Eur. Spine J. 11, 519–526. doi:10.1007/s00586-002-0397-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Smit, T. H., Odgaard, A., and Schneider, E. (1997). Structure and function of vertebral trabecular bone. Spine 22, 2823–2833. doi:10.1097/00007632-199712150-00005

PubMed Abstract | CrossRef Full Text | Google Scholar

Spitzer, V., Ackerman, M. J., Scherzinger, A. L., and Whitlock, D. (1996). The visible human male: A technical report. J. Am. Med. Inf. Assoc. 3, 118–130. doi:10.1136/jamia.1996.96236280

PubMed Abstract | CrossRef Full Text | Google Scholar

Stavness, I. (2010). Byte your tongue: A computational model of human mandibular-lingual biomechanics for biomedical applications. Dissertation. Vancouver: The University of British Columbia, Electrical and Computer Engineering.

Google Scholar

Stavness, I., Hannam, A. G., Lloyd, J. E., and Fels, S. (2010). Predicting muscle patterns for hemimandibulectomy models. Comput. Methods Biomech. Biomed. Engin 13, 483–491. doi:10.1080/10255841003762034

PubMed Abstract | CrossRef Full Text | Google Scholar

Stavness, I., Lloyd, J. E., and Fels, S. (2012). Automatic prediction of tongue muscle activations using a finite element model. J. Biomechanics 45, 2841–2848. doi:10.1016/j.jbiomech.2012.08.031

PubMed Abstract | CrossRef Full Text | Google Scholar

Stavness, I., Lloyd, J. E., Payan, Y., and Fels, S. (2011). Coupled hard-soft tissue simulation with contact and constraints applied to jaw-tongue-hyoid dynamics. Int. J. Numer. Meth. Biomed. Engng. 27, 367–390. doi:10.1002/cnm.1423

CrossRef Full Text | Google Scholar

Stokes, I. A., Gardner-Morse, M. G., and Henry, S. M. (2011). Abdominal muscle activation increases lumbar spinal stability: analysis of contributions of different muscle groups. Clin. Biomech. (Bristol, Avon) 26, 797–803. doi:10.1016/j.clinbiomech.2011.04.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Stokes, I. A., and Gardner-Morse, M. G. (1995). Lumbar spine maximum efforts and muscle recruitment patterns predicted by a model with multijoint muscles and joints with stiffness. J. Biomechanics 28, 173–186. doi:10.1016/0021-9290(94)E0040-A

PubMed Abstract | CrossRef Full Text | Google Scholar

Stokes, I. A., and Gardner-Morse, M. G. (1999). Quantitative anatomy of the lumbar musculature. J. Biomechanics 32, 311–316. doi:10.1016/S0021-9290(98)00164-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Subramani, A. V., Whitley, P. E., Garimella, H. T., and Kraft, R. H. (2020). Fatigue damage prediction in the annulus of cervical spine intervertebral discs using finite element analysis. Comput. Methods Biomech. Biomed. Engin 23, 773–784. doi:10.1080/10255842.2020.1764545

PubMed Abstract | CrossRef Full Text | Google Scholar

Tafazzol, A., Arjmand, N., Shirazi-Adl, A., and Parnianpour, M. (2014). Lumbopelvic rhythm during forward and backward sagittal trunk rotations: combined in vivo measurement with inertial tracking device and biomechanical modeling. Clin. Biomech. (Bristol, Avon) 29, 7–13. doi:10.1016/j.clinbiomech.2013.10.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Takahashi, I., Kikuchi, S., Sato, K., and Sato, N. (2006). Mechanical load of the lumbar spine during forward bending motion of the trunk-a biomechanical study. Spine 31, 18–23. doi:10.1097/01.brs.0000192636.69129.fb

PubMed Abstract | CrossRef Full Text | Google Scholar

Tayebi, S., Gutierrez, A., Mohout, I., Smets, E., Wise, R., Stiens, J., et al. (2021). A concise overview of non-invasive intra-abdominal pressure measurement techniques: from bench to bedside. J. Clin. Monit. Comput. 35, 51–70. doi:10.1007/s10877-020-00561-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Todros, S., Cesare, N. de, Concheri, G., Natali, A. N., and Pavan, P. G. (2020). Numerical modelling of abdominal wall mechanics: the role of muscular contraction and intra-abdominal pressure. J. Mech. Behav. Biomed. Mater 103, 103578. doi:10.1016/j.jmbbm.2019.103578

PubMed Abstract | CrossRef Full Text | Google Scholar

Turbucz, M., Pokorni, A. J., Szőke, G., Hoffer, Z., Kiss, R. M., Lazary, A., et al. (2022). Development and validation of two intact lumbar spine finite element models for in silico investigations: comparison of the bone modelling approaches. Appl. Sci. 12, 10256. doi:10.3390/app122010256

CrossRef Full Text | Google Scholar

Valente, G., Crimi, G., Vanella, N., Schileo, E., and Taddei, F. (2017). nmsBuilder: freeware to create subject-specific musculoskeletal models for OpenSim. Comput. Methods Programs Biomed. 152, 85–92. doi:10.1016/j.cmpb.2017.09.012

PubMed Abstract | CrossRef Full Text | Google Scholar

van Ramshorst, G. H., Salih, M., Hop, W. C. J., van Waes, O. J. F., Kleinrensink, G.-J., Goossens, R. H. M., et al. (2011). Noninvasive assessment of intra-abdominal pressure by measurement of abdominal wall tension. J. Surg. Res. 171, 240–244. doi:10.1016/j.jss.2010.02.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Sint Jan, S. (2007). Color atlas of skeletal landmark definitions: Guidelines for reproducible manual and virtual palpations. Edinburgh, New York: Churchill Livingstone.

Google Scholar

Vette, A. H., Yoshida, T., Thrasher, T. A., Masani, K., and Popovic, M. R. (2011). A complete, non-lumped, and verifiable set of upper body segment parameters for three-dimensional dynamic modeling. Med. Eng. Phys. 33, 70–79. doi:10.1016/j.medengphy.2010.09.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, K., and Gasser, T. (1997). Alignment of curves by dynamic time warping. Ann. Stat. 25, 1251–1276. doi:10.1214/aos/1069362747

CrossRef Full Text | Google Scholar

Wang, W., Wang, D., Groote, F. de, Scheys, L., and Jonkers, I. (2020). Implementation of physiological functional spinal units in a rigid-body model of the thoracolumbar spine. J. Biomechanics 98, 109437. doi:10.1016/j.jbiomech.2019.109437

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilke, H.-J., Neef, P., Hinz, B., Seidel, H., and Claes, L. E. (2001). Intradiscal pressure together with anthropometric data – A data set for the validation of models. Clin. Biomech. 16, 111–126. doi:10.1016/S0268-0033(00)00103-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Wong, K. W. N., Luk, K. D. K., Leong, J. C. Y., Wong, S. F., and Wong, K. K. (2006). Continuous dynamic spinal motion analysis. Spine 31, 414–419. doi:10.1097/01.brs.0000199955.87517.82

PubMed Abstract | CrossRef Full Text | Google Scholar

Zajac, F. E. (1989). Muscle and tendon: properties, models, scaling, and application to biomechanics and motor control. Crit. Rev. Biomed. Eng. 17, 359–411.

PubMed Abstract | Google Scholar