# General Computational Model for Human Musculoskeletal System of Spine

Authors: Kyungsoo Kim; Yoon Hyuk Kim; SuKyoung Lee
Journal: Journal of Applied Mathematics
 (2012)

---

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


---

## Body

## 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 N vertebrae and M 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 ith vertebral body center is given as a node by pi, i=1,2,…,N. Let Fkm and PCSAk be the kth muscle force and the physiological cross-sectional area of the kth muscle for k=1,2,…,M. Assume that there are Mi muscles acting on ith vertebra among M trunk muscles and Fi,jm, j=1,2,…,Mi, denotes the jth muscle force vector starting from the attachment point in ith vertebra. Let pi,j be the position vector of the attachment point of jth muscle acting on ith 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(2.1)[F1ms⋮FNmsM1ms⋮MNms]=K⋅[d1d2⋮dN],
where Fims, Mims, and di, denote the motion segment force, the motion segment moment, and the displacement vectors at ith vertebral body center, respectively. K represents the stiffness matrix describing linear elasticity of the spine model. The stiffness matrix K of the motion segment can be obtained from experimental studies such as [23, 24]. The displacement vector di at ith node consists of the translation components, di,kt, k=1,2,…,K, and the rotation components, di,lr, l=1,2,…,L, where K and L are the number of translational and rotational degrees of freedom at each node, respectively.Then, for given external forcesFie and moments Mie applied at ith vertebral body center, the static equilibrium equations at the vertebral nodes can be formulated by(2.2)∑j=1MiFi,jm-Fims+Fie=0,i=1,2,…,N,(2.3)∑j=1Miri,j×Fi,jm-Mims+Mie=0,i=1,2,…,N,
where ri,j=pi-pi,j for all i and j, 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,Fijt, at ith vertebra is calculated iteratively: for i=1,2,…,N,(2.4)Fijt=∑j=1MiFi,jm+Fie+Fi-1jt=Fims+Fi-1jt
with F0jt=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 vectorci at ith node be(2.5)ci=pi-pi-1‖pi-pi-1‖,i=1,2,…,N,
under the assumption that p0=0, which indicates the direction of ith beam element.Then,Fijt can be decomposed into two perpendicular compressive forces Fic=(Fijt·ci)ci and shear force Fis=Fijt-(Fijt·ci)ci at ith node for i=1,2,…,N (Figure 1) as(2.6)Fijt=(Fijt⋅ci)ci+(Fijt-(Fijt⋅ci)ci)=Fic+Fis.
The resultant joint moment Mijt at ith node for i=1,2,…,N is the same to the motion segment moment Mims.Figure 1
Decomposition of joint forceFijt into the compressive force Fic and the shear force Fis.


## 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 N vertebrae and M 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 ith vertebral body center is given as a node by pi, i=1,2,…,N. Let Fkm and PCSAk be the kth muscle force and the physiological cross-sectional area of the kth muscle for k=1,2,…,M. Assume that there are Mi muscles acting on ith vertebra among M trunk muscles and Fi,jm, j=1,2,…,Mi, denotes the jth muscle force vector starting from the attachment point in ith vertebra. Let pi,j be the position vector of the attachment point of jth muscle acting on ith 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(2.1)[F1ms⋮FNmsM1ms⋮MNms]=K⋅[d1d2⋮dN],
where Fims, Mims, and di, denote the motion segment force, the motion segment moment, and the displacement vectors at ith vertebral body center, respectively. K represents the stiffness matrix describing linear elasticity of the spine model. The stiffness matrix K of the motion segment can be obtained from experimental studies such as [23, 24]. The displacement vector di at ith node consists of the translation components, di,kt, k=1,2,…,K, and the rotation components, di,lr, l=1,2,…,L, where K and L are the number of translational and rotational degrees of freedom at each node, respectively.Then, for given external forcesFie and moments Mie applied at ith vertebral body center, the static equilibrium equations at the vertebral nodes can be formulated by(2.2)∑j=1MiFi,jm-Fims+Fie=0,i=1,2,…,N,(2.3)∑j=1Miri,j×Fi,jm-Mims+Mie=0,i=1,2,…,N,
where ri,j=pi-pi,j for all i and j, 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,Fijt, at ith vertebra is calculated iteratively: for i=1,2,…,N,(2.4)Fijt=∑j=1MiFi,jm+Fie+Fi-1jt=Fims+Fi-1jt
with F0jt=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 vectorci at ith node be(2.5)ci=pi-pi-1‖pi-pi-1‖,i=1,2,…,N,
under the assumption that p0=0, which indicates the direction of ith beam element.Then,Fijt can be decomposed into two perpendicular compressive forces Fic=(Fijt·ci)ci and shear force Fis=Fijt-(Fijt·ci)ci at ith node for i=1,2,…,N (Figure 1) as(2.6)Fijt=(Fijt⋅ci)ci+(Fijt-(Fijt⋅ci)ci)=Fic+Fis.
The resultant joint moment Mijt at ith node for i=1,2,…,N is the same to the motion segment moment Mims.Figure 1
Decomposition of joint forceFijt into the compressive force Fic and the shear force Fis.


## 3. Formulation of Optimization Scheme

### 3.1. Assumptions for Physiology

In this study, the displacement vectordi at ith vertebral body center for i=1,2,…,N and the kth muscle force Fkm for k=1,2,…,M 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(3.1)f(d1,…,dN,F1m,…,FMm)=w1∑i=1N‖Fic‖2+w2∑i=1N‖Fis‖2+w3∑i=1N‖Mijt‖2+w4∑k=1M|Fkm|2+w5∑k=1M|FkmPCSAk|3,
where w1,w2,w3,w4, and w5 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(3.2)f(d1,…,dN,F1m,…,FMm)=w1∑i=1N‖Fic‖2+w2∑i=1N‖Fis‖2+w3∑i=1N‖Mijt‖2,(3.3)f(d1,…,dN,F1m,…,FMm)=∑k=1M|Fkm|2,
or(3.4)f(d1,…,dN,F1m,…,FMm)=∑k=1M|FkmPCSAk|3.
In addition, the follower load concept can be formulated by a constraint as(3.5)‖Fic‖≤α‖Fis‖,
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(3.6)f(d1,…,dN,F1m,…,FMm)s. t.(1)
Fm-K·d+Fe=0, where
(3.7)Fm=[∑j=1M1F1,jm⋮∑j=1MNFN,jm∑j=1M1r1,j×F1,jm⋮∑j=1MNrN,j×FN,jm],d=[d1⋮⋮⋮dN],Fe=[F1e⋮FNeM1e⋮MNe],
and K is the stiffness matrix defined in (2.1);(2)
∥Fic∥≤α∥Fis∥,i=1,2,…,N;(3)
|Fkm/PCSAk|≤σ, k=1,2,…,M where σ is the maximum muscle stress;(4)
0≤|di,kt|≤di,k,maxt and 0≤|di,lr|≤di,l,maxr, i=1,2,…,N where di,k,maxt and di,l,maxr denote the upper bounds of kth translation component and of lth rotation component of di.


## 3.1. Assumptions for Physiology

In this study, the displacement vectordi at ith vertebral body center for i=1,2,…,N and the kth muscle force Fkm for k=1,2,…,M 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(3.1)f(d1,…,dN,F1m,…,FMm)=w1∑i=1N‖Fic‖2+w2∑i=1N‖Fis‖2+w3∑i=1N‖Mijt‖2+w4∑k=1M|Fkm|2+w5∑k=1M|FkmPCSAk|3,
where w1,w2,w3,w4, and w5 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(3.2)f(d1,…,dN,F1m,…,FMm)=w1∑i=1N‖Fic‖2+w2∑i=1N‖Fis‖2+w3∑i=1N‖Mijt‖2,(3.3)f(d1,…,dN,F1m,…,FMm)=∑k=1M|Fkm|2,
or(3.4)f(d1,…,dN,F1m,…,FMm)=∑k=1M|FkmPCSAk|3.
In addition, the follower load concept can be formulated by a constraint as(3.5)‖Fic‖≤α‖Fis‖,
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(3.6)f(d1,…,dN,F1m,…,FMm)s. t.(1)
Fm-K·d+Fe=0, where
(3.7)Fm=[∑j=1M1F1,jm⋮∑j=1MNFN,jm∑j=1M1r1,j×F1,jm⋮∑j=1MNrN,j×FN,jm],d=[d1⋮⋮⋮dN],Fe=[F1e⋮FNeM1e⋮MNe],
and K is the stiffness matrix defined in (2.1);(2)
∥Fic∥≤α∥Fis∥,i=1,2,…,N;(3)
|Fkm/PCSAk|≤σ, k=1,2,…,M where σ is the maximum muscle stress;(4)
0≤|di,kt|≤di,k,maxt and 0≤|di,lr|≤di,l,maxr, i=1,2,…,N where di,k,maxt and di,l,maxr denote the upper bounds of kth translation component and of lth rotation component of di.


## 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 (N=7). Here, 117 pairs of trunk muscles were considered (M=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 K was obtained from previous experimental studies [23, 24].In this test, (3.2) was used for the objective function. The weight factors w1,w2, and w3 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 Table1. 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 forceNumber of activated muscles0%–20%720%–40%740%–60%660%–80%280%–100%20Total42


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



---

*Source: 484759-2012-03-01.xml*