1 Introduction

As low back injuries remain to be the most frequently occurring injury in occupation-related settings, a comprehensive understanding of spinal loads, muscle forces and the mechanisms behind maintaining trunk stability is essential to mitigate low back injury.

Early biomechanical research for simulating lumbar spine loading in lifting activities resemble basic static models in two dimensions [20, 21]. Recent studies have developed 3D models using dynamic analysis and more detailed musculoskeletal structures of the torso [5, 15]. However, as models became more complex, the number of unknowns significantly exceeds the known parameters. Thus, to resolve the redundancy problem and determine individual muscle forces, different musculoskeletal modelling approaches have been proposed including equivalent muscle models [17], optimization models [7, 15], electromyography (EMG)-assisted models [27], finite element models [26], and multi-body models [22]. However, these models are only supported by commercially available software or customized programs. Furthermore, most of the models are generic models with limited DOFs, which lack the capability of simulating subject-specific motions performed by individuals.

OpenSim is an open-source musculoskeletal modelling software [10]. A detailed lumbar spine model, including the pelvis, sacrum, lumbar vertebrae, and torso, in OpenSim was firstly developed by Christophy et al. [9]. More recently, a thoracolumbar spine model, including the torso and upper limbs, was developed in OpenSim by Bruno et al. [7]. However, this model was only used to simulate static and isometric conditions.

Therefore, the purpose of the current study was to develop an enhanced OpenSim model with whole-body structures including detailed musculoskeletal representations of the trunk, upper, and lower extremities. The model was evaluated by analysing kinematics and kinetics under dynamic lifting conditions.

2 Methods

2.1 Model Development

A 3D multi-segment, musculoskeletal, whole-body model was developed, including 49 DOFs and 258 muscle-tendon units. The model consists of components for the lower limbs, trunk, and upper limbs.

Lower Limbs.

The lower limbs, containing rigid bodies defining the pelvis, thighs, shanks and feet, were sourced from a full-body model (Gait2354) [14] created in OpenSim (Stanford University, Stanford, US). The skeletal portion of the model is modelled as a set of rigid bodies with interconnecting joints. The shank consists of bone geometries for the tibia and fibula, and the foot contains geometries for the calcaneus, navicular, cuboid, cuneiforms and metatarsals. The kinematic relationship between lumbar 5 and the sacrum is based on a study by Anderson and Pandy [2], while the hip joints are represented by a ball and socket joint. The knee joint is based on an earlier model by Delp et al. [11] which formulated transformations for the femur, tibia and patella as functions of the knee angle based on a simple, planar model of the knee in Yamaguchi and Zajac [29]. The pelvic frame of reference is located at the midpoint between the left and right anterior superior iliac spines. In the default position, the pelvic frame of reference aligns with the global coordinate system. It has identified that in the neutral position the pelvis tilts at approximately 13° [6, 19, 28]. This tends to lead to an off-set between the hip flexion angles when comparing experimental data to results found in the literature. In order to simplify this problem, the orientation of the pelvis was modified in the default position such that it exhibits a 13 degree tilt. The lower limbs consist of 48 musculotendon actuators representing 40 muscles. Details of muscles are described in the muscle model part.

Lumbar and Torso.

The lumbar spine model is modified from a base lumbar model developed by Christophy et al. [9]. The model consists of rigid bodies for the pelvis, sacrum, lumbar vertebrae, and torso consisting of bone geometries for the thoracic spine and ribcage.

The motion of lumbar bodies is represented by six DOFs: flexion-extension, lateral bending, axial rotation, and three translational directions. Each lumbar body consists of an independent coordinate system which contributes to the overall motion of the spine. The net bending of the trunk around three orthogonal directions are distributed throughout the lumbar bodies by describing their individual range of motion as a fraction of the net motion. That is, the rotation of each of the lumbar vertebrae is assumed to be linear functions of the trunk flexion–extension, axial rotation, and lateral bending [9]. The coefficients of the linear functions were obtained through a kinematic study of the lumbar spine during lifting [1].

The dimensions of the base lumbar model were modified to match the dimensions of the torso in the lower limb model. The scale factors were calculated using landmarks in the anterior-posterior, superior-inferior and medial-lateral directions. The distance between the xiphoid process and the spinous process of the 9th thoracic vertebra were calculated for the anterior-posterior direction. The distance between the spinous process of the 1st thoracic vertebra and the spinous process of the 5th lumbar vertebra was used for the superior-inferior direction. Markers were placed on the most lateral points of the 6th rib for the medial-lateral direction.

The scaled lumbar model was incorporated into the lower limb model. The distance between the 5th lumbar vertebra and the sacrum were determined by the ratio of the vertebra height to the height of the intervertebral disc. In Gilad and Nissan [13], the relative height index was found using four measurements, as illustrated in Fig. 1, where b and d measure. The relative height index was calculated by dividing the average of the posterior and anterior measurements as defined in the Eq. 1.

Fig. 1.
figure 1

Measurements for determining the relative height index.

Full size image
$$ Relative\,Height\,Index = \frac{(g + h)}{(b + d)} $$
(1)

The vertebral height and g and h measure the intervertebral spacing.

Gilad and Nissan [13] reported a mean height index of 0.402 for the 5th lumbar vertebra. Using the positions of the markers in the OpenSim model, g had a value of 0.020324 and h had a value of 0.021. Thus, by taking the average of intervertebral spacing the distance was estimated to be 0.0083 m. Once the position of the 5th lumbar vertebra was correctly defined, all other bodies in the lumbar model, including the rest of the lumbar vertebrae and the upper torso, were added to the new model.

Upper Limbs.

Upper limbs were included based on a modified and expanded gait model similar to Gait2354 [14]. In our model, the muscles in the upper limbs were excluded and the components purely consist of skeletal structures because the complexity of these muscles makes it difficult to be modelled reliably.

Muscle Model.

The musculature includes the main muscle groups of the lower limb and trunk consisting of 258 muscle fascicles. The main muscle groups included in the model are the erector spinae, rectus abdominis, internal obliques, external obliques, quadratus lumborum, multifidus, psoas major, iliacus, rectus femoris, vastus intermedius, sartorius, gracilis, adductor magnus, pectineus, gluteus medius/maximus, tensor fasciae latae, quadratus femoris, fixme gem, piriformis, biceps femoris, medial gastrocnemius, soleus, and tibialis anterior/posteriors (Fig. 2). The latissimus dorsi was excluded because it has less contribution of trunk stabilisation [9].

Fig. 2.
figure 2

Anterior, mediolateral, and posterior views of the model.

Full size image

The muscle models in OpenSim are based on the Hill-type model [16]. An updated Hill-type muscle model was included in the developed model [25]. Thelen’s muscle model allows muscles to be fully defined by physiological parameters of the muscle architecture. In order to prevent the model from encountering singularities, it now allows users to define a lower bound for muscle activation as well as maximum pennation angles.

Changes were also made such that the fibre length does not reach a lower bound and that energy is conserved throughout the simulation [25].

Anthropometric Properties.

The anthropometric properties of the base models were derived from different sample populations and therefore reflect the characteristics of different subjects. This means that mass and inertial properties are inconsistent throughout the new model. A solution to this problem was to redistribute the overall mass of the generic model to each component in the newly developed model [3].

2.2 Model Evaluation

To evaluate the capacity of the model for estimation of spinal loading, an experiment was performed under lifting tasks. The experiment protocol was approved by University of Auckland’s Human Research Ethics committee.

Customized Marker-Set.

Amarker-set was formulated, mapping the marker positions during motion-capture. As there are currently no standard marker-sets available for researching lifting mechanics, a marker-set based on the Cleveland Clinic marker-set [24] was designed in the current study.

Experimental Protocol.

Eight Vicon MX infrared cameras (Oxford Metric, Oxford, UK) were used to collect the marker positions with a sampling frequency of 100 Hz. Two force platforms (Bertec Corporation, Worthington, Ohio, USA) with a sampling rate of 1,000 Hz were used to measure ground reaction forces. The experiment was conducted on eight healthy, young male adults (age: 21 ± 2.49 years, height: 1.74 ± 0.08 m, weight: 70.76 ± 7.56 kg). Reflective markers of 20 mm in diameter for the clavicle and sternum and 14 mm for all other anatomical landmarks were used. The participants were asked to lift a loaded crate from ground level to a table placed directly in front of them. The crate had dimensions of 40 × 35 × 25 cm (width × length × height) and was loaded to make a total weight of either 7 kg or 12 kg. The table had a height of 72 cm and was positioned at a distance of 40 cm from the participants. The participants carried out each task at a self-selected pace, yet all were completed within two seconds. Participants used a squatted lifting technique with some trunk flexion.

Post-Processing.

Using the tools available in OpenSim, the biomechanical model was scaled to match the anthropometry of each participant based on measured marker positions. Kinematic analysis was then performed to obtain the positions and orientations of the body segments. Then inverse dynamics and static optimisation analysis were carried out to obtain the muscle forces and joint torques. Subsequently, joint reaction analysis enabled calculation of joint reaction forces at each lumbar vertebra. The compression forces and shear forces of each lumbar joint were examined for the lifting motion trials. The outcomes were presented for peak loads averaged over three repeated trials for each of the eight participants. A paired t-test was performed to compare the difference between the two lifting conditions.

3 Results and Discussion

The L5/S1 and L4/L5 joints exhibited the largest forces with very similar values, whereas the L3/L4 joint showed the smallest forces during lifting motions. The largest of the L5/S1 joint may associate with disc degeneration [8]. Very similar trends were observed for the compressive and anterior shear forces (Fig. 3); the first peak occurring immediately after lifting off and the second peak occurring when the crate was lowering on the table. Shear forces change directions, from posterior to anterior, at the L3/L4 joint. This is a result of the spine’s lordosis – the natural curve of the spine. Shear forces act in the anterior direction in lower vertebrae and in the posterior direction for upper vertebrae.

Fig. 3.
figure 3

3D lumbar spinal force patterns after time-normalised for each lumbar level throughout lifting motion with two lifting weights. Positive values indicate anterior and right lateral direction.

Full size image

These results have shown trends similar to what is expected from basic knowledge of lumbar spine loading. Thus, it encourages the feasibility of applying the developed model to research in lumbar spine loading and future opportunities to investigate low back pain and injuries. When lifting a heavier weight, the peak compressive and anterior and lateral shear forces showed higher lumbar spinal loading than with a 7 kg weight for most joints (Fig. 4).

Fig. 4.
figure 4

Peak 3D lumbar spinal forces for each lumbar level throughout lifting motion. A paired t-test was performed. * indicates a significant (p < 0.05) difference between the two lifting conditions.

Full size image

The current results were generally in accordance with previous literature [3,4,5, 12]. Recently this model has been adopted other experimental study and showed good agreement with previous reported lumbar loading and trunk muscle forces [18]. We also compared the simulation results with reported measured intradiscal compressive forces [23] and showed a good match.

One of the major assumptions involved in the developed model is the approximation of human structures into rigid segments. For example, the kinematics of the thoracic and cervical spine was not considered and was, instead, modelled as a single rigid structure. Adding kinematics to these structures can provide more accurate representations of vertebral movement. Another one is the assumption of simple pivotal joints. Although translational motions of the spine were included, there are still mechanisms of the joint that have not been modelled. In addition, passive structures provide friction and tension during extreme flexion and extension motions. Future work could be done to improve the model by incorporating the passive structures. In addition, we assume that the rotation of each of the lumbar vertebrae is linear functions of the lumbar flexion–extension, axial rotation, and lateral bending. Although this approach may introduce some uncertainties of the kinematic data, it’s more uncertain to place small reflective markers on processes of each lumbar processes as the skin movement artefacts will make the measurement unreliable during dynamic movement.

Although the model presented in this study has been designed for the purpose of investigating low back loading, the generality of the model means that it can also be applied to other research questions. Motions that involve whole-body dynamic movement can be simulated and analysed using this model. An example would be gait analysis for patients with cerebral palsy, which includes much more trunk movement than gait for those not affected by cerebral palsy. Due to the accessibility of the model, the opportunities for its application and development are endless.

4 Conclusion

This study provides an advanced and accessible tool for investigating lumbar spine loading. The developed model will support research involved in predicting risk levels for a variety of manual handling tasks. The developed model includes detailed musculoskeletal structures. Methods for performing experiments for occupational lifting tasks have been presented. The results from the lifting experiment are similar to the results from some of the literature. Further validation is needed to verify the accuracy of the model when applied to a specific research question.