**1. Introduction**

When humans evolved to adopt an upright position, the spine became a central structure of the human biomechanical system. As such, it can be a source of pain that can significantly impact a person's quality of life. Intervertebral disc (IVD) degeneration and herniation are possible causes of back pain. Research has shown that both aging [1] and overload [2,3], often experienced by ambitious athletes, can contribute to the degeneration process. However, the relationship between abnormal loading and degeneration is not fully understood. In addition, IVD changes that are visible in magnetic resonance imaging (MRI) have not proven to be clear evidence of back pain [4]. To better understand the mechanisms underlying back pain and identify potential solutions, it is important to study the spine and its related structures in a holistic manner, considering the interactions of motion or posture, pain, and IVD biomechanics. However, the difficulty of directly observing these structures in vivo makes this task challenging. Numerical simulations can be useful in modeling and analyzing the mechanics of the spine, to identify factors that contribute to spinal health problems and potential interventions.

Multibody simulation (MBS) is widely used to gain insights into the healthy and pathological biomechanics of the spine from a macroscopic perspective. MBS models

**Citation:** Nispel, K.; Lerchl, T.; Senner, V.; Kirschke, J.S. Recent Advances in Coupled MBS and FEM Models of the Spine—A Review. *Bioengineering* **2023**, *10*, 315. https://doi.org/10.3390/ bioengineering10030315

Academic Editors: Christina Zong-Hao Ma, Zhengrong Li, Chen He, Aurélien Courvoisier and Elena A. Jones

Received: 31 December 2022 Revised: 13 February 2023 Accepted: 22 February 2023 Published: 1 March 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

study joint reaction forces, muscle forces and muscle activation patterns in combination with respective movements or positions. There are several such models that include intervertebral joints with three rotational [5–19] or six translational [20–22] degrees of freedom (DoFs) while neglecting passive structures, such as ligaments and facet joints. Using force elements as additions to joints allows individual stiffness definitions for all given DoFs [15–20]. For all models, including rotational joints in at least one dimension, the position of joints is a central issue, as it defines the centers of rotation and therefore, directly influences joint kinematics. The fixed centers of rotation in these studies were set either in the center of the IVD [6,12,17], or according to Pearcy and Bogduk [23], in the posterior half of the upper endplate of the inferior vertebra of each motion segment [9,15,16,24].

Finite element method (FEM) spine models aim to simulate deforming bodies in a detailed manner to reveal their inner stresses. Models may include the whole spine, a segment or only an IVD. Most models contain manually created geometries for IVDs and vertebrae [25], but recent approaches use automated segmentation algorithms to derive patient-specific geometries [26,27]. For IVDs, the current gold standard is a biphasic, poroelastic and nonlinear model including an anulus fibrosus (AF) and a nucleus pulposus (NP) component. While bone is usually modeled with linear material properties [28], or as a rigid body [29], material properties for the IVD commonly include hyperelastic material models, such as the Neo-Hooke and Mooney–Rivlin models [30]. As IVD degeneration is believed to affect IVD biomechanics [31,32], several approaches also deal with degenerationdependent material properties [33,34]. A key evaluation criterion of FEM spine models is their ability to describe arbitrary material behavior of the IVD while numerically reporting all its biomechanical functions, such as load transfer and bulging events [35].

In summary, MBS models are suitable for investigations of the overall kinematics and kinetics of the trunk, or even the full body. However, detailed information on flexible body deformation and internal stress cannot be provided. Complex, heterogeneous structures, such as the IVD, or even entire functional spinal units (FSU) combining all stabilizing structures (IVD, ligaments, facet joints), are reduced to a resultant mechanical response to external deformation, often in reference to an approximated center of rotation. In contrast, FEM models are beneficial for simulating detailed structures and analyzing inner stresses and deformations. However, detailed FEM meshes may cause large computational costs, and thus, models often include simplifications. Additionally, defining realistic boundary conditions (BC) for FEM models remains a main challenge due to the lack of possibilities for in vivo measurements. Until now, no satisfactory, broad investigations are available containing this data [22].

Coupling MBS and FEM takes advantage of the strengths of each method and offers great potential to analyze the deformations and stresses of IVDs while considering the overall kinematics and kinetics of the trunk. How coupling is implemented can be subdivided into two approaches: One is to calculate acting forces or present deformations completely before starting the other simulation (unidirectional data flow). The second approach considers acting forces or displacements anew in every increment of an ongoing simulation, generating a constant, bidirectional data flow. Unidirectional coupling commonly includes an inverse dynamic MBS simulation with resulting forces and moments, which are subsequently used as loading and BCs in an FEM model [36–43]. Less common approaches use FEM model displacements and stresses to define properties of force elements, occasionally called bushing components, in MBS simulations [44]. Few approaches apply order-reduction techniques to the FEM model before integrating it into an MBS model [45], which reduces computational costs in linear models, but does not support parameter variation after the reduction is carried out. In general, unidirectional coupling is suitable for linearly deforming FEM bodies because the positions of the body's particles do not change substantially, providing the possibility of directly applying the calculated forces of the MBS simulation subsequently to the respective locations in the FEM bodies. Bidirectional coupling is preferential whenever large deformations are present [46,47]. In that

case, the deformations of the flexible bodies have a direct impact on multibody kinematics, hindering the option of consecutive simulations.

Hence, by using co-simulation, we can examine the pathological spine in a holistic way to analyze the interactions of various factors, such as motion, posture, back pain and the biomechanics of the IVD. In this work, we review advances in coupled MBS and FEM simulations of the spine. Specifically, we excluded reduced order models and instead focused on the type of data exchange, whether unidirectional or bidirectional, as this is essential for the simulation of large deformations and constant interactions of different components in the spine.

#### **2. Methods**

Coupling MBS and FEM models has become a common practice in spine modeling during the last decade. However, no terminology has been established yet to differentiate between what are here called unidirectional and bidirectional co-simulation. To achieve comprehensive coverage of studies using co-simulation, the following keywords were searched in several combinations in different databases: coupled, spine, multibody, MBS, FEM, hybrid, finite element, co-simulation, combined. Resulting articles were manually reviewed and included only if they: (i) contained a type of coupling between MBS and FEM and (ii) were published after 2009 or included models that were still used in studies published after 2009.

#### **3. Unidirectional Co-Simulation of the Spine**

We first reviewed MBS and FEM models of the spine using unidirectional coupling, which was characterized by the execution of one simulation and the integration of the results into another simulation. Considering the order of execution, three distinct methods of coupling were identified: MBS execution and transfer of data to an FEM simulation; FEM and transfer of data to MBS; and a threefold data transfer from FEM to MBS and to FEM again. Figure 1 visualizes the representative simulation models with the different coupling methods.

**Figure 1.** Schematic overview of unidirectional coupling methods. Arrows on the bottom indicate the respective coupling methods. (**a**) Representative FEM model of a FSU containing a detailed IVD, as found in Karajan et al. [44]. Resulting FEM displacements were converted to MBS bushing element parameters. (**b**) Typical MBS model with a joint representing the IVD. Results from the simulation, such as muscle forces or time-dependent displacements of vertebrae or joints, were subsequently included in an FEM simulation. (**c**) Representative FEM model of a FSU or a larger section of the spine. With the BCs defined according to the results of (**b**), these models were often used to calculate IDPs, stress distributions in the IVDs and load-sharing mechanisms.
