## 1. Introduction

According to the World Health Organization (WHO), back pain of the lower spine has reached epidemic proportions, being reported by about 80% of people at some time in their life [1]. Hoy notes that lower back pain causes more global disability than any other condition [2]. The causes of back pain are manifold. Reasons may be degenerative changes of intervertebral discs, spinal stenosis, spondylolisthesis, spinal stenosis, traumatic injury, skeletal irregularities, and osteoporosis. But there is also a correlation between increased body weight and low back pain [3–5]. These overloads may result in sequelae such as instability or dysesthesia with narrowing of the spinal canal and compression of the associated nerve roots may be a consequence. If a conservative treatment does not achieve a pain relief anymore, only a surgically induced stabilization like a therapeutic option remains. Mayer [6] indicates that no validated results exist about the success of these operations over a period of two years, and thus provide clear evidence of reduced adjacent segment degeneration. Furthermore, statements about the sustainability of clinical and radiological results with respect to a reconstruction of the segmental mobility are still missing. A modern non-invasive method for assessing the effect of medical operational measures is the load determination in various spinal structures by means of computer modeling.

To capture the load distribution of the lumbar spine, highly developed so-called finite-element (FE) models [7–11] can be used which allow a very detailed calculation of force and torque transmission—both in the rigid bones and in the elastic proportions as intervertebral discs, ligaments, and implants. The complexity of the models, however, requires high computing times for each load case. Is the aspect of the holistic approach of the human movement in the focus of interest, the multibody simulation (MBS) is a suitable simulation method? Due to the highly efficient short computation times, even complex motion sequences can be simulated.

## 2. Basics of model building and simulation

Computer modeling is a widely used method for determining the biomechanical behavior of a system. A distinction is made between the FE modeling and MBS modeling. Both simulation methods are used to aid engineers and researchers in investigating and analyzing mechanical and mechatronic systems in various fields such as the car, railway, engine, aerospace industries, and medical applications. The difference between the FE modeling and MBS modeling lies in the basic model structure, and thus in the field of application. For analyzing highly sophisticated problems, the FE modeling is the appropriate modeling method. Within this mathematical method, based on the numerical solution of partial differential equations, the system is divided into a finite number of simple geometric elements, the finite elements [12]. An FE model is composed of geometric elements starting with points in space, which are linked to each other by so-called edges. Three points define a triangular plane. The corners are called vertices. Connected triangles result in faces which have a contour defined by a sequence of vertices, which is called a polygon. For visualization purposes, the faces can be colored and textured. This information is attached to the vertices, as it is common in computer graphics (**Figure 1**).

For each junction of elements, the so-called nodes, boundary, and transition conditions are defined in accordance with a specific material law [13]. In this way, highly detailed analyses of material properties can be carried out by an FE model. However, FE modeling requires a relatively long computational time. According to Berkley [14], the model's accuracy increases with the number of finite elements that are used to describe the geometry. But each additional element also means additional computational time. Solving such large numbers of FE equations leads to an enormously time-consuming calculation as in [15, 16]. The accuracy of the system and the expected computing time must therefore be carefully matched. A much faster calculation method is the MBS. The MBS describes the interaction of individual bodies with each other and with their environment, that is, the forces and torques acting on the bodies and between the bodies as well as the resulting movements [17].

The modeled bodies are assumed to be rigid and thus not deformable. The different bodies are linked through massless joints and kinematic constraints, which allow certain relative motions and restrict others [12]. Each body has a specific mass and a moment of inertia. This mass and inertia characteristics are associated with a single point, the center of gravity of the body. Thereby, every rigid body is provided with its own coordinate system (reference system). Based on this coordinate system, the exact location in space can be determined with respect to the inertial frame.

The degrees of freedom (DoF) of a rigid body are described by independent coordinates. A rigid body has a maximum of six DoF. It can move along the *x*-, *y*-, and *z*-axis (**Figure 2**, left), and rotate about these axes (**Figure 2**, right).

Possible translations and rotations relative to the inertial system can be realized through the implementation of appropriate joint types. If a model is composed of two or more rigid bodies, they are connected by a massless joint. The joints provide certain DoF. Depending on the definition of joint type, the joint

allows bodies or a body and a reference frame to be connected to each other with no DoF,

provides free movement to the body, with six DoF,

applies kinematic excitations to a body, as a rheonomic joint.

Furthermore, it is possible that forces and torques can be transferred between two body markers, which are generally located on different bodies, by a so-called force element. The type of interaction depends on the implemented force law. According to [19], the DoF of a model are determined by a number of independent state variables, which characterize the movement of the body. Furthermore, the movement of the MBS system depends on the inertial properties of its bodies, masses, inertia tensors, and its centers of gravity. The kinematics of the model is defined by a system of coupled differential equations. These equations of motion give the relation between motion and forces. For this reason, they can be solved in two ways: the forward or direct dynamics and the inverse dynamics. In the case of forward dynamics, the movement is determined by known internal forces or torques, whereas the reconstruction of the internal forces or torques from movements and external forces, like the gravity, is called inverse dynamics [20].

It should be noted that the presented modeling styles MBS and FE can be used independently of each other, but can also be combined with each other.

## 3. MBS lumbar spine models

The aim of our biomechanical computer modeling is to consider the characteristics of a musculoskeletal system through the use of knowledge from the fields of mechanics, anatomy, and physiology in the model in an appropriate manner, in order to obtain as accurately as possible a realistic simulation of the biomechanical behavior of the system. In the following, various application examples of a lumbar spine model that takes all spinal structures with their specific material properties into account are presented.

### 3.1. Structure of the spine model

The MBS models of the lumbar spine consist of vertebrae L1–L5, os sacrum and os ilium. These bony parts are connected by joints with appropriate DoF and ligament structures, which are attached to characteristic points of the skeletal parts. The facet joints are realized as three-dimensional (3D) contact areas, so that the acting contact forces avoid the penetration of two corresponding joint surfaces. All the individual structures are modeled with specific material properties in order to simulate their realistic mechanical behavior.

### 3.2. Surface generation

The 3D surface models of the vertebral bodies are based on computed tomography (CT) data sets of artificial vertebrae of Europeans. These data sets are generated through the use of segmentation and visualization plugins to visualize the data sets and make them available for simulation as Computer Aided Design (CAD) data sets (**Figure 3**) [21, 22].

### 3.3. Modeling of the biomechanical behavior of the spinal structures

The surface models are extended to biomechanical models by adding specific biomechanical properties to all spinal structures. The biomechanical property of the translational movement of the intervertebral discs is described by a force law, which is composed of a geometry-based stiffness term and a damping term [18, 23]. Most movements of the intervertebral discs do not consist of purely translational movements but also of intervertebral rotations. If the disc is deflected from its initial position, a deflection-dependent reaction torque will be developed [24].

As already said, in the model also the anterior and posterior longitudinal ligament (ALL/PLL), the ligamentum flavum (LF), the capsular ligament (CL), the interspinous ligament (ISL), the supraspinous ligament (SSL) as well as the intertransverse ligament (ITL) are implemented. Because a ligament can be defined only between two points, broad ligament structures are realized by a bundle of several fibers. The biomechanical properties of the ligaments are represented by characteristic curves of appropriate literature [25].

In the following models, the facet joints are simulated by a 3D contact surface whose dimensions and orientation are directed to the inclinations of the facet joints of the surface model. In this way, the modeling considers the different individual dimensions and specific orientations of the facet joints in all lumbar levels. Is exclusively a detailed analysis of the load behavior of the facet joints in focus, the facet surfaces are simulated by a 3D layer of cartilage (see Section 4). It should be noted that a detailed modeling of the facet surfaces is accompanied by a large increase of computation time. Therefore, prior to the model building a specific question should be defined, in order to decide whether such a detailed facet joint modeling is needed.

A further important structure to stabilize the human body in upright position is the musculature. According to [7], the four muscle groups, left and right musculus erector spinae and left and right rectus abdominis muscle, are modeled. The force representing the musculus erector spinae is varied so, that additional torque on the upper endplate of the vertebra L1 is no longer necessary to maintain the lumbar spine in a state of equilibrium. The model parameters used for muscles are taken from [7, 26]. These muscular structures are implemented at the moment just in the models of Chapter 4.1. A comparison of the different spine models with the corresponding muscle models is performed (see Section 4).

### 3.4. Validation and sensitivity analysis of the model

For the validation process, there is the difficulty of developing a suitable method, which confirms the accuracy of the modeling. An established method compares the simulation results with results from accepted publications. But it has to be mentioned that not all parameters which may have a significant influence on the result are always published. If such a factual circumstance is known, which is at the same time the model limitation, this has to be considered in the discussion of the model validity.

The model validation was performed by comparing the simulation results with FE results and in-vivo data as found in references [7, 27–29]. A detailed presentation of the validation is to be taken from reference [24]. Another difficulty lies in the selection of suitable input parameters from the literature. The values of these parameters differ significantly in some cases [30]. The validation and modeling cannot be seen as complete, and are successively continued to develop MBS models that get closer to the reality.

## 4. Practical application examples of simulation

To gain an insight into the practical applications of computer modeling in the field of biomechanical modeling of human structures, selected examples of different spinal simulation cases are shown below. In the first examples, the effects of different spinal alignments and obesity in adults and children on the lumbar spine are discussed. Subsequently practical examples of possible use of computer models in medicine will be introduced.

### 4.1. Effects of different spine alignments in the standing position

Due to different physical constitution and daily physical stress, an individual characteristic alignment of the spine is developed in the course of life. Therefore, the double-S shape of the human spine is subject to a wide range of variations [31–34]. To investigate the effects of such different spinal curvatures on the intradiscal pressure, five models of the lumbar spine with different lordosis angles are created (**Figure 4**). In this work, the lordosis angle is defined by the line through the upper endplate of L1 and the endplate of os sacrum.

Upright standing of a normal-weight person is chosen as simulation case because it represents the natural load case. The models are applied with the weight of the upper body. The force application point is located in the center of gravity. This gravitational force causes a deformation of the intervertebral discs, and a corresponding reaction force in the intervertebral discs is built. Therefore, the average intervertebral disc force, calculated from the different values of corresponding functional spinal units (FSUs) of the five models, increases from the lowest FSU to the uppermost FSU. In this context, the term FSU is understood as the smallest physiological motion unit consisting of two adjacent vertebrae, the corresponding intervertebral disc and all adjoining ligament [25].

**Figure 5** shows the intradiscal pressure of all functional spine units of the different aligned models in comparison. Because of the different spine alignments, it could be assumed that also the intradiscal pressures of the various models should greatly differ from each other. But this assumption is not immediately evident from **Figure 5**. However, the percentage difference of the maximum and minimum intradiscal pressure value of an FSU is determined, and the percentage difference of the FSUs L5-Sac and L4-L5 is 8% and of the FSU L2-L3 even 13%. For the intervertebral discs, the resulting percentages of L3-L4 and L1-L2 are 3% and 5% respectively.

To sum up, it can be stated that the individual spinal alignment has an impact on the load distribution of the intervertebral disc.

### 4.2. Effects of overweight on the spinal biomechanics

The WHO called being overweight and obesity the leading chronic health problems [35]. In the world more than 1.9 billion adults, 18 years and older, were overweight in 2014 and of these over 600 million were obese [36]. Overweight and obesity are contributing factors for many ailments and can lead to the development of chronic diseases, like diabetes, cardiovascular disease, coronary heart disease, or osteoarthritis [37–39]. While there are a large number of studies on the effects of these factors on the cardiovascular system and the psyche of the concerned person [40, 41], the potential consequences of orthopedic damage, particularly of the spine, associated with obesity are hardly known yet. In [3–5], a direct correlation between back pain and obesity is described. The exact load changes within the various spinal structures have not been sufficiently explored yet. With the help of computer modeling, the effects of obesity on the kinematics and kinetics of various spine structures can be analyzed. In the following the potential effects of obesity of an adult and an adolescent human are illustrated by specific MBS models exemplary.

Example 1: Simulation of the effects on the spinal structures of body weight of a normal-weight and an overweight adult male

The modeling was carried out in three steps: step 1 is the creation of an MBS model of the lumbar spine, step 2 the creation of a suitable surface models of the two weight classes, including the determination of their anthropometric characteristics, and step 3 the fusion of the MBS lumbar spine model with the surface model [42]. In this example, the force application point is also located in the center of gravity. Taking into account the body weight, body height (1.85 m), sex (male), age (37), and the original two surface models, a normal-weight man (75 kg) and an obese man (127 kg, grade II) are created with the help of an open-source software for the creation of human 3D surfaces. By means of the mass-distribution model of Zatsiorskj [43], the mass distribution of the individual body segments is determined. All other model parameters are identical in both models. By fusing the surface models with the detailed biomechanical model of the lumbar spine, two new simulation models are built up, so that the effects of a normal-weight and obese man on the spinal structures can be analyzed.

Due to the body weight, the intervertebral discs are deformed. In the case of the normal-weight person, in general the intervertebral deformations in all FSUs are very small (**Figure 6**). The cranial located intervertebral discs are more deformed than the caudal intervertebral discs. Comparing the deformation values of the obese with those of normal weight, it can be seen that they are about 1.7 times greater. The validation of the deformation values is difficult, because we are not aware of studies, whose study design corresponds exactly to ours. Brinkmann investigated the fatigue fracture of human lumbar vertebra under cyclic axial compression and presented an example of force versus deformation curve of a specimen with material properties characterized as “tough” [44]. The force versus deformation curve indicates that a force of 500 N causes a deformation of about 0.025 cm and a force of 750 N a deformation of 0.044 cm. Because the simulated acting weight force of the upper body segments is in the case of normal body weight approximately 460 N, and in the case of overweight approximately 780 N and the deformations are in the range of 0.025 cm–0.04 cm (normal weight) and between 0.042 cm and 0.068 cm (obese), it is seen that the magnitude of this values corresponds to those of Brinkmann.

Each intervertebral disc develops a disc force, which depends on the deformation and deformation velocity of the corresponding FSU. **Figure 7** shows the force component acting perpendicular to the disc. It can be clearly seen that the intervertebral discs of the obese are much more stressed than those of normal weight. The intervertebral disc force of the obese is, as well as in the above-discussed case of deformation, 1.7 times higher than the intervertebral disc force of the normal weighted. It could therefore be concluded that the deformation value, as opposed to the deformation velocity (see Section 3.3), has a great impact on the intervertebral disc force.

**Figure 8** shows the intersegmental rotations of the intervertebral discs. Whereas all intervertebral discs of the obese person perform flexions, the intervertebral discs of the normal-weight person are deflected in different directions. The lowest two intervertebral discs of normal-weight person perform flexion. The FSU L4-L3 and L2-L1 rotate backwards and the uppermost FSU L2-L1 rotates forward. Particularly evident is the high rotation values of the obese person. The FSU L1-L2 performs with about 17°, the largest rotational movement. In the literature [45, 46], the maximum range of motion (RoM) of FSU L1-L2 is specified by values between 2 and 13°. Consequently, our simulation values differ from these values. But the difference is relativized by the therein given standard deviations of about 2.5°. Further, it has to be analyzed whether the model of the obese person sufficiently reflects the reality, because identical biomechanical parameters are used in both simulation models. It is necessary to investigate in further work whether the biomechanical properties of the spinal structures are adapted in reality to the body weight, in order to avoid such rotations and prevent overloading or degenerative damage. Moreover, in this model, the stabilizing muscles are not considered, which could prevent such rotation.

The movements of the FSUs can affect the biomechanical behavior of posterior located facet joints (**Figure 9**). It is striking that the loads of the FSU L1-L2 respectively L2-L3 are not loaded. Further it can be seen, that in the case of obesity the loads decrease to cranial. A possible reason may be, in this section, the pure ventral directed rotations. The two corresponding joint surfaces of the facets are distracted, and thus no contact force is build. To achieve a complete clarification of this, a model specification in the form of a sophisticated 3D facet joints modeling is required (see Section 5).

Example 2: Simulation of the effects on the spinal structures of body weight of a normal-weight, an overweight, and an obese child

More than 42 million children under the age of five were overweight in 2013 worldwide [35]. According to Deckelbaum [47], not only in the USA the number of overweight children and adolescents has doubled in the last two to three decades but also the similar doubling rates are being observed worldwide. Therefore, the WHO called overweight and obesity as the leading chronic health problem [35].

Three MBS models of the child's lumbar spine were created to quantify the effects of normal-weight, overweight, and obese children to the kinematics and transmitted forces and torques of the different spinal structures. With the help of these 3D models, dynamic movements and static situations can be simulated. For simulation and load calculation of the effects of different weight classes, the masses of the body segments are adapted to these weight classes [48]. Therefore, the first lumbar vertebra L1 is loaded with the weight of the upper body segments of a normal weight, an overweight, and an obese child in different simulations. The total mass fractions of the segment parts is *m*_{normal} = 22.91 kg for normal-weight child, *m*_{overweight} = 27.73 kg for overweight child, and *m*_{obese} = 32.56 kg for obese child. In the case of overweight and obesity, the intervertebral discs are more heavily loaded in comparison to the normal-weight child. In both cases, the intervertebral discs of the FSU L4-L5 are most loaded (**Figure 10**). The slightest load undergoes the FSU L1-L2. The increase in body weight has a direct influence on the loading structure of the intervertebral discs. If the normal-weight child would become overweight, the loads of the intervertebral discs would rise by an average of 1.5 times. In the case of obesity, the intervertebral discs are more than 2.3 times higher loaded.

Through the increased body weight, the intervertebral discs perform pure flexion (**Figure 11**). Here, the amount of flexion of the obese child is mostly more than doubled compared with the movement of the normal-weight child. Particularly obvious is the difference in the FSU L3-L4. In the case of obesity, this FSU is deflected more than 3.5 times larger from its initial position (normal weight).

Through the forward rotations of the intervertebral discs, the approach and origin points of the posterior ligaments move away from each other. As a result, the corresponding ligaments are stretched and develop in an opposite direction acting as reaction force (**Figure 12**). In the case of an obese child, the posterior ligaments are 2 to 2.5 times more loaded than under normal body weight. Under obesity, a particularly extreme ligament load occurs in the FSU L3-L4. This more than 4 times larger load may result from the above-described extensive intersegmental rotations in these FSU.

It seems to be evident that by the increase of body weight, the loads of the internal structures rise to the same extent. This correlation, which can be predicted even before the simulation, was confirmed by the simplest load case, the upright standing. Further investigations, such as the analysis of load distribution in different structures of an overweight or obese child during everyday activities or during highly dynamic movements, like they occur in sports, can give additional insights into the effects of obesity on the musculoskeletal system.

### 4.3. Application possibilities of biomechanical computer models in medicine

Because the MBS features very short computation times and surgical planning is increasingly computer-based, preoperative simulations can be used to predict the effects of different surgical methods and to identify the best possible surgical option. Furthermore, in future a transfer of topographic and kinematic simulation data of implants in 3D planning and navigation procedures is conceivable. In this process, coordinates of the implants from the computer model and the 3D model data of the spine may be transferred to the navigation system in the appropriate data format. It should be noted that the MBS modeling of the biomechanical properties of the intervertebral disc represents an initial approach. Therefore, the implementation of FE parts is indispensable to simulate patient-specific biomechanical properties.

To demonstrate the medical application possibilities of the simulation, the effect of a spinal fusion under different load cases is exemplified by an appropriate computer model (**Figure 13**). For that the mechanical properties of the MBS model were adjusted, so that no residual movement in FSU L4-L5 is possible. Thus, two models of the lumbar spine were created [49]. As in previous models, the upright position is used as the load case in this simulation, and also an external load of 600 N and 700 N is applied at the center of the vertebral endplate of L1.

Considering only the disc force of the model without fused FSU L5-L4, it can be seen that the pressures in all FSUs increase with higher external force (**Figure 14**). Comparing the pressure of the intervertebral discs of both models, the disc loads of the FSU L1-L2 are almost the same. A deviation results for the other functional units. The pressure in the FSU’s L5-Sac to L2-L3, of the model with fused FSU L4-L5, is lower than without fusion.

In this case, a spinal fusion affects especially the loads of the facet joints. **Figure 15** clearly shows that after fusion, the forces in the facet joints are particularly higher in the upper FSU. Through such an increased load situation degenerative changes can be caused in the facet joints in long term.

Comparing the intersegmental rotation of the discs (**Figure 16**) with and without fused L4-L5 disc, it can be seen that the discs of the fused model perform less deflection in the FSUs L5-Sac to L2-L3. The lower lumbar spine thus has lower mobility than in the non-fused state. This means that after an intervertebral disc stiffening, the lower lumbar spine would have a lower mobility than in a non-fused state.

In summary, it can be stated that a fusion can lead to load redistribution. In this simulation case, especially the posterior structures that are located above the fused FSUs are stressed more than the ventrally located structures. A further example in the field of medicine is studying the effects of different body weights on selected surgical procedures.

## 5. Meaningfulness and limitation of MBS modeling

It is clearly pointed out that modeling has to be understood as an approximation to the reality. In near future, it will probably not be possible that all influencing factors can be simulated completely as in reality. Often, input parameters based on literature data representing the average of a specific cohort. Because each person has a very specific anthropometric and morphological characteristic, the standard deviations of such investigations may be relatively large. For example, the standard deviation for geometric parameters that describe the sagittal alignment of the spine is partially enormous. The deviation from average of the pelvic incidence is partly more than 20% [31], of the pelvic tilt 53% [32] respectively 70% [33], and of the total lumbosacral lordosis 19% [34], etc. But based on sensitivity analysis, the exact influence of individual parameters can be determined. However, because of the complexity of the models, it is very important to know the exact configuration of the model and the limitations of its input parameters. Therefore, the modeling and the validation process should to be understood as an evolving process and will be advanced in future research.

Opportunities for improvement and expansion of the different model components will be summarized in the following:

In the presented models, the biomechanics of the intervertebral disc is defined by a force law, which can be understood as an initial approach. To precise the biomechanical properties of the annulus fibrosis and nucleus pulposus, a 3D hybrid model consisting of MBS and FE units will be built. The aim is to analyze whether the implementation of FE-disc features contributes to a significant gain of knowledge concerning load distribution of the intervertebral discs and the adjacent spinal structures.

Currently, a highly accurate modeled individual cartilaginous contact layer for all facet joints is in development. These individually formed layers of cartilage (**Figure 17**) allow us to do a 3D calculation of the resulting forces.

An extremely important spinal structure that describes the dynamics of the musculoskeletal system is the spinal musculature. Up to now, the muscle groups of the erector spinae and the rectus abdominis implemented in the models are very rudimentary. In focus of current work are corresponding muscle models that describe the dynamics of muscle contraction in an appropriate manner [50, 51].

The aim of our research is the creation of patient-specific 3D hybrid computer models of the human spine for preoperative planning in neurosurgical and orthopedic stabilizing spinal operations. Known and alternative surgical procedures are simulated and a comparison of these possible surgical procedures is presented to determine the biomechanical effects. The simulation is realized highly efficient with short computation times, which will allow a later use in real systems. That means, the 3D coordinates of the optimized implant position can be exported out of the individual models as data sets and can be transferred, to position the implants very accurately, to a navigation system.