**4. Bidirectional Co-Simulation of the Spine**

Bidirectional co-simulation models benefit from an iterative data exchange that requires less manual intervention and more accurately accounts for large deformations of, e.g., IVDs. Approaches using this type of coupling are listed below, again focusing on how data are exchanged between MBS and FEM models.

Monteiro et al. [46] realized a bidirectional data flow to explore intersomatic fusion biomechanics. In their study, the MBS framework included specific reference points, for which the kinematic data, namely, displacements and rotations, were calculated. The results were transferred as initial data to the FEM simulation, which calculated the reaction

forces and moments. Those were transmitted back to the MBS software and served as new starting points for another forward dynamic iteration. The process was managed by a co-simulation module, which used a gluing algorithm (an algorithm for bidirectional coupling of numerical models) to communicate between the MBS and FEM model of a C5-C6 or C6-C7 segment. This gluing algorithm was an adaptation of an algorithm developed in 2001 by Tseng et. al. [70], which in turn was based on the coordinate split (CS) technique by Yen et al. [71]. Wang et al. [72] complemented the gluing algorithm with an interfacing communication platform to make it suitable for practical applications. This approach made the single submodels black boxes, which could be coupled with one of three distinct algorithms distinguished by the data provided by the MBS system. Monteiro et al. applied the algorithm in which the MBS system (coordinator) provided kinematic data, as it was more convenient in the case of a forward dynamic analysis with given displacements. The general environment was developed in Abaqus and Apollo, a multibody system dynamics (MSD) simulator based on the Adams–Moulton method. The MBS algorithm was implemented in Fortran code. The co-simulation module core was incorporated into the Apollo code, and the co-simulation partner was Abaqus. The geometry of the IVD was derived from images and divided into NP and AF based on ratios from the literature. Material parameters were chosen as viscoelastic and quasiincompressible with a hydrostatic NP and a fiber-drawn AF [37]. IVDs non-adjacent to the fused vertebrae were modeled as bushing components for efficacy reasons. Seven ligaments were introduced as viscoelastic, nonlinear elements. For the linkage of the MBS and FEM system, Monteiro et al. implemented so-called co-simulation elements as RPs. The possibility of more than one RP per linkage was mentioned as an option to consider more complex deformations. In his work, however, one master RP and one slave RP were introduced, either to the center of the top side of the IVD (slave) or to the center of the bottom side of the IVD (master). Refer to Figure 2a for a graphical representation of these RPs. The master RP belonged to the constrained master body, while the slave RP drove the deformation of the model. The information flux therefore consisted of the three basic stations: First, the kinematic data of the RPs were analyzed and stored. Second, the FEA was launched by taking into account the kinematic data of the RPs. Third, the kinetic results of step two were processed by the coordinator software. The validation of the model with a sagittal moment of 1.5 Nm applied to the head showed realistic results, confirming the compatibility of the MBS software with the FEM analysis considering the modeling of the spine.

Another coupling method was implemented by Dicko et al. [73] in combination with a composite lumbar spine model. The algorithm worked by dividing the lumbar spine model into particles with either rigid- or flexible-body characteristics. The vertebrae were represented by rigid bodies, each comprising a rigid coordinate system. The IVDs were represented as a discretized, meshed FEM body containing both rigid and flexible particles: Endplate particles were modeled as rigid-body particles; thus, they remained the same distance from each other over time. The remaining part of the IVD consisted of flexible particles, each experiencing an independent deformation (Figure 2b). The advantage of this approach was the reduction of unknowns, as FEM nodes could be attached to rigid-body particles using Lagrange multipliers. The movement of these attached nodes was then fully prescribed by the rigid body movement of the vertebrae, resulting in a reduction in the size of the equation system. A multimapping step reunited all particles before FEM parameters such as inertia and material properties were applied. The authors state that the model delivers accurate results without penalizing precision. The method was inspired by Stavness et al., who already implemented it in their software, ArtiSynth (Vancouver, BC, Canada) [74]. The software defines deformable bodies by representing their nodes as three DoF particles. Together with the other type of dynamic components—six-DoF rigid bodies the model can be formulated as an ordinary differential Equation (ODE) and solved by a semi-implicit integrator. While particle-based approaches such as ArtiSynth are particularly well-suited for coupled biomechanical simulations, FEM models are approximated by a

lumped mass matrix and a linear co-rotation, and other rotation effects are neglected. Linear FEM representations might not be able to adequately represent large deformations, as undergone by the IVD.

In 2021, Remus et al. [47] published a passive spine model created with ArtiSynth combining rigid vertebral bodies and deformable IVDs as FEM bodies. In Remus' work, data were segmented from the VHM [67] and smoothed. Auxiliary vertebral bodies were additionally derived and acted as interfaces to the IVDs. Vertebral bodies and endplates were not differentiated. Facet joints were modeled in the shape of cylindrical rigid bodies. As FEM components, the IVDs were modeled with a Yeoh material model for the NP and a Mooney–Rivlin material model for the four lamellae of the AF, which was composed with multi-point springs linking the external nodes of the lamellae. Ligaments were modeled as multi-point springs or axial springs. The researchers validated their model in a quasi-static framework with multiple load cases and respective kinematics of in vitro literature data and numerical data. They calculated intradiscal pressure by using the negative mean of normal stresses of all FEM nodes in the NP. Values showed high alignment with in vitro literature data of IVDs at all vertebrae levels. Still, no muscular components were integrated into the model. In subsequent studies, the model was used to study the impact of degenerative changes of the IVD on the axis of rotation by altering its mechanical properties [75] and the effects of a simplified intra-abdominal pressure [76] with integrated muscles as active components. Both studies showed reasonable results, as the authors stated.

Refer to Table 2 for an overview of the reviewed bidirectional co-simulation studies. In summary, two main approaches have been used to couple MBS and FEM spine models bidirectionally: a gluing algorithm providing constant data exchange at certain RPs, managed by a co-simulation engine, and a particle-based approach dividing the model into rigid and flexible particles and solving the resulting ODE with a semi-implicit integrator.

**Figure 2.** Graphical representation of bidirectional coupling methods found in the literature. (**a**) Coupling algorithm implemented by Monteiro et al. [46]. Two RPs, pictured as black dots, served as an interface between the flexible IVD and the rigid vertebrae. (**b**) Coupling method realized by Dicko et al. [73], which was inspired by an algorithm of Stavness et al. [74]. The model was divided into two types of particles. Rigid particles are illustrated as black dots and flexible particles as circles.


**Table 2.** Simulation studies using bidirectional co-simulation to investigate the spine with information on the execution order, transferred data, the software structure and the source of the model geometry.

#### **5. Limitations and Challenges**

Most co-simulation models of the spine use a unidirectional coupling approach, transferring data singularly from one simulation to the other by applying both simulation methods consecutively. MBS spine models can profit by defining joint stiffness parameters based on FEM simulations of the IVD. FEM models of the spine or spinal components can be improved when muscle or ligament forces and moments are implemented as BCs. However, when providing a time-dependent input, the time-dependency is influenced by the deformation properties of a material or component. Equal mechanics cannot be expected when comparing FEM and MBS IVD representations due to their different modeling approaches. The input accomplished by the results of the MBS initially carried out is therefore only partly suitable for being incorporated in a subsequent FEM simulation. This limitation is demonstrated in the effort that has been put into synchronizing the respective models of the unidirectional co-simulations. The manual adaption steps that become necessary may provoke inaccuracies and take much time. Kumaran et al. [64], for example, identified an issue with the interaction of forces and displacement in their simulation and stated that the muscle forces did not produce the desired motion of the vertebrae. They then implemented the simulation with a given displacement rather than a correct interaction of forces and displacement. The same synchronization difficulties were experienced by Khoddam-Khorasani et al. [40] and Liu et al. [41], who both mentioned the need for a manual, iterative process to achieve convergence between the FEM and MBS model or to account for factors such as compressive and shear stiffness. As in Liu et al. [41], the concept of a FL is frequently used to apply summarized loads in the direction of the spinal curvature to provide realistic loading conditions. However, it neither accounts for time-dependent changes in loads, nor dissolves the need of trial and error procedures [36]. In sum, unidirectional co-simulation is often associated with lower accuracy and convergence issues.

To overcome these limitations, a few recent studies have implemented bidirectional co-simulations of the spine. By updating the deformation values of the IVDs and the resulting positions of the vertebrae in every increment, updated reaction moments and forces of muscles, tendons and ligaments can be considered. Of the research groups using bidirectional co-simulation models, all authors found that their models were able to deliver accurate results [46,47,73]. However, two main methods were identified to bidirectionally couple the MBS and FEM solver. Monteiro et al. used an interface approach consisting of two linking RPs, at which the kinematic and kinetic data were exchanged constantly. Although this constant exchange of data solves many of the problems encountered in unidirectional co-simulation, a limiting factor could have been that only one reference node represented the interface between the vertebra and its endplate in this study. The pressure distribution on the IVD thus needed to be derived from one single node, which may have resulted in lowered accuracy. As already reported by Monteiro et al., the implementation of multiple RPs per linkage would be a reasonable adaption to account for more complex deformations.

The particle-based approach divides the whole model into two types of particles, rigid and flexible ones [47,73,74]. Thus, the interface between the flexible FEM and the rigid MBS bodies consists of more than one RP. A mapping step combines the distinct particle sets into one model consisting of a single ODE, which is solved by an integrator. A limiting factor in this approach is the use of a lumped mass matrix and linear co-rotation definitions (neglecting other rotation effects) for the FEM component, which has been associated with less adequate representation of large deformations, as they appear in the IVD. [77]

Despite these limitations, bidirectional co-simulation models of the spine can provide a holistic understanding of the spine because they consider both the overall kinematics with the muscular and gravitational forces and moments, and the detailed mechanics of the IVDs with their deformations and stress distributions.

#### **6. Conclusions and Future Directions**

FEM is broadly used in detailed analyses of internal stresses in deforming bodies but lacks computational efficiency. MBS is more efficient, but can only achieve a certain degree of detail when it comes to deformations of model components. A combination of both methods, not only in a unidirectional way, but in a bidirectional manner of data exchange, can provide both accuracy and efficiency.

Future studies will likely include widespread use of bidirectional co-simulation models to understand and predict the behavior of the spine. Including automated segmentation algorithms such as the one implemented by Sekuboyina et al. [26] could accelerate two things: individual, more adequate diagnoses due to patient-specific geometries, and clinical investigations containing large cohorts. Those detailed, personalized simulations of large cohorts could be used to better understand the underlying mechanisms of pathological changes and the biomechanics of overload situations in ambitious athletes, or to predict injuries before they occur.

**Author Contributions:** Conceptualization, K.N. and J.S.K.; writing—original draft preparation, K.N.; writing—review and editing, T.L., V.S. and J.S.K.; supervision, J.S.K.; project administration, J.S.K.; funding acquisition, J.S.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program, Grant No.: 101045128—iBack-epic—ERC-2021-COG.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** We thank Andreas Zwölfer for the constructive and helpful exchange.

**Conflicts of Interest:** J.S.K. is a Co-Founder of Bonescreen GmbH. All other 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.

#### **Abbreviations**

The following abbreviations are used in this manuscript:

