Journal of Applied MathematicsVolume 2012 (2012), Article ID 484759, 8 pagesdoi:10.1155/2012/484759

Research Article

General Computational Model for Human Musculoskeletal System of Spine

Kyungsoo Kim,¹ Yoon Hyuk Kim,2,3 and SuKyoung Lee⁴

¹Department of Mathematics, Kyonggi University, Suwon 443-760, Republic of Korea²Department of Mechanical Engineering, Kyung Hee University, Yongin 446-701, Republic of Korea³e-Spine Center, Kyung Hee University, Yongin 446-701, Republic of Korea⁴Department of Computer Science, Yonsei University, Seoul 120-749, Republic of Korea

Received 30 September 2011; Accepted 14 November 2011

Academic Editor: Chang-Hwan Im

Copyright © 2012 Kyungsoo Kim et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 , . Let and be the th muscle force and the physiological cross-sectional area of the th muscle for . Assume that there are muscles acting on th vertebra among trunk muscles and , , 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 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, , , and the rotation components, , , 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 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 , with .

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 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 (Figure 1) as The resultant joint moment at th node for 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 and the th muscle force for 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 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 or In addition, the follower load concept can be formulated by a constraint as 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

s. t.(1), where and is the stiffness matrix defined in (2.1);(2); (3), where is the maximum muscle stress;(4) and , where and 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 (). Here, 117 pairs of trunk muscles were considered (): 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.

Ratio of muscle force to maximum muscle force
Number of activated muscles

0%–20%
7

20%–40%
7

40%–60%
6

60%–80%
2

80%–100%
20

Total

42

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).

References

A. Nachemson, “Lumbar intradiscal pressure,” in The Lumbar Spine and Back Pain, M. I. V. Jayson, Ed., chapter 9, pp. 191–203, 1987.

A. Schultz, “Loads on the lumbar spine,” in The Lumbar Spine and Back Pain, M. I. V. Jayson, Ed., chapter 10, pp. 205–214, 1987.

J. J. Crisco, M. M. Panjabi, I. Yamamoto, and T. R. Oxland, “Euler stability of the human ligamentous lumbar spine. Part II: experiment,” Clinical Biomechanics, vol. 7, no. 1, pp. 27–32, 1992.

A. G. Patwardhan, R. M. Havey, K. P. Meade, B. Lee, and B. Dunlap, “A follower load increases the load-carrying capacity of the lumbar spine in compression,” Spine, vol. 24, no. 10, pp. 1003–1009, 1999.

J. J. Crisco and M. M. Panjabi, “The intersegmental and multisegmental muscles of the lumbar spine: a biomechanical model comparing lateral stabilizing potential,” Spine, vol. 16, no. 7, pp. 793–799, 1991.

H. J. Wilke, S. Wolf, L. E. Claes, M. Arand, A. Wiesend, and T. Bendix, “Stability increase of the lumbar spine with different muscle groups: a biomechanical in vitro study,” Spine, vol. 20, no. 2, pp. 192–198, 1995.

H. J. Wilke, A. Rohlmann, S. Neller, F. Graichen, L. Claes, and G. Bergmannt, “A Novel approach to determine trunk muscle forces during flexion and extension: a comparison of data from an in vitro Experiment and in vivo measurements,” Spine, vol. 28, no. 23, pp. 2585–2593, 2003.

A. G. Patwardhan, R. M. Havey, A. J. Ghanayem et al., “Load-carrying capacity of the human cervical spine in compression is increased under a follower load,” Spine, vol. 25, no. 12, pp. 1548–1554, 2000.

M. Beck, “Die Knicklast des einseitig eingespannten, tangential gedrückten Stabes,” Zeitschrift für Angewandte Mathematik und Physik, vol. 3, no. 6, pp. 476–477, 1952.

G. Herrmann, “Stability of equilibrium of elastic systems subjected to nonconservative forces,” Applied Mechanics Reviews, vol. 20, pp. 103–108, 1967.

I. A. Stokes and M. Gardner-Morse, “Lumbar spine maximum efforts and muscle recruitment patterns predicted by a model with multijoint muscles and joints with stiffness,” Journal of Biomechanics, vol. 28, no. 2, pp. 173–186, 1995.

I. A. F. Stokes and M. Gardner-Morse, “Lumbar spinal muscle activation synergies predicted by multi-criteria cost function,” Journal of Biomechanics, vol. 34, no. 6, pp. 733–740, 2001.

A. G. Patwardhan, K. P. Meade, and B. Lee, “A frontal plane model of the lumbar spine subjected to a follower load: implications for the role of muscles,” Journal of Biomechanical Engineering, vol. 123, no. 3, pp. 212–217, 2001.

A. Shirazi-Adl and M. Parnianpour, “Load-bearing and stress analysis of the human spine under a novel wrapping compression loading,” Clinical Biomechanics, vol. 15, no. 10, pp. 718–725, 2000.

N. Arjmand and A. Shirazi-Adl, “Model and in vivo studies on human trunk load partitioning and stability in isometric forward flexions,” Journal of Biomechanics, vol. 39, no. 3, pp. 510–521, 2006.

Y. H. Kim and K. Kim, “Numerical analysis on quantitative role of trunk muscles in spinal stabilization,” JSME International Journal Series C, vol. 47, no. 4, pp. 1062–1069, 2004.

Y. H. Kim and K. Kim, “Computational modeling of spine and trunk muscles subjected to follower force,” Journal of Mechanical Science and Technology, vol. 21, no. 4, pp. 568–574, 2007.

K. Kim, Y. H. Kim, and S. K. Lee, “Increase of load-carrying capacity under follower load generated by trunk muscles in lumbar spine,” Journal of Engineering in Medicine, vol. 221, no. 3, pp. 229–235, 2007.

K. Kim and Y. H. Kim, “Role of trunk muscles in generating follower load in the lumbar spine of neutral standing posture,” Journal of Biomechanical Engineering, vol. 130, no. 4, Article ID 041005, 2008.

K. Kim, Y. H. Kim, and S. Lee, “Shear force allowance in lumbar spine under follower load in neutral standing posture,” Acta of Bioengineering and Biomechanics, vol. 12, no. 4, pp. 47–51, 2010.

K. Kim, Y. H. Kim, and S. Lee, “Investigation of optimal follower load path generated by trunk muscle coordination,” Journal of Biomechanics, vol. 44, no. 8, pp. 1614–1617, 2011.

I. A. F. Stokes and M. Gardner-Morse, “Quantitative anatomy of the lumbar musculature,” Journal of Biomechanics, vol. 32, no. 3, pp. 311–316, 1999.

M. M. Panjabi, R. A. Brand, and A. A. White, “Three dimensional flexibility and stiffness properties of the human thoracic spine,” Journal of Biomechanics, vol. 9, no. 4, pp. 185–192, 1976.

M. G. Gardner-Morse, J. P. Laible, and I. A. F. Stokes, “Incorporation of spinal flexibility measurements into finite element analysis,” Journal of Biomechanical Engineering, vol. 112, no. 4, pp. 481–483, 1990.

A. A. White III and M. M. Panjabi, Clinical Biomechanics of the Spine, Lippincott Williams & Wilkins, Philadelphia, Pa, USA, 1990.

A. B. Schultz, “Biomechanical analyses of loads on the lumbar spine,” in The Lumbar Spine, pp. 160–171, W. B. Saunders, Philadelphia, Pa, USA, 1990.

H. J. Wilke, P. Neef, M. Caimi, T. Hoogland, and L. E. Claes, “New in vivo measurements of pressures in the intervertebral disc in daily life,” Spine, vol. 24, no. 8, pp. 755–762, 1999.