Table of Contents Author Guidelines Submit a Manuscript
Journal of Applied Mathematics
Volume 2012 (2012), Article ID 484759, 8 pages
http://dx.doi.org/10.1155/2012/484759
Research Article

General Computational Model for Human Musculoskeletal System of Spine

1Department of Mathematics, Kyonggi University, Suwon 443-760, Republic of Korea
2Department of Mechanical Engineering, Kyung Hee University, Yongin 446-701, Republic of Korea
3e-Spine Center, Kyung Hee University, Yongin 446-701, Republic of Korea
4Department 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 [57], 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, 1121]. Although the follower load concept was considered in [14, 1621], 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𝑚𝑠𝐌𝑁𝑚𝑠𝐝=𝐊1𝐝2𝐝𝑁,(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 𝐩0=𝟎, 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 𝐌𝑖𝑚𝑠.

484759.fig.001
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𝑁𝑖=1𝐅𝑐𝑖2+𝑤2𝑁𝑖=1𝐅𝑠𝑖2+𝑤3𝑁𝑖=1𝐌𝑖𝑗𝑡2+𝑤4𝑀𝑘=1||𝐹𝑚𝑘||2+𝑤5𝑀𝑘=1||||𝐹𝑚𝑘𝑃𝐶𝑆𝐴𝑘||||3,(3.1) where 𝑤1,𝑤2,𝑤3,𝑤4, and 𝑤5 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𝑁𝑖=1𝐅c𝑖2+𝑤2𝑁𝑖=1𝐅𝑠𝑖2+𝑤3𝑁𝑖=1𝐌𝑖𝑗𝑡2,(3.2)𝑓𝐝1,,𝐝𝑁,𝐹𝑚1,,𝐹𝑚𝑀=𝑀𝑘=1||𝐹𝑚𝑘||2,(3.3) or𝑓𝐝1,,𝐝𝑁,𝐹𝑚1,,𝐹𝑚𝑀=𝑀𝑘=1||||𝐹𝑚𝑘𝑃𝐶𝑆𝐴𝑘||||3.(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,,𝐝𝑁,𝐹𝑚1,,𝐹𝑚𝑀(3.6)

s. t.(1)𝐅𝑚𝐊𝐝+𝐅𝑒=𝟎, where 𝐅𝑚=𝑀1𝑗=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, 1922]. 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 𝑤1,𝑤2, and 𝑤3 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.

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

References

  1. A. Nachemson, “Lumbar intradiscal pressure,” in The Lumbar Spine and Back Pain, M. I. V. Jayson, Ed., chapter 9, pp. 191–203, 1987. View at Google Scholar
  2. 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. View at Google Scholar
  3. 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. View at Publisher · View at Google Scholar · View at Scopus
  4. 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. View at Publisher · View at Google Scholar · View at Scopus
  5. 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. View at Publisher · View at Google Scholar · View at Scopus
  6. 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. View at Publisher · View at Google Scholar · View at Scopus
  7. 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. View at Publisher · View at Google Scholar · View at Scopus
  8. 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. View at Publisher · View at Google Scholar · View at Scopus
  9. 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. View at Publisher · View at Google Scholar · View at Scopus
  10. G. Herrmann, “Stability of equilibrium of elastic systems subjected to nonconservative forces,” Applied Mechanics Reviews, vol. 20, pp. 103–108, 1967. View at Google Scholar
  11. 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. View at Publisher · View at Google Scholar
  12. 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. View at Publisher · View at Google Scholar · View at Scopus
  13. 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. View at Publisher · View at Google Scholar · View at Scopus
  14. 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. View at Publisher · View at Google Scholar · View at Scopus
  15. 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. View at Publisher · View at Google Scholar · View at Scopus
  16. 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. View at Publisher · View at Google Scholar · View at Scopus
  17. 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. View at Publisher · View at Google Scholar · View at Scopus
  18. 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. View at Publisher · View at Google Scholar · View at Scopus
  19. 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. View at Publisher · View at Google Scholar · View at Scopus
  20. 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. View at Google Scholar
  21. 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. View at Publisher · View at Google Scholar
  22. 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. View at Publisher · View at Google Scholar · View at Scopus
  23. 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. View at Publisher · View at Google Scholar · View at Scopus
  24. 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. View at Publisher · View at Google Scholar · View at Scopus
  25. A. A. White III and M. M. Panjabi, Clinical Biomechanics of the Spine, Lippincott Williams & Wilkins, Philadelphia, Pa, USA, 1990.
  26. 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. View at Google Scholar
  27. 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. View at Publisher · View at Google Scholar · View at Scopus