*3.1. MBS* → *FEM*

A unidirectional data transfer from MBS to FEM was the most common approach in unidirectional coupling. Muscle forces, reaction forces or displacements were calculated by an MBS simulation and integrated as boundary conditions (BC) in a subsequent FEM simulation.

Esat et al. [37,38] established a unidirectional MBS and FEM co-simulation of the cervical spine to calculate intradiscal pressure and AF stresses at large impact accelerations. Geometries for rigid MBS bodies were inspired by human anatomy but not specifically derived from medical images. Bodies were created in computer-aided design (CAD) software and imported in Nastran (MSC Software, Garching, Germany) before ligaments and facet joint contacts were integrated. IVD joints were modeled with nonlinear viscoelastic bushing components. One hundred and forty muscles with their respective time-dependent forces were transcribed using a dynamic calculation framework (Virtual Muscle 3.1.5 [48]). The FEM model was then implemented with independent dimensions and locations of IVDs, which were derived from the general quantitative anatomy of the spine in the software Marc/Mentat (MSC Software, Garching, Germany). NP and AF were defined using linear material parameters. Two impact loadings, namely, 15 g frontal and 8.5 g rear-end impacts, were implemented at the neck of the MBS model. The resulting loads—two sagittal forces and one sagittal moment—at the interface points for each IVD, were measured in the MBS model and implemented as time-dependent force BCs in the subsequently set-up FEM model to predict von Mises stresses in the AF and the intradiscal pressure in the NP.

Du et al. [49] investigated the impact of ejection out of, e.g., a plane on the IVDs using undirectional co-simulation. The MBS model included the entire human body seated in a chair and consisted of 16 rigid elements, with the spine divided only into the neck, upper body and lower body. Joint properties were taken from the literature. The joint between the upper and lower torso was replaced by a spring with a stiffness level of 900 N/mm. Loading was applied through an accelerative load peak of 15 g. The nonlinear FEM model of the thoracolumbar–pelvis complex (T9-S1) included the following material properties: vertebrae and pelvis as isotropic homogeneous elastic material, IVD divided into NP and AF, with the NP containing brick elements of an incompressible, hyperelastic Mooney–Rivlin material model and the AF containing hyperelastic ground substance and 3D cable elements as fibers with nonlinear stress–strain behavior, along with seven ligaments as tension-only cable elements. Attachment points and cross-sectional areas of the ligaments were obtained from Agur et al. [50]. The model was validated statically and quasi-dynamically. The MBS model was used to calculate the translation and rotation of freely chosen reference points (RP), namely, the hip joint and a point at the center of the superior endplate at T9, which were subsequently implemented as time-dependent BCs in the FEM model. For both models, Hypermesh (Altair Engineering Corp., Troy, MI, USA) was used.

Henao et al. [51] simulated surgery procedures to predict potential spinal cord damage. To do so, they implemented a complete, patient-specific FEM spine model that included the spinal cord. A previously developed MBS model of the spine [52] contained 6-DoF springs as IVD and ligament representatives. It was used to calculate the displacements of respective vertebrae, which were then incorporated as BC in the FEM model. No information was given on the stiffness parameters of the spine and how the authors ensured consistency with the FEM IVD model. The study was rated useful for surgery planning, as spinal cord injuries could be predicted by the model. Individualized parameters, especially in the area of the IVD, would likely improve prediction quality.

Honegger et al. [53] investigated the IVD stresses during the sit-to-stand (STS) transfer of lower-limb amputees using a unidirectional co-simulation. The lower-body MBS model was built in OpenSim and included 294 Hill-type muscles and five lumbar vertebrae, which were complemented with three DoF bushing elements containing stiffness values found in the literature [20]. Simulating the STS motion in the MBS model resulted in the lumbar pelvic rhythm, namely, the respective time-dependent joint angles of the lumbar spine and the muscle forces. In a lumbar FEM model taken from Campbell et al. [54], these values were introduced as time-dependent inputs, together with joint contact forces and moments. Results included, among others, the time-dependent AF stresses, facet joint loads and

IDPs during the STS. Comparing the IDPs calculated based on the MBS joint force and the IDPs found in the FEM model showed greater agreement between MBS values and in vivo data [55]. The authors stated that personalized data would possibly improve the muscle force and joint-load estimations from the MBS model, and thus have an impact on the FEM model and final simulation results.

A research group bulit around Shirazi-Adl and Arjmand has been developing numerical models of the spine since 1985 [56] and has started to advance, in parallel, two distinct model versions, which have been known by various names throughout the years. For overview reasons, the two models are herein marked as "I" and "II".


Version II was not solved with a classic MBS solver, but included clear characteristics of an MBS model, such as the involvement of mainly rigid bodies and their interconnections by joint-like components. While version II could not represent detailed deformations, it was able to perform large numbers of nonlinear analyses [57]. The researchers refined these models into a passive, osteoligamentous (I) and an active, muscular part (II) and considered their respective roles in a coupled, biomechanical system [57–59]. They implemented their passive FEM model in Abaqus (Dassault Systèmes, Paris, France), consisting of single rigid vertebrae (L1–L5) and one rigid body combining T1–T12. The IVDs were created by extruding the endplate geometries and were infused with rebar elements as AF fibers. The active model was based on a Fortran code. For the coupling process, virtual springs were attached to each vertebra in the passive model (I) to allow the transmission of shear forces and axial, lateral and sagittal moments. Hence, the load sharing between the passive structures and the muscles could be controlled.

In 2002, Shirazi-Adl et al. [60] incorporated a novel, kinematic-based Matlab (the Mathorks Inc., Natick, Massachusetts, USA) algorithm into the simplified model (II) to interactively calculate optimized solutions for muscle- and passive reaction forces. Refer to [61] for a detailed data flow chart. Note, optimization algorithms are commonly used in MBS models for muscle force estimations as well. Beam joints in the active model (II) represented the overall nonlinear stiffness levels of the spinal segments—namely, IVDs, facets and ligaments—with a nonlinear load–displacement curve, which was defined in an iterative process towards yielding the best agreement with the detailed FEM results (I) [57] and was direction-dependent [60]. Ten muscles, five global and five local, were included with six distinct fascicles [62]. The detailed FEM model (I) was based on CT images and consisted of 6 vertebrae (L1-S1) as rigid bodies, five IVDs divided ino AF and NP, 10 contact facet surfaces and 9 sets of ligaments [40]. Fourteen AF lamellae were considered with rebar fibers and a linear elastic ground substance, along with an incompressible NP. Satisfactory agreement considering predicted rotations under flexion, extension and torque loading validated the two respective models against each other [60]. Given IVD cross sectional areas, values in the FEM model (I) were 17% smaller than the ones in the MS model (II) [36].

Azari et al. [36] modified the models by applying the idea of a follower load (FL), which is a method for including simplified muscle forces in passive FEM models (I) to achieve more realistic load conditions [40]. The FL's line of action was aligned with the lumbar spine's curvature, passing the endplates and vertebrae at the approximate center. The FEM model was based on the detailed, passive FEM model (I), but included only the L4–L5 segment. The MS model (II) was adopted mainly unchanged from earlier work by the research group [56,57,60–62]. To estimate stresses and strains in IVDs; ligaments and facets; and load sharing among these structures under realistic loading conditions, Azari et al. applied gravity loads and muscle forces from an MS trunk model (II) to the passive FEM model (I) in several static positions. Twelve static positions were simulated based on the availability of the results of Wilke et al. [55] and the kinematic data gathered by

Arjmand et al. [62]. A resultant force including all gravity loads and muscle forces with upper insertions at and above L4 was calculated using the MS muscle force predictions. After the passive FEM model was rotated manually per IVD to fit the kinematic position of the twelve respective MS models, the research group tried two approaches: directly applying muscle forces estimated by the MS model in the center of the L4 body, and determining a substituting FL by trial and error to yield a similar IDP as derived by the previous approach. The FL was applied as a unilateral, pre-compressed spring between the L4 and L5 centroid.

To analyze scoliotic spine loads and their growth patterns, Kamal et al. [39] established a combined MS and FEM simulation based on subject-specific upright positional data. An MS model data of the bony spine components, including the pelvis, hip and ribcage, was derived from subject-specific CT images using mimics and was subsequently aligned into an estimated upright position with the aid of optical imaging tools. IVD centroids were considered the centers of rotation (CoRs) for the modeled spine section T12–L1 and set up with three nonlinear rotational stiffness values according to Hajihosseinali et al. [17]. Simulations were run statically with all rigid bodies fully constrained. Muscle attachment points for more than 160 muscles were calculated using Matlab. A static optimization algorithm then computed muscle forces and reaction moments, namely, individual forces F*<sup>m</sup>* and reaction moment M*CoR*, for each muscle. Reaction forces were considered as the sum of gravitational forces and muscle forces. The resulting forces and moments in a static equilibrium condition calculated in the MS model were applied to the vertebral growth plate FEM bodies as distributed pressures and shear stresses. Note, IVDs were solely considered as joints, and the biomechanical focus lay on the growth plates. FEM simulations were carried out with Abaqus (Dassault Systèmes, Paris, France). The resulting stresses in the growth plates were comparable with results from the literature and may be helpful in predicting growth patterns in scoliotic spines.

Further studies included a hybrid model applying muscle forces calculated in the beam-joint-based MS model (II) to the updated, passive T12-S1 FEM model (I) [40]. The MS model consisted of 56 muscles and was solved with an optimization algorithm minimizing the sum of the cubed muscle stresses. The outcome served as a basis for evaluating the loads on the beam joints. Simulations included eight static tasks for which the muscle forces were estimated by connector elements between muscle insertion points in the MS model. Together with the pelvic rotations gathered from in vivo experiments, these forces were substituted into the FEM model. Final deformed positions for the eight tasks were found to vary between the MS and FEM model, and the authors subsequently applied an iterative process, in which the compression-dependent stiffness properties of the MS beam-joints were adapted and resulting forces were again prescribed into the FEM model until the solution was convergent. The hybrid nature of the approach was thus the adaption of the beam-joint parameters based on a passive FEM model deformed position, which was defined actively by MS results. Kinematics of the MS and FEM models differed by less than 1 mm in the final state.

In order to develop a self-defined gold standard of spine modeling, Rajaee et al. [63] further developed the hybrid MS and FEM model introduced by Khoddam-Khorasani et al. [40] towards a coupled model, which enhances the hybrid model by incorporating the musculature of the MS model (II) in the FEM model (I). The manual approach of iteratively updating the nonlinear beam-joint stiffness was thus avoided. The resulting model hence consists of the rigid vertebrae defined earlier in the passive FEM model (I), the detailed geometries of the IVDs and all muscles and ligaments defined in the MS model. In the resulting FEM model, muscle forces are predicted by a procedure similar to the optimization procedure described above [60], and due to the novel method, depend on detailed IVD FEM model deformations. Static positions were simulated to evaluate the performance of the model compared to earlier versions and in vivo data. Stress distributions in IVDs were not available in the study, but compression forces and muscles forces could be compared between different model versions.

To sum up, multiple simulation approaches have been executed using MBS simulation results in FEM models. Transferred data included reaction forces, muscle forces or displacements of vertebrae on different spine levels. The vast majority of MBS studies rely on experimental data to define joint properties of the IVDs. However, one study investigated the potential of FEM simulations to define MBS model properties.
