General Computational Model for Human Musculoskeletal System of Spine

Abstract

A general computational model of the human lumbar spine and trunk muscles including optimization formulations was provided. For a given condition, the trunk muscle forces could be predicted considering the human physiology including the follower load concept. The feasibility of the solution could be indirectly validated by comparing the compressive force, the shear force, and the joint moment. The presented general computational model and optimization technology can be fundamental tools to understand the control principle of human trunk muscles.

1. Introduction

The human lumbar spine can support large loads during daily activities such as standing, walking, running, and lifting, where the loads are up to several thousand Newtons [1, 2]. However, it has been reported in experimental studies [3, 4] that an intact ligamentous lumbar spine buckled at the load less than 100 N when a load was applied at the superior end in the vertical direction. Although the trunk muscles have been known to play an important role to withstand external loads [5–7], the principle of trunk muscle activation to obtain such load-carrying capacity of the spine has not been elucidated. Recent experimental studies [4, 8] have demonstrated that the load-carrying capacity of the human spine significantly increased as the load applied to the spine was transferred along a path that approximates its curvature, which is called a follower load path originated from the field of mechanical engineering to solve the problems associated with the stability of columns [9, 10] since the 1950s. In the follower load case, a nearly compressive force was produced in the spine with a small shear force. The follower load concept is a possible principle of muscle activation pattern.

It is not easy to directly investigate trunk muscle activations because there have been difficulties in the in vivo measurements of the activated muscle forces, and the responses of the lumbar spine, such as the intradiscal pressure, resultant joint forces and moments at each vertebral joint. Thus, the computational modeling of the human musculoskeletal system is indispensable to predict the forces of the muscles and the responses of the lumbar spine. Several computational models of the lumbar spine and trunk muscles have been developed to estimate the trunk muscle forces [5, 11–21]. Although the follower load concept was considered in [14, 16–21], it is necessary to improve the generality of model to reflect the physiological conditions of the human spine. In this study, a general computational model of the human lumbar spine and trunk muscles including optimization formulations was provided to predict muscle forces based on the follower load. A three-dimensional numerical example was tested to validate the given model.

2. Preliminaries

2.1. Finite Element Model of the Spine and Trunk Muscles

In this paper, the fundamental definitions and notations were based on [17]. A part of the human spine consisting of 𝑁 vertebrae and 𝑀 trunk muscles is considered. Each spinal motion segment consisting of vertebra-intervertebral disc vertebra is modeled as a linear elastic beam element located at the vertebral body centers. Position vector of the 𝑖th vertebral body center is given as a node by 𝐩_𝑖, 𝑖=1,2,…,𝑁. Let 𝐹^𝑚𝑘 and 𝑃𝐶𝑆𝐴_𝑘 be the 𝑘th muscle force and the physiological cross-sectional area of the 𝑘th muscle for 𝑘=1,2,…,𝑀. Assume that there are 𝑀_𝑖 muscles acting on 𝑖th vertebra among 𝑀 trunk muscles and 𝐅^{𝑚𝑖,𝑗}, 𝑗=1,2,…,𝑀_𝑖, denotes the 𝑗th muscle force vector starting from the attachment point in 𝑖th vertebra. Let 𝐩_𝑖,𝑗 be the position vector of the attachment point of 𝑗th muscle acting on 𝑖th vertebra. Geometric data such as vertebral positions and locations of muscle attachment points can be obtained from published anatomical data of the human spine and muscles [11, 22].

2.2. Static Equilibrium Equations

Let us assume that the spinal system is in static equilibrium. The displacements including translations and rotations of each beam element are related with the forces and the moments acting on the vertebral body centers. Let us call these forces and moments motion segment forces and motion segment moments, respectively. The relation between the motion segment forces and the motion segment moments, and the displacements at vertebral nodes, could be defined as⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣𝐅^1𝑚𝑠𝐅⋮_𝑁𝑚𝑠𝐌^1𝑚𝑠𝐌⋮_𝑁𝑚𝑠⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦=𝐊⋅⎡⎢⎢⎢⎢⎢⎢⎣𝐝₁𝐝₂𝐝⋮_𝑁⎤⎥⎥⎥⎥⎥⎥⎦,(2.1) where 𝐅^𝑖𝑚𝑠, 𝐌^𝑖𝑚𝑠, and 𝐝_𝑖, denote the motion segment force, the motion segment moment, and the displacement vectors at 𝑖th vertebral body center, respectively. 𝐊 represents the stiffness matrix describing linear elasticity of the spine model. The stiffness matrix 𝐊 of the motion segment can be obtained from experimental studies such as [23, 24]. The displacement vector 𝐝_𝑖 at 𝑖th node consists of the translation components, 𝑑^{𝑡𝑖,𝑘}, 𝑘=1,2,…,𝐾, and the rotation components, 𝑑^{𝑟𝑖,𝑙}, 𝑙=1,2,…,𝐿, where 𝐾 and 𝐿 are the number of translational and rotational degrees of freedom at each node, respectively.

Then, for given external forces 𝐅^𝑒𝑖 and moments 𝐌^𝑒𝑖 applied at 𝑖th vertebral body center, the static equilibrium equations at the vertebral nodes can be formulated by^𝑀𝑖_𝑗=1𝐅^{𝑚𝑖,𝑗}−𝐅^𝑖𝑚𝑠+𝐅^𝑒𝑖=𝟎,𝑖=1,2,…,𝑁,(2.2)_𝑀𝑖_𝑗=1𝐫_𝑖,𝑗×𝐅^{𝑚𝑖,𝑗}−𝐌^𝑖𝑚𝑠+𝐌^𝑒𝑖=𝟎,𝑖=1,2,…,𝑁,(2.3) where 𝐫_𝑖,𝑗=𝐩_𝑖−𝐩_𝑖,𝑗 for all 𝑖 and 𝑗, represents the moment arm of the muscle force.

2.3. Resultant Joint Force and Resultant Joint Moment

The resultant joint force at each vertebra is the sum of all the muscle forces, the applied external forces, and the force transmitted from the upper vertebra. Hence, the resultant joint force, 𝐅^𝑖𝑗𝑡, at 𝑖th vertebra is calculated iteratively: for 𝑖=1,2,…,𝑁,𝐅^𝑖𝑗𝑡=^𝑀𝑖_𝑗=1𝐅^{𝑚𝑖,𝑗}+𝐅^𝑒𝑖+𝐅^{𝑗𝑡𝑖−1}=𝐅^𝑖𝑚𝑠+𝐅^{𝑗𝑡𝑖−1}(2.4) with 𝐅^0𝑗𝑡=𝟎.

The follower load path direction at each node was defined in order to decompose the resultant joint force into the compressive force and the shear force. Let the compressive force direction vector 𝐜_𝑖 at 𝑖th node be𝐜_𝑖=𝐩^𝑖−𝐩_𝑖−1‖‖𝐩_𝑖−𝐩_𝑖−1‖‖,𝑖=1,2,…,𝑁,(2.5) under the assumption that 𝐩₀=𝟎, which indicates the direction of 𝑖th beam element.

Then, 𝐅^𝑖𝑗𝑡 can be decomposed into two perpendicular compressive forces 𝐅^𝑐𝑖=(𝐅^𝑖𝑗𝑡⋅𝐜_𝑖)𝐜_𝑖 and shear force 𝐅^𝑠𝑖=𝐅^𝑖𝑗𝑡−(𝐅^𝑖𝑗𝑡⋅𝐜_𝑖)𝐜_𝑖 at 𝑖th node for 𝑖=1,2,…,𝑁 (Figure 1) as𝐅^𝑖𝑗𝑡=𝐅^𝑖𝑗𝑡⋅𝐜_𝑖𝐜_𝑖+𝐅^𝑖𝑗𝑡−𝐅^𝑖𝑗𝑡⋅𝐜_𝑖𝐜_𝑖=𝐅^𝑐𝑖+𝐅^𝑠𝑖.(2.6) The resultant joint moment 𝐌^𝑖𝑗𝑡 at 𝑖th node for 𝑖=1,2,…,𝑁 is the same to the motion segment moment 𝐌^𝑖𝑚𝑠.

Figure 1 

Decomposition of joint force 𝐅𝑖𝑗𝑡 into the compressive force 𝐅𝑐𝑖 and the shear force 𝐅𝑠𝑖.

3. Formulation of Optimization Scheme

3.1. Assumptions for Physiology

In this study, the displacement vector 𝐝_𝑖 at 𝑖th vertebral body center for 𝑖=1,2,…,𝑁 and the 𝑘th muscle force 𝐹^𝑚𝑘 for 𝑘=1,2,…,𝑀 were unknowns. Since the number of unknowns is substantially larger than that of equilibrium equations (2.2) and (2.3), an optimization scheme is necessary to predict nodal displacements and muscle forces. To formulate the optimization scheme, requirements from human physiology for the spine must be considered. First, the compressive forces, the shear forces, and the joint moments at vertebral body centers should be minimized in order to avoid injuries or damages to soft tissues such as intervertebral discs or ligaments in the spine region [25]. The square sum of muscle forces and the cubic sum of muscle stresses should be minimized in order to increase the efficiency of muscle activation and to decrease the fatigue of muscles, respectively [11, 12], where the muscle stress is defined by the ratio of the muscle force to the physiological cross-sectional area of muscle. Finally, the follower load concept to minimize the shear force in comparison to the compressive force at each vertebra should be considered [25]. These multiple issues can be formulated in a multiobjective cost function as𝑓𝐝₁,…,𝐝_𝑁,𝐹^𝑚1,…,𝐹^𝑚𝑀=𝑤_1𝑁_𝑖=1‖‖𝐅^𝑐𝑖‖‖²+𝑤_2𝑁_𝑖=1‖‖𝐅^𝑠𝑖‖‖²+𝑤_3𝑁_𝑖=1‖‖𝐌^𝑖𝑗𝑡‖‖²+𝑤_4𝑀_𝑘=1||𝐹^𝑚𝑘||²+𝑤_5𝑀_𝑘=1||||𝐹^𝑚𝑘𝑃𝐶𝑆𝐴_𝑘||||³,(3.1) where 𝑤₁,𝑤₂,𝑤₃,𝑤₄, and 𝑤₅ are weight factors.

To find the relevant solution of the given optimization problem, the weight factors must be selected based on the quantitative relationships between physiological characteristics. However, there is little information regarding the quantitative relationships due to difficulties in experimental measurement of such in vivo characteristics. Thus, it is necessary to reduce the number of terms in the objective function under feasible assumptions. Since (2.2), (2.3), and (2.4) indicate that the resultant joint forces and moments are dependent on the muscle forces, and the square sum of muscle forces and the cubic sum of muscle stresses are calculated from the muscle forces, the objective function can be modified as𝑓𝐝₁,…,𝐝_𝑁,𝐹^𝑚1,…,𝐹^𝑚𝑀=𝑤_1𝑁_𝑖=1‖‖𝐅^c_𝑖‖‖²+𝑤_2𝑁_𝑖=1‖‖𝐅^𝑠𝑖‖‖²+𝑤_3𝑁_𝑖=1‖‖𝐌^𝑖𝑗𝑡‖‖²,(3.2)𝑓𝐝₁,…,𝐝_𝑁,𝐹^𝑚1,…,𝐹^𝑚𝑀=^𝑀_𝑘=1||𝐹^𝑚𝑘||²,(3.3) or𝑓𝐝₁,…,𝐝_𝑁,𝐹^𝑚1,…,𝐹^𝑚𝑀=^𝑀_𝑘=1||||𝐹^𝑚𝑘𝑃𝐶𝑆𝐴_𝑘||||³.(3.4) In addition, the follower load concept can be formulated by a constraint as‖‖𝐅^𝑐𝑖‖‖≤𝛼‖‖𝐅^𝑠𝑖‖‖,(3.5) where 𝛼 is a restriction coefficient for the follower load concept. The physiologically feasible upperbounds for the displacements of vertebrae and the muscle forces must also be presented.

The optimization scheme can then be formulated as follows.

Minimize 𝑓𝐝₁,…,𝐝_𝑁,𝐹^𝑚1,…,𝐹^𝑚𝑀(3.6)

s. t.(1)𝐅^𝑚−𝐊⋅𝐝+𝐅^𝑒=𝟎, where 𝐅^𝑚=⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣^𝑀1𝑗=1𝐅^𝑚1,𝑗⋮_{𝑀𝑁𝑗=1}𝐅^{𝑚𝑁,𝑗}^𝑀1𝑗=1𝐫_1,𝑗×𝐅^𝑚1,𝑗_{𝑀𝑁𝑗=1}⋮𝐫_𝑁,𝑗×𝐅^{𝑚𝑁,𝑗}⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦,𝐝=⎡⎢⎢⎢⎢⎢⎢⎢⎢⎣𝐝₁𝐝⋮⋮⋮_𝑁⎤⎥⎥⎥⎥⎥⎥⎥⎥⎦,𝐅_𝑒=⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣𝐅^𝑒1𝐅⋮_𝑒𝑁𝐌^𝑒1𝐌⋮_𝑒𝑁⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦,(3.7) and 𝐊 is the stiffness matrix defined in (2.1);(2)‖𝐅^𝑐𝑖‖≤𝛼‖𝐅^𝑠𝑖‖,𝑖=1,2,…,𝑁; (3)|𝐹^𝑚𝑘/𝑃𝐶𝑆𝐴_𝑘|≤𝜎, 𝑘=1,2,…,𝑀 where 𝜎 is the maximum muscle stress;(4)0≤|𝑑^{𝑡𝑖,𝑘}|≤𝑑^{𝑡𝑖,𝑘,max} and 0≤|𝑑^{𝑟𝑖,𝑙}|≤𝑑^{𝑟𝑖,𝑙,max}, 𝑖=1,2,…,𝑁 where 𝑑^{𝑡𝑖,𝑘,max} and 𝑑^{𝑟𝑖,𝑙,max} denote the upper bounds of 𝑘th translation component and of 𝑙th rotation component of 𝐝_𝑖.

4. Numerical Tests

A three-dimensional problem of the spine from T12 to S1 is tested to confirm the developed computational model and the formulation of the optimization scheme predicting the muscle forces (𝑁=7). Here, 117 pairs of trunk muscles were considered (𝑀=234): 5 longissimus pars lumborum, 4 iliocostalis pars lumborum, 12 longissimus pars thoracis, 8 iliocostalis pars thoracis, 11 psoas, 5 quadratus lumborum, 6 external oblique, 6 internal oblique, 1 rectus abdominus, 12 thoracic multifidus, 20 lumbar multifidus, 6 interspinales, 10 intertransversarii, and 11 rotatores. The anatomical data at the initial position of the vertebrae, muscle attachments, and physiological cross-sectional areas were obtained from the literature and medical images [11, 19–22]. The stiffness matrix 𝐊 was obtained from previous experimental studies [23, 24].

In this test, (3.2) was used for the objective function. The weight factors 𝑤₁,𝑤₂, and 𝑤₃ are supposed to be 3, 3, and 1, respectively, since 3 N of force and 1 Nmm of moment are considered equally based on the presumed safe limits of intervertebral loads being approximately 3000 N for forces and 9000 Nmm for moments as shown in [12]. The restriction coefficient 𝛼 was selected to be 0.25 based on [20] and the maximum muscle stress 𝜎 was assumed to be 0.46 MPa based on [26]. The upperbounds of all translation component and rotation component were 20.0 mm and 10.0°, respectively. An upright standing posture was considered for the external loading as 300 N of the upper body weight, 3 Nm of the resulting flexion moment applied to T12, and a vertebral weight of 10 N was added to each lumbar vertebra from L1–L5.

The muscle force distribution satisfying the formulated optimization problem was obtained using MATLAB (MathWorks Inc., USA). The number of activated muscles according to the ratio of muscle force to maximum muscle force was summarized in Table 1. The maximum compressive force and the maximum shear force were 691.1 N and 172.8 N while the maximum joint moment was 2271 Nmm. The previous in vivo studies reported that the maximum compressive force, shear force, and joint moment were about 650 N, 190 N, and 8400 Nmm, respectively, in the upright standing posture [1, 14, 15, 27]. The validity of our results seems to be indirectly achieved since the models in [1, 14, 15, 27] were not exactly same to our model.

Table 1 

The number of activated muscles according to the ratio of muscle force to maximum muscle force.

5. Conclusion

In this study, a general computational model of the human lumbar spine and trunk muscles including optimization formulations was provided. For a given condition, the trunk muscle forces could be predicted. The feasibility of the solution could be indirectly validated by comparing the compressive force, the shear force, and the joint moment. The presented general computational model and optimization technology can be fundamental tools to understand the control principle of human trunk muscles.

Acknowledgments

This research was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (2010-0005167), and 2009 National Agenda Project (NAP) funded by Korea Research Council of Fundamental Science & Technology (P-09-JC-LU63-C01).