Abstract

Despite development of accurate musculoskeletal models for human lumbar spine, the methods for prediction of muscle activity patterns in movements lack proper association with corresponding sensorimotor integrations. This paper uses the directional information of the Jacobian of the musculoskeletal system to orchestrate adaptive critic-based fuzzy neural controller modules for controlling a complex nonlinear redundant musculoskeletal system. The proposed controller is used to control a 3D 3-degree of freedom (DOF) musculoskeletal model of trunk, actuated by 18 muscles. The controller is capable of learning to control from sensory information, without relying on pre-assumed model parameters. Simulation results show satisfactory tracking of movements and the simulated muscle activation patterns conform to previous EMG experiments and optimization studies. The proposed controller can be used as a computationally inexpensive muscle activity generator to distinguish between neural and mechanical contributions to movement and for study of sensory versus motor origins of motor function and dysfunction in human spine.

Introduction

Understanding the biomechanics of the human lumbar spine is essential for assessment of risk of injury in healthy individuals and in low back pain patients [1]. Calculation of muscle activation patterns in different isometric tasks or 3D movements are especially of interest from a clinical perspective [2,3] and to study the motor disorders that cause biomechanical dysfunction [4].

Biomechanical approaches to resolve kinetic redundancy (mapping of few joint torques to many muscle activation levels) usually incorporate experimental electromyography (EMG) or use optimality assumptions that a cost function is being minimized during control [5–8]. Other computational methods to predict muscle activation patterns, e.g., neural networks-based methods [9,10] or computed muscle control [11], are commonly based on similar assumptions [1]. The use of artificial neural networks for generalization of experimental EMG [12,13], or for predictive models based on static equilibrium and within-network competition [14] by Nussbaum et al. are notable. The kinematic-based finite-element models [15] are accurate in terms of musculoskeletal biomechanics, but do not address the role of neural control mechanisms in trunk movements. Experimental and clinical data suggest that the complexity of muscle recruitment patterns cannot always be explained by simple optimality principles [8]. However, the methods that give better conformity with physiological data have been only used for simple models in upper limb and suffer from computational costs [16–20].

In this context, intelligent control methods can provide solutions for including the role of sensory-motor integration and learning in control [21]. An example is the application of dynamic neural networks for control of a simplified trunk model [22]. Other studies for incorporating the sensorimotor physiology into control are limited to static postures [23,24]. Notably, a gradient descent method with reference to muscle moment arms and joint reaction force for solving isometric exertions has been reported [23]; however, the method is computationally expensive and is limited to static postures.

In this study, we show how the sensory information can provide an estimate of the Jacobian of the musculoskeletal system which can then be used to orchestrate excitation-inhibition-based controller units. This array of single-input single-output (SISO) controller units makes a multi-input multi-output (MIMO) controller which is consequently used to control the musculoskeletal system by generating the required muscle activity levels for a given desired trajectory. In what follows, the main components and the architecture of the “controller” will be explained. Using a 3D musculoskeletal “model” of the human torso, the capability of the controller to control complex nonlinear systems, robustness and the ability to learn to control from sensory information is tested. More importantly, the conformity of the generated muscle activity levels to previous EMG [25,26] or state of the art optimization-based studies [6] are evaluated.

Methods

Musculoskeletal Model.

The skeletal and muscular parts of the model are primarily chosen to be simple and easy to interpret, yet not very distant from the actual anatomical and physiological properties [6].

Skeletal System.

The whole trunk is considered as a 3D rigid body that allows realistic modeling of origin and insertion points of muscles and generation of moment about all 3 axes by muscles. It is supposed to move around a frictionless ball and socket joint at the level of the lumbosacral (L5/S1) joint. The kinematics and dynamics of this 3D rigid body is defined by a 3D inverted pendulum with three rotational DOF. The orientation of inverted pendulum is described by three Euler angles which are obtained by rotation about x, y and z axes and describe lateral flexion, axial rotation and flexion/extension, respectively [27]. The governing equation of the constrained inverted pendulum is written as
$I~ω→·+ω→×(I~ω→)=τ→+a→×(f→+g→)Θ→·=B~(Θ→)ω→$
(1)

where $I~$, $ω→$, $τ→$, $a→$, $f→$, and $g→$ are the moment of inertia matrix about the joint axes, angular velocity vector, external applied torque vector, position vector of the center of gravity of body, external applied force vector, and gravity force vector, all expressed in a body coordinate system that rotates with the trunk. $Θ→=[θ1 θ2 θ3]T$ is the vector of Euler angles and $B~$ the transformation matrix between angular velocity vector and the rate of Euler angles vector [27]. All the model parameters and their values are listed in the Appendix, Table A1 (See the “Supplemental Data” tab for this paper on the ASME Digital Collection.)

Musculature.

Five back and abdominal muscle groups were considered to contribute to 3D movements and resisting external forces [28]. These groups of muscles, including rectus abdominis (RA), erector spinae (ES), internal oblique (IO), external oblique (EO), and latissimus dorsi (LD), were modeled using 18 vectors with Hill-type muscle model [29]. Additionally, the passive force of three muscle groups, including pars lumborum, psoas, and quadratus lumborum were included in the model. In order to model the muscles, anatomical data including muscle insertions and origins, physiological cross section area (PCSA) and rest length of muscles were obtained from literature [30]. Figure 1 illustrates the spatial distribution and alignment of the five actively modeled muscle groups, using labeled anatomical drawings. The passive force [31], force-length [32], and force-velocity [33] relationships were approximated, as presented by the following equations:

Fig. 1
Fig. 1
Close modal
$fp(l)=exp(-10.671+7.675(l/l0))$
(2)
$fl(l)=5.144-29.14(l/l0)+55.92(l/l0)2-40.80(l/l0)3+9.868(l/l0)4$
(3)
$fv(l·)=0.14330.1074+exp(-1.409sinh(3.2(l·/l·max)+1.6))$
(4)
where $fp(l)$, $fl(l)$, and $fv(l·)$ are passive force, force-length, and force-velocity functions and whole muscle force is given by
$f=(α·fl(l)·fv(l·)+fp)fmax while fmax=PCSA·σ$
(5)

And $α$, $l$, $l·$, $l0$, and $l·max$ are muscle activation level (bounded between 0 and 1), muscle length, rate of change of muscle length, resting length and maximum velocity of muscle, respectively. A recommended value $l·max=l0/0.1s$ was used [29]. The value of $fmax$ is the maximum force of muscle which is the product of maximum muscle stress $σ$ and PCSA of each muscle. Muscle stress was taken to be 55 $N/cm2$ [6] as previously described for general muscle models [29,34], and for trunk muscles [15]. Inspection of muscles' length in different postures according to muscle resting lengths [30] indicated that in upright standing posture most muscles are nearly in their resting lengths, while they approach their minimum and maximum lengths as system takes postures in extreme ranges of motion [35].

Controller.

The problem of tracking control of the described musculoskeletal model is a MIMO control problem including 3 controlled variables and 18 control signals. The constituent elements of the controller are explained as follows.

Critic-Based Fuzzy Neurocontroller Units.

The SISO controller “units” in this study are adaptive fuzzy neurocontrollers whose weights are instantaneously “tuned” by an analogue critic signal [36,37]. These excito-inhibitory units use an adaptive gradient descent algorithm to continuously tune controller weights to penalize it for the critic signals (see Fig. 2). The structure and tuning rules of each controller unit is explained in the Appendix A1 (See the “Supplemental Data” tab for this paper on the ASME Digital Collection.). The total error function to be minimized is defined as:

Fig. 2
Fig. 2
Close modal
$E=Ee+Eα=12ke(h1e+h2e·)2+12kα(abs(α))2$
(6)
where $e$ and $e·$ are error and rate of error of tracking (the difference between the model's actual kinematics and the desired kinematics), $h1$, and $h2$ scale factors to bring the products in [−1,1] range, and $ke$, and $kα$ are the weights that determine the priorities of components in error signal formation. $α$ is the activation level of controlled actuator (muscle). The error function consists of $Ee$ which is an indicator of control error, and $Eα$ which penalizes the controller unit for control effort signal. This $Eα$ component plays an essential role in control of redundant systems. By taking into account that scalar $j=∂θ/∂α$ represents the Jacobian of the SISO controlled system and the control effort is a positive value (activation levels of muscles are positive values), the overall weight tuning is computed by the following equation:
$Δwi=η(ke·h1·re·j-kα·rα)·∂α∂wi$
(7)

where $wi$ is the weight for the ith neuron of the neurocontroller (See Appendix A1 in the “Supplemental Data” tab for this paper on the ASME Digital Collection), $η$ is the tuning rate which determines the rate of change of weights and $re=h1e+h2e·$, and $rα=abs(α)$.

It should be noted that for SISO control, a very rough estimate of j is usually adequate. In practice, even the sign of j works. To consider the bounded control signals, conditional weight tuning is used to ensure that the controller output is always in the feasible range. The derived tuning rules are used for online tuning of controller weights.

Directional Encoding (DE).

To perform the control task, parallel operation of several controller units, in a multi-agent framework is required [37]. This coordination for the role of each controller unit in the array is managed by extended formulation and use of the elements of the Jacobian matrix in tuning rules, which is referred to as directional encoding framework. This leads to 18 separate neurocontroller units (one for each muscle) that operate in parallel (Fig. 3), and are tuned simultaneously. The elements of the Jacobian matrix, defined in Eq. (8), are devised to generalize the controllers learning rules

Fig. 3
Fig. 3
Close modal
$J~=∂Θ→∂α→$
(8)
where $α→$ is the vector of activation levels of all muscles. The generalized control error for the MIMO controller is written as
$Ee=Ee1+Ee2+Ee3=12ke·(re12+re22+re32)=12ke·((h1e1+h2e·1)2+(h1e2+h2e·2) 2+(h1e3+h2e·3) 2)$
(9)
where $Ee1$, $Ee2$, and $Ee3$ are the error functions that are used to penalize the controllers for error in first, second, and third Euler angles, respectively. The coefficients for $Eei$ (error functions corresponding to $θi$) are chosen to be equal for simplicity. As a result, the generalized tuning rule for the controller weights is:
$Δwmi=η(ke·h1·(Jm1·re1+Jm2·re2+Jm3·re3)-kα·rα)·∂αm∂wmi$
(10)
where indices $m$ and $i$ are related to mth unit and ith neuron. To minimize the interaction of controller units on another, we may use cosine tunings. This means that contribution of each SISO controller unit to error compensation is weighted as a function of instantaneous Jacobian of controlled system: when the control output of unit (and hence the generated force for corresponding muscle) is in the direction of instantaneous error the controller unit will fire maximally and by deviation from the error direction the weight of the controller is decreased by a cosine function. This line of action direction of muscle is a function of instantaneous geometrical state or configuration of musculoskeletal system and the directionally encoded error will be the error along each muscle's line of action. From mathematical point of view, Eq. (10) can be modified in the following form:
$Δwmi=η(ke·h1·(jm·rem)-kα·rα)·∂αm∂wmi$
(11)

where $rem$ and $jm$ are the directionally encoded error reward signal and scalar Jacobian term for mth controller. It can be seen by inspection that $rem$ can be computed by dot-product of unit vector of mth column of the Jacobian matrix and the error or reward vector.

Equation (11) is the weight tuning rule for mth controller unit with directional encoding and can be used with the appropriate architecture shown in Fig. 3.

In order to geometrically interpret the new mathematical formulation, the elements of the Jacobian matrix can be written as a function of model parameters. As mentioned before, a rough estimate of Jacobian matrix is adequate for the proposed control application. By taking into account that the directions of applied moments and angles of rotation are the same due to the chosen convention, Eqs. (3) and (4) describe positive functions and maximum force in Eq. (5) equals $(PCSA·σ)$, an estimate of the Jacobian matrix may be used
$Jmj∝dmj·PCSAm$
(12)

where $∝$ is the proportionality symbol, $dmj$ is the corresponding member of moment arm matrix of muscles about center of rotation of the pendulum (representing the moment arm of mth muscle about the jth direction), and $PCSAm$ is the physiological cross-sectional area of the corresponding muscle, respectively. As Eq. (12) shows, the $Jmj$ is approximated by $(dmj·PCSAm)$ or simply by $dmj$. That is, the error vector $e→$ is in fact projected along the moment arm vector, assigning the appropriate components of the error vector to the relevant actuators.

The resultant synergy, implicitly generated by this formulation, is a very close match to the neural synergy discussed by Neilson and Neilson [38,39] where the contribution weight of a muscle to the generated joint moment is determined by the sensory estimate of the moment arm. Here the same sensory estimate is used to determine the rate of change of the contribution weight of the muscle activity. According to the components of the Jacobian matrix, this synergy can be brought to the central nervous system (CNS) from sensory information.

Supervised Learning.

The sensory information from muscle spindles, as the most important mechanoreceptors, reflects the muscle length and rate of change of length [40]. The state of limb is found by state estimator network, whose design details are out of scope of this paper. As a result, the muscle length and joint state matrices are provided for the controller. On the other hand, the moment arm matrix in Eq. (12) can also be obtained by Eq. (13), based on principle of virtual work [41]
$dmj=∂lm∂θj$
(13)

Equation (13) shows that a simple differentiation is needed to compute $dmj$ from $lm$ and $θj$. Consequently the moment arm matrix $dmj$, as the main constituent of $J~$ matrix, is available to the system, and can be learned by a neural network through a supervised “learning” algorithm. In this work, an array of static feed-forward neural networks is used to learn the required mapping. Each network consists of a 3-node sigmoid input layer for three Euler angles, 21-node sigmoid hidden layer and a 1-node linear output layer, trained via offline Levenberg–Marquat algorithm for faster training and avoidance of relevant training problems [42].

Controller Architecture.

Figure 3 illustrates the multi-agent neurofuzzy architecture, based on directional encoding of error signals. The Jacobian matrix introduces the directional information of actuators (muscles) to the controller. Furthermore, the scalar $jm$ obtained from the Jacobian matrix affects the cooperation between agents. We may succinctly consider the role of each neurofuzzy controller (NFC) block as excitation-inhibition according to this directional and weight information.

The role of the proposed control architecture in the overall motor control hierarchy is shown in Fig. 4. The arrangement of modules is considerably simplified for this study compared to previous literature [16,17].

Fig. 4
Fig. 4
Close modal

It is also noteworthy that the learning process for the Jacobian network is a rather long-term learning, and is separate from the online fast weight tuning in neurocontroller units.

Simulation Cases.

In order to evaluate the performance of the proposed control architecture, study the effect of learning on tracking performance, and study the simulated motor patterns and corresponding characteristics in movements, several simulations were performed in various planes. The movements were chosen to be in normal range of motion of lumbar spine and cover the most common types of movement in this region [35]. Flexion, axial rotation, lateral bending, as well as combinations of these movements were simulated in the form of point to point movements or oscillatory maneuvers. The movements were simulated from neutral position to arbitrary positions. The minimum-jerk trajectory was used for discrete movement trajectories [16], hence the job of controller is to solve for muscle activation levels for given desired kinematics. The time step for all simulation cases were set to 0.01 s. The tuning rate $η$ for Critic-based fuzzy neurocontroller agents were chosen to be 50. The initial tuning phase takes place in the beginning of movement. After initial formation of weights, the tuning continues and controller agents can accomplish the control task properly. Simulations use an already trained Jacobian neural network, except the simulation in Sec. 3.2, where the effect of the Jacobian neural network training is studied.

Controller stability is assessed by applying perturbations to the biomechanical system: the abrupt increase of system mass by 50% at 0.15 s of the movement and return to its initial value at 0.75 s of movement; and application of −100 Nm moment during the [0.45 s, 0.52 s] time interval.

Results

Controller Performance.

To analyze the controller performance several simulations were performed after complete offline training of the Jacobian neural network. The initial online tuning of fuzzy neurocontrollers took only about 2 s (200 iterations). The presented results show the performance after initial tuning phase. A summary of simulation results within the spine range of motion are listed in Table 1. The mean, maximum and steady state tracking errors were below 1.11 deg, 1.87 deg, and 2.46 deg, respectively. Figure 5 shows the reference angle, actual angle and tracking error during flexion movement from upright 0 deg to 55.4 deg flexion. The low tracking errors in the flexion movement was slightly affected by movement duration: in a faster movement with 0.8 s duration RMS tracking error was increased by 16%, while in a slower 1.3 s movement, the tracking error was reduced by 13%. In both cases, steady state error was only affected by less than 0.5%. Tracking a nonsinusoidal reference oscillation (tested by adding a 3 Hz sinusoid with 8% of the 1 Hz harmonic amplitude) increased the RMS tracking error by 26%.

Fig. 5
Fig. 5
Close modal
Table 1

Summary of controller's tracking performance in various movement simulations. Duration of movement in point to point movements and half period of oscillatory movements is 1.0 s in all cases. Simulations are performed for movements within normal range of motion of lumbar spine, where $θ1$ = lateral bending, $θ2$ = axial rotation, and $θ3$ = extension/flexion.

CasePlanes of movementMovement typeInitial state (deg)Final state (deg)RMS tracking error deg (%)Maximum error deg (%)Steady state error deg (%)
1PlanarOscillatory$θ1=0θ2=0θ3=0$$θ1=0θ2=0θ3=-55.4$0.93 (1.68)1.49 (2.70)N/A
2PlanarPoint-to-point$θ1=0θ2=0θ3=0$$θ1=0θ2=0θ3=-55.4$0.61 (1.10)1.21 (2.18)0.58 (1.06)
3PlanarOscillatory$θ1=0θ2=0θ3=0$$θ1=-22.8θ2=0θ3=0$0.33 (1.43)0.52 (2.30)N/A
4PlanarPoint-to-point$θ1=0θ2=0θ3=0$$θ1=-22.8θ2=0θ3=0$0.23 (1.03)0.35 (1.52)0.49 (2.16)
5PlanarOscillatory$θ1=0θ2=0θ3=0$$θ1=0θ2=-14.1θ3=0$0.35 (2.46)0.49 (3.50)N/A
6PlanarPoint-to-point$θ1=0θ2=0θ3=0$$θ1=0θ2=-14.1θ3=0$0.37 (2.66)0.63 (4.51)0.82 (5.85)
7Three-dimensionalOscillatory$θ1=0θ2=0θ3=0$$θ1=-22.8θ2=-14.1θ3=0$0.48 (1.77)0.72 (2.68)N/A
8Three-dimensionalPoint-to-point$θ1=0θ2=0θ3=0$$θ1=-22.8θ2=-14.1θ3=0$0.57 (2.14)0.93 (3.49)1.09 (4.06)
9Three-dimensionalOscillatory$θ1=0θ2=0θ3=0$$θ1=0θ2=-14.1θ3=-55.4$1.17 (2.06)1.77 (3.10)N/A
10Three-dimensionalPoint-to-point$θ1=0θ2=0θ3=0$$θ1=0θ2=-14.1θ3=-55.4$1.04 (1.83)1.52 (2.66)1.65 (2.94)
11Three-dimensionalOscillatory$θ1=0θ2=0θ3=0$$θ1=-22.8θ2=0θ3=-45$1.01 (1.69)1.68 (2.81)N/A
12Three-dimensionalPoint-to-point$θ1=0θ2=0θ3=0$$θ1=-22.8θ2=0θ3=-45$1.11 (1.85)1.87 (3.13)2.46 (4.10)
13Three-dimensionalOscillatory$θ1=0θ2=0θ3=0$$θ1=-14.1θ2=-14.1θ3=-14.1$0.73 (3.00)1.20 (4.91)N/A
14Three-dimensionalPoint-to-point$θ1=0θ2=0θ3=0$$θ1=-14.1θ2=-14.1θ3=-14.1$0.83 (3.38)1.43 (5.86)1.93 (7.91)
CasePlanes of movementMovement typeInitial state (deg)Final state (deg)RMS tracking error deg (%)Maximum error deg (%)Steady state error deg (%)
1PlanarOscillatory$θ1=0θ2=0θ3=0$$θ1=0θ2=0θ3=-55.4$0.93 (1.68)1.49 (2.70)N/A
2PlanarPoint-to-point$θ1=0θ2=0θ3=0$$θ1=0θ2=0θ3=-55.4$0.61 (1.10)1.21 (2.18)0.58 (1.06)
3PlanarOscillatory$θ1=0θ2=0θ3=0$$θ1=-22.8θ2=0θ3=0$0.33 (1.43)0.52 (2.30)N/A
4PlanarPoint-to-point$θ1=0θ2=0θ3=0$$θ1=-22.8θ2=0θ3=0$0.23 (1.03)0.35 (1.52)0.49 (2.16)
5PlanarOscillatory$θ1=0θ2=0θ3=0$$θ1=0θ2=-14.1θ3=0$0.35 (2.46)0.49 (3.50)N/A
6PlanarPoint-to-point$θ1=0θ2=0θ3=0$$θ1=0θ2=-14.1θ3=0$0.37 (2.66)0.63 (4.51)0.82 (5.85)
7Three-dimensionalOscillatory$θ1=0θ2=0θ3=0$$θ1=-22.8θ2=-14.1θ3=0$0.48 (1.77)0.72 (2.68)N/A
8Three-dimensionalPoint-to-point$θ1=0θ2=0θ3=0$$θ1=-22.8θ2=-14.1θ3=0$0.57 (2.14)0.93 (3.49)1.09 (4.06)
9Three-dimensionalOscillatory$θ1=0θ2=0θ3=0$$θ1=0θ2=-14.1θ3=-55.4$1.17 (2.06)1.77 (3.10)N/A
10Three-dimensionalPoint-to-point$θ1=0θ2=0θ3=0$$θ1=0θ2=-14.1θ3=-55.4$1.04 (1.83)1.52 (2.66)1.65 (2.94)
11Three-dimensionalOscillatory$θ1=0θ2=0θ3=0$$θ1=-22.8θ2=0θ3=-45$1.01 (1.69)1.68 (2.81)N/A
12Three-dimensionalPoint-to-point$θ1=0θ2=0θ3=0$$θ1=-22.8θ2=0θ3=-45$1.11 (1.85)1.87 (3.13)2.46 (4.10)
13Three-dimensionalOscillatory$θ1=0θ2=0θ3=0$$θ1=-14.1θ2=-14.1θ3=-14.1$0.73 (3.00)1.20 (4.91)N/A
14Three-dimensionalPoint-to-point$θ1=0θ2=0θ3=0$$θ1=-14.1θ2=-14.1θ3=-14.1$0.83 (3.38)1.43 (5.86)1.93 (7.91)

Predicted Muscle Activity Patterns.

Figure 6 shows the generated activation levels for muscles, during the point-to-point movement of Fig. 5. The effect of movement duration on the predicted muscle activity patterns is shown in Fig. 7. It is shown that while all the activity bursts are scaled with movement duration, the second agonist burst tends to vanish in slower movements. Similar muscle activation pattern for the corresponding extension movement is shown in Fig. 8(a).

Fig. 6
Fig. 6
Close modal
Fig. 7
Fig. 7
Close modal
Fig. 8
Fig. 8
Close modal

Effect of Learning on Control.

The Effect of the Jacobian learning-error, on the control performance is shown in Fig. 9. It is shown that as the Jacobian information is consolidated, the tracking control of MIMO controller improves.

Fig. 9
Fig. 9
Close modal

Controller Robustness.

The RMS and maximum tracking error, due to presence of two perturbations increased to 0.75 deg and 1.54 deg (23% and 27% increase with respect to unperturbed condition, respectively).

Discussion

Controller Performance.

The proposed control strategy was applied to control a nonlinear, redundant unknown MIMO system with only sensory information available to the controller for learning. We briefly recap the operation of the controller, as depicted in Fig. 4. The proposed controller receives the desired trajectory (kinematics) as the input to the controller and sends the motor commands to the musculoskeletal model, accordingly. To properly do this job, the controller gains only partial knowledge of the model's dynamics (the Jacobian) in the training stage from the sensory information received from the model. After the learning is complete, the only knowledge the controller receives from the controlled system during control is its current states. It should be noted that the training and tuning is not specific to the input kinematics and the system could track many different kinematics with the same training (the system's Jacobian). As the results in Fig. 5 and Table 1 show, the controller was capable of satisfactory tracking in mono- and multiplanar movements after learning and initial tuning.

The performance beyond normal range of movement of lumbar spine degrades. This degradation (especially in axial rotation), can first be attributed to extremely nonlinear solution in activation-states space. This plant nonlinearity is mainly caused by the spatial orientation of muscle passive force. Second, none of the muscles are considered as the pre-emptive source of axial moment generation and the simple cosine tuning rule gives only little specificity to the activated muscles for each direction. That is, the accompanying torque in other directions, generated per unit of axial torque is considerable. This issue is not significant for moments generated in other directions.

Predicted Muscle Activities.

The controller uses an estimation of the Jacobian of the system from sensory information along with excito-inhibitory controller units to generate muscle activities. In fact, the successful tracking of the system by the control architecture finds a solution to the kinetically redundant inverse dynamics problem. While the subspace of the solutions in the muscle activity space can be generally solved with different additional constraints [5], the controller gives a solution that dynamically minimizes tracking error with minimum co-activation and with a distribution of muscle activity that changes in proportion to muscles' sensory involvement (roughly proportional to the moment arms). However, this emerged synergy is not hardwired and is determined by the requested task kinematics (motion profile, movement direction, amplitude and duration) through the control error. The simulations in this study were limited to deterministic conditions (assuming a fully trained Jacobian network), yielding similar predictions in different simulation runs of exactly the same movements. By extending the model to a stochastic one, variability or signal noise in the sensory information as well as in the actual and desired movement position and velocity can be simulated. This shall lead to muscle activations that all generate the same motion but with a slight trial-to-trial variability of the predicted EMG.

Due to large variability in the normal recruitment patterns of trunk muscles, an exact validation of the simulation results is not possible [1]. However, the temporal patterns and the magnitude of the simulated activation signals are similar to optimization-based predictions (Figs. 10(a) and 10(b) [6] and to experimental EMG profiles (Figs. 8 and 10(a)10(c) [25,26,43]. This similarity is promising and informs of the potentials of the proposed simulation method. More anatomically accurate models of the musculature can also improve the predicted muscle activities. As an example, modeling a multifunctional flat abdominal muscle such as Internal Oblique muscle by only two vectors (Fig. 1, based on the data from Ref. [30]) can be improved. This muscle should primarily contribute to trunk flexion due to its overall positive flexion moment arm. However, the slightly negative flexion moment arm of one of the modeled vectors (see Fig. 1, right), leads to the prediction of its activity during the deceleration phase of flexion (Figs. 6 and 8).

Fig. 10
Fig. 10
Close modal

The three-burst activation pattern (agonist-antagonist-agonist)2 [44,45] can be distinguished in the results. Additionally, the second predicted agonist burst tends to vanish in the slower movements (Fig. 7). This three-burst pattern is the common activation pattern in fast reaching movements [44,45], but can be seen only in some specific faster trunk flexion movements [43]. The predicted activation patterns, especially in faster movements, have good conformity with flexion EMG patterns of Oddsson and Thorstensson [43] (reproduced in Fig. 10(c)). For further validation, RMS filtered EMG activity of main trunk muscles (Rectus Abdominis and Erector Spinae), during extension movement, are shown in Fig. 8(b) (Data from Ref. [25] and reported in Ref. [26]). The sequence of activation, the magnitude and shape of muscle activities as well as their approximate maximum values, are in good conformity. Minor differences (such as differences in tonic activation levels at the start and end of movement) are due to biomechanical simplifications in the model, such as lack of multiple DOF afforded to each functional motion segment in the lumbar spine and not including the passive resistance of discs and ligaments.

We may also compare the simulation results with the activation levels generated by inverse dynamics and static optimization methods, which are the most widely used computational techniques to estimate the muscles activation patterns [1]. A comparison with optimization (with sum of squared activation as cost function) of Ref. [6] is presented in Fig. 10. The optimization results are computed using the same musculoskeletal dynamic model with all the main considered muscle groups present, as well as some extra muscle vectors included in modeling. Although the musculature is slightly different in the optimization study, there is very good agreement in the phase of activity of agonists and antagonists. The third burst of activity (AG2) in the generated activation by directional encoding based neurofuzzy controller (DE-NFC), cannot be predicted by optimization or stability based optimization. The oscillatory pattern of activity, including the third burst, observed at the latest movement period and after the movement is finished is a mechanism used by the CNS to accomplish fast and stable movements. Further discussion on this physiologically observed phenomenon [46] can be found elsewhere [16].

The focus on the simulation of muscle activity patterns in a complex musculoskeletal system such as trunk, can be considered a step forward compared to the previous applications of control theory and intelligent control for the study of human neuromotor control in simpler biomechanical systems [47–50].

Joint Forces.

We also calculated the compressive joint reaction force at L5/S1 level in flexion movement and compared it with the result of McGill et al. [51]. The model predicts the minimum force of 667 N in upright position and peak force of 3896 N at 0.84 s, while McGill et al. [51] estimate them to be 1068 N (46% difference) and 3461 N (12% difference). In comparison to the reported values of Bazrgari et al. [52], 500 N (approximated) and 2974 N, there is 28% and 27% differences, respectively. Also for static postures in flexion (60 deg), our model gives the compressive force of 2920 N compared to that of Arjmand et al. [53], 1939 N (40% difference). In all cases, the overall trends of variation of joint reaction force matches the cited literature.

Flexion relaxation phenomenon as an important characteristic in spine biomechanics expected at deep flexion [54] can be partly spotted in the results, which confirms the validity of the biomechanics.

Limitations.

The main purpose of this study was to include the potential sensorimotor integration mechanisms in modeling of spine movements and motor command generation. Future works may contain a complete hierarchy of motor control blocks (Fig. 4). Incorporation of co-contraction modulation into the model, and more elaborate Jacobian estimation from sensory information are other suggested improvements. More realistic biomechanical modeling, including musculotendon model [29,34] and passive tissues, and addition of reflexes and corresponding delays to the system shall provide a more realistic behavior of the overall neuromusculoskeletal system. The current model does not account for electromechanical delays. Consequently, in interpretation of the timing of the predicted muscle activities with respect to the movement kinematics, a time-shift of tens of milliseconds should be included for correction. This delay can be best modeled when other control and reflex delays in the system are also modeled.

Significance and Applications.

The main difference of our approach compared to optimization-based models is that we use no explicit optimization assumption or procedure. The muscle activity is predicted based on the error minimization guided by sensory information, and the penalty on α only prohibits co-activation but does not necessarily lead to an optimal pattern. Hence, the emerging solution is suboptimal or near-optimal. This pattern is generated by the controller units and by the guiding sensory information that contribute to the temporal and spatial characteristics of the activation patterns, respectively. Importantly, there is no systematic way in optimization-based methods to account for temporal and spatial features of the predicted muscle activity as a function of sensory information. The simulation results also show that while predictions are close to optimal, they can be notably different from the globally optimal solutions in some aspects (i.e., three-burst patterns rather than two-burst patterns predicted by optimization).

The comparison of the joint reaction forces and muscle activity patterns show agreement with previous modeling and experimental studies for flexion and extension movements (Figs. 8 and 10), i.e., two important and clinically relevant movements [3] among the repertoire of trunk movements. However, we want to emphasize that the main purpose of this study was to develop a method and a model that can link the potential sensorimotor integration mechanisms to musculoskeletal biomechanics, and not an accurate biomechanical modeling per se. Yet, the method does provide advantages over classic optimization techniques and in comparisons to experimental data (Figs. 8 and 10): when there is a mismatch between our method and optimization (second agonist burst in the three-burst pattern), it is the suggested model that gives a better conformity to experimental data. However, our method eventually has limitations regarding physiological predictions, as optimization and all other methods do (see Sec. 4.4 on limitations).

The main advantage of the proposed method is the ability to account for sensory versus motor contributions to muscle activity. Specifically, it is possible to model a motor deficiency, i.e., reduction of force generating capacity of a muscle without affecting the sensory role of the muscle and the controller's operation. Similarly, it is possible to simulate a sensor deficiency, i.e., reduction in the reported length change of a muscle, without altering the biomechanics of the system. Our preliminary simulations in isometric tasks (unpublished) show that these two types of deficits lead to abnormal but considerably distinct activation patterns. This can potentially lead to a useful investigation tool for the study of different pathological patterns of muscle activity in low back pain patients.

Conclusions

The proposed control method can be a promising candidate as a low-level motor command generator, as it has physiologically valid features of movement control, i.e., the three-burst EMG pattern and integration of sensory information. The notable ability of the system to learn to control or adapt according to sensory and motor characteristics of the system, makes it an attractive tool for study of sensory versus motor origins of motor function and dysfunction.

Acknowledgment

We thank all the members of the Biomechanics Laboratory at Sharif University of Technology for their invaluable suggestions, as well as the anonymous reviewers for the useful comments on this manuscript.

2

Definition of agonist/antagonist is based on the acceleration phase of movement.

References

1.
Reeves
,
N. P.
, and
Cholewicki
,
J.
,
2003
, “
Modeling the Human Lumbar Spine for Assessing Spinal Loads, Stability, and Risk of Injury
,”
Crit. Rev. Biomed. Eng
,
31
(
1–2
), pp.
73
139
.10.1615/CritRevBiomedEng.v31.i12.30
2.
Mousavi
,
S. J.
,
Olyaei
,
G. R.
,
Talebian
,
S.
,
Sanjari
,
M. A.
, and
Parnianpour
,
M.
,
2009
, “
The Effect of Angle and Level of Exertion on Trunk Neuromuscular Performance During Multidirectional Isometric Activities
,”
Spine
,
34
(
5
), pp.
E170
E177
.10.1097/BRS.0b013e31818aec05
3.
Larivière
,
C.
,
Gagnon
,
D.
, and
Loisel
,
P.
,
2000
, “
The Comparison of Trunk Muscles EMG Activation Between Subjects With and Without Chronic Low Back Pain During Flexion–Extension and Lateral Bending Tasks
,”
J. Electromyogr. Kinesiol.
,
10
(
2
), pp.
79
91
.10.1016/S1050-6411(99)00027-9
4.
Ebenbichler
,
G. R.
,
Oddsson
,
L. I.
,
Kollmitzer
,
J.
, and
Erim
,
Z.
,
2001
, “
Sensory-Motor Control of the Lower Back: Implications for Rehabilitation
,”
Med. Sci. Sports Exercise
,
33
(
11
), pp.
1889
1898
.10.1097/00005768-200111000-00014
5.
Rashedi
,
E.
,
Khalaf
,
K.
,
Nassajian
,
M. R.
,
Nasseroleslami
,
B.
, and
Parnianpour
,
M.
,
2010
, “
How Does the Central Nervous System Address the Kinetic Redundancy in the Lumbar Spine? Three-Dimensional Isometric Exertions With 18 Hill-Model-Based Muscle Fascicles at the L4–L5 Level
,”
Proc. Inst. Mech. Eng., Part H
,
224
(
3
), pp.
487
501
.10.1243/09544119JEIM668
6.
Zeinali-Davarani
,
S.
,
Hemami
,
H.
,
Barin
,
K.
,
,
A.
, and
Parnianpour
,
M.
,
2008
, “
Dynamic Stability of Spine Using Stability-Based Optimization and Muscle Spindle Reflex
,”
IEEE Trans. Neural Syst. Rehabil. Eng.
,
16
(
1
), pp.
106
118
.10.1109/TNSRE.2007.906963
7.
Li
,
G.
,
Kaufman
,
K. R.
,
Chao
,
E. Y.
, and
Rubash
,
H. E.
,
1999
, “
Prediction of Antagonistic Muscle Forces Using Inverse Dynamic Optimization During Flexion/Extension of the Knee
,”
J. Biomech. Eng.
,
121
(
3
), pp.
316
322
.10.1115/1.2798327
8.
Amarantini
,
D.
,
Rao
,
G.
, and
Berton
,
E.
,
2010
, “
A Two-Step EMG-and-Optimization Process to Estimate Muscle Force During Dynamic Movement
,”
J. Biomech.
,
43
(
9
), pp.
1827
1830
.10.1016/j.jbiomech.2010.02.025
9.
Hou
,
Y.
,
,
J. M.
,
Karwowski
,
W.
,
Marras
,
W. S.
, and
Davis
,
K.
,
2007
, “
Estimation of the Dynamic Spinal Forces Using a Recurrent Fuzzy Neural Network
,”
IEEE Trans. Syst. Man Cybern., Part B: Cybern.
,
37
(
1
), pp.
100
109
.10.1109/TSMCB.2006.881298
10.
Lee
,
W.
,
Karwowski
,
W.
,
Marras
,
W. S.
, and
Rodrick
,
D.
,
2003
, “
A Neuro-Fuzzy Model for Estimating Electromyographical Activity of Trunk Muscles Due to Manual Lifting
,”
Ergonomics
,
46
(
1–3
), pp.
285
309
.10.1080/00140130303520
11.
Thelen
,
D. G.
,
Anderson
,
F. C.
, and
Delp
,
S. L.
,
2003
, “
Generating Dynamic Simulations of Movement Using Computed Muscle Control
,”
J. Biomech.
,
36
(
3
), pp.
321
328
.10.1016/S0021-9290(02)00432-3
12.
Nussbaum
,
M. A.
,
Chaffin
,
D. B.
, and
Martin
,
B. J.
,
1995
, “
A Back-Propagation Neural Network Model of Lumbar Muscle Recruitment During Moderate Static Exertions
,”
J. Biomech.
,
28
(
9
), pp.
1015
1024
.10.1016/0021-9290(94)00171-Y
13.
Nussbaum
,
M. A.
, and
Chaffin
,
N. B.
,
1996
, “
Evaluation of Artificial Neural Network Modelling to Predict Torso Muscle Activity
,”
Ergonomics
,
39
(
12
), pp.
1430
1444
.10.1080/00140139608964562
14.
Nussbaum
,
M. A.
,
Martin
,
B. J.
, and
Chaffin
,
D. B.
,
1997
, “
A Neural Network Model for Simulation of Torso Muscle Coordination
,”
J. Biomech.
,
30
(
3
), pp.
251
258
.10.1016/S0021-9290(96)00138-8
15.
Arjmand
,
N.
, and
,
A.
,
2006
, “
Model and in Vivo Studies on Human Trunk Load Partitioning and Stability in Isometric Forward Flexions
,”
J. Biomech.
,
39
(
3
), pp.
510
521
.10.1016/j.jbiomech.2004.11.030
16.
,
R.
, and
Wise
,
S. P.
,
2005
,
The Computational Neurobiology of Reaching and Pointing: A Foundation for Motor Learning
,
MIT Press
,
Cambridge, MA
.
17.
Jordan
,
M. I.
, and
Wolpert
,
D. M.
,
1999
, “
Computational Motor Control
,”
The Cognitive Neuroscience
,
M.
Gazzaniga
, ed.,
MIT Press
,
Cambridge, MA
, pp.
601
620
.
18.
Li
,
W.
,
2006
,
Optimal Control for Biological Movement Systems
,
Ph.D. thesis, University of California
,
San Diego
, CA.
19.
Nelson
,
W. L.
,
1983
, “
Physical Principles for Economies of Skilled Movements
,”
Biol. Cybern.
,
46
(
2
), pp.
135
147
.10.1007/BF00339982
20.
Todorov
,
E.
,
2005
, “
Stochastic Optimal Control and Estimation Methods Adapted to the Noise Characteristics of the Sensorimotor System
,”
Neural Comput.
,
17
(
5
), pp.
1084
1108
.10.1162/0899766053491887
21.
Arbib
,
M. A.
,
2003
,
The Handbook of Brain Theory and Neural Networks
,
MIT Press
,
Cambridge, MA
.
22.
Kim
,
J.
, and
Hemami
,
H.
,
1998
, “
Coordinated Three-Dimensional Motion of the Head and Torso by Dynamic Neural Networks
,”
IEEE Trans. Syst. Man Cybern., Part B: Cybern.
,
28
(
5
), pp.
653
666
.10.1109/3477.718516
23.
Pomero
,
V.
,
Lavaste
,
F.
,
Imbert
,
G.
, and
Skalli
,
W.
,
2004
, “
A Proprioception Based Regulation Model to Estimate the Trunk Muscle Forces
,”
Comput. Methods Biomech. Biomed. Eng.
,
7
(
6
), pp.
331
338
.10.1080/1025584042000327115
24.
Perez
,
M. A.
,
1999
, “
Empirical Evaluation of Models Used to Predict Torso Muscle Recruitment Patterns
,” M.Sc. thesis, Virginia Polytechnic Institute and State University, Blacksburg, VA.
25.
Ross
,
E. C.
,
1991
,
The Effect of Resistance Level on Muscle Coordination Patterns and Truncal Velocity, Acceleration, and Deceleration During Isoinertial Trunk Extension
,
Ph.D. thesis, New York University
, New York.
26.
Ross
,
E. C.
,
Parnianpour
,
M.
, and
Martin
,
D.
,
1993
, “
The Effects of Resistance Level on Muscle Coordination Patterns and Movement Profile During Trunk Extension
,”
Spine
,
18
(
13
), pp.
1829
1838
.10.1097/00007632-199310000-00019
27.
Hemami
,
H.
,
2002
, “
Towards a Compact and Computer-Adapted Formulation of the Dynamics and Stability of Multi Rigid Body Systems
,”
J. Autom. Control
,
12
(
1
), pp.
64
70
.10.2298/JAC0201064H
28.
Schultz
,
A. B.
,
,
G. B. J.
,
,
K.
,
Ortengren
,
R.
,
Nordin
,
M.
, and
Bjork
,
R.
,
1982
, “
Analysis and Measurement of Lumbar Trunk Loads in Tasks Involving Bends and Twists
,”
J. Biomech.
,
15
(
9
), pp.
669
675
.10.1016/0021-9290(82)90021-5
29.
Zajac
,
F. E.
,
1989
, “
Muscle and Tendon: Properties, Models, Scaling, and Application to Biomechanics and Motor Control
,”
Crit. Rev. Biomed. Eng.
,
17
(
4
), pp.
359
411
.
30.
Cholewicki
,
J.
, and
McGill
,
S. M.
,
1996
, “
Mechanical Stability of the in Vivo Lumbar Spine: Implications for Injury and Chronic Low Back Pain
,”
Clin. Biomech. (Bristol Avon)
,
11
(
1
), pp.
1
15
.10.1016/0268-0033(95)00035-6
31.
McGill
,
S. M.
, and
Norman
,
R. W.
,
1986
, “
Partitioning of the L4-L5 Dynamic Moment Into Disc, Ligamentous, and Muscular Components During Lifting
,”
Spine
,
11
(
7
), pp.
666
678
.10.1097/00007632-198609000-00004
32.
Nussbaum
,
M. A.
, and
Chaffin
,
D. B.
,
1998
, “
Lumbar Muscle Force Estimation Using a Subject-Invariant 5-Parameter EMG-Based Model
,”
J. Biomech.
,
31
(
7
), pp.
667
672
.10.1016/S0021-9290(98)00055-4
33.
Hatze
,
H.
,
1977
, “
A Myocybernetic Control Model of Skeletal Muscle
,”
Biol. Cybern.
,
25
(
2
), pp.
103
119
.10.1007/BF00337268
34.
Cheng
,
E. J.
,
Brown
,
I. E.
, and
Loeb
,
G. E.
,
2000
, “
Virtual Muscle: A Computational Approach to Understanding the Effects of Muscle Properties on Motor Control
,”
J. Neurosci. Methods
,
101
(
2
), pp.
117
130
.10.1016/S0165-0270(00)00258-2
35.
Van Herp
,
G.
,
Rowe
,
P.
,
Salter
,
P.
, and
Paul
,
J. P.
,
2000
, “
Three-Dimensional Lumbar Spinal Kinematics: A Study of Range of Movement in 100 Healthy Subjects Aged 20 to 60+ Years
,”
Rheumatology
,
39
(
12
), pp.
1337
1340
.10.1093/rheumatology/39.12.1337
36.
Jazbi
,
S. A.
,
Marjovi
,
A.
,
Lucas
,
C.
, and
Ghafoorifard
,
M. H.
,
1998
, “
Critic Based Adaptive Fuzzy Controller for SRM
,”
Proceedings of 5th International Conference on Soft Computing and Information/Intelligent Systems
. Vol.
2
, Oct. 16–20, World Scientific, pp.
692
695
.
37.
Abdi
,
J.
,
Khalili
,
G. F.
,
Fatourechi
,
M.
,
Lucas
,
C.
, and
Khaki Sedigh
,
A.
,
2004
, “
Control of Multivariable Systems Based on Emotional Temporal Difference Learning Controller
,”
Int. J. Eng. Trans. Basics
,
17
(
4
), pp.
363
376
.
38.
Neilson
,
P. D.
, and
Neilson
,
M. D.
,
2005
, “
Motor Maps and Synergies
,”
Hum. Mov. Sci.
,
24
(
5–6
), pp.
774
797
.10.1016/j.humov.2005.09.008
39.
Neilson
,
P. D.
, and
Neilson
,
M. D.
,
2005
, “
An Overview of Adaptive Model Theory: Solving the Problems of Redundancy, Resources, and Nonlinear Interactions in Human Movement Control
,”
J. Neural Eng.
,
2
(
3
), pp.
S279
S312
.10.1088/1741-2560/2/3/S10
40.
Prochazka
,
A.
,
1999
, “
Quantifying Proprioception
,”
Peripheral and Spinal Mechanisms in the Neural Control of Movement
,
M.
Binder
, and
J.
McDonnaugh
, eds.,
Elsevier
,
Amsterdam, NY
, pp.
133
142
.
41.
Dariush
,
B.
,
Parnianpour
,
M.
, and
Hemami
,
H.
,
1998
, “
Stability and a Control Strategy of a Multilink Musculoskeletal Model With Applications in FES
,”
IEEE Trans. Biomed. Eng.
,
45
(
1
), pp.
3
14
.10.1109/10.650346
42.
Haykin
,
S. S.
,
1999
,
Neural Networks: A Comprehensive Foundation
,
Prentice-Hall
,
.
43.
Oddsson
,
L.
, and
Thorstensson
,
A.
,
1987
, “
Fast Voluntary Trunk Flexion Movements in Standing: Motor Patterns
,”
Acta Physiol. Scand.
,
129
(
1
), pp.
93
106
.10.1111/j.1748-1716.1987.tb08044.x
44.
Sergio
,
L. E.
,
Hamel-Paquet
,
C.
, and
,
J. F.
,
2005
, “
Motor Cortex Neural Correlates of Output Kinematics and Kinetics During Isometric-Force and Arm-Reaching Tasks
,”
J. Neurophysiol.
,
94
(
4
), pp.
2353
2378
.10.1152/jn.00989.2004
45.
Berardelli
,
A.
,
Hallett
,
M.
,
Rothwell
,
J. C.
,
Agostino
,
R.
,
Manfredi
,
M.
,
Thompson
,
P. D.
, and
Marsden
,
C. D.
,
1996
, “
Single-Joint Rapid Arm Movements in Normal Subjects and in Patients With Motor Disorders
,”
Brain
,
119
(
Pt 2
), pp.
661
674
.10.1093/brain/119.2.661
46.
Britton
,
T. C.
,
Thompson
,
P. D.
,
Day
,
B. L.
,
Rothwell
,
J. C.
,
Findley
,
L. J.
, and
Marsden
,
C. D.
,
1994
, “
Rapid Wrist Movements in Patients With Essential Tremor. The Critical Role of the Second Agonist Burst
,”
Brain J. Neurol.
,
117
(
Pt 1
), pp.
39
47
.10.1093/brain/117.1.39
47.
Zatsiorsky
,
V. M.
,
Li
,
Z.-M.
, and
Latash
,
M. L.
,
1998
, “
Coordinated Force Production in Multi-Finger Tasks: Finger Interaction and Neural Network Modeling
,”
Biol. Cybern.
,
79
(
2
), pp.
139
150
.10.1007/s004220050466
48.
Kim
,
J.
, and
Hemami
,
H.
,
1995
, “
Control of a One-Link Arm by Burst Signal Generators
,”
Biol. Cybern.
,
73
(
1
), pp.
37
47
.10.1007/BF00199054
49.
Mussa-Ivaldi
,
F. A.
,
Hogan
,
N.
, and
Bizzi
,
E.
,
1985
, “
Neural, Mechanical, and Geometric Factors Subserving Arm Posture in Humans
,”
J. Neurosci.
,
5
(
10
), pp.
2732
2743
.
50.
Golkhou
,
V.
,
Parnianpour
,
M.
, and
Lucas
,
C.
,
2005
, “
Neuromuscular Control of the Point to Point and Oscillatory Movements of a Sagittal Arm With the Actor-Critic Reinforcement
,”
Comput. Methods Biomech. Biomed. Eng.
,
8
(
2
), pp.
103
113
.10.1080/10255840500167952
51.
McGill
,
S. M.
,
Norman
,
R. W.
, and
Cholewicki
,
J.
,
1996
, “
A Simple Polynomial That Predicts Low-Back Compression During Complex 3-D Tasks
,”
Ergonomics
,
39
(
9
), pp.
1107
1118
.10.1080/00140139608964532
52.
Bazrgari
,
B.
,
,
A.
,
Trottier
,
M.
, and
Mathieu
,
P.
,
2008
, “
Computation of Trunk Equilibrium and Stability in Free Flexion-Extension Movements at Different Velocities
,”
J. Biomech.
,
41
(
2
), pp.
412
421
.10.1016/j.jbiomech.2007.08.010
53.
Arjmand
,
N.
,
Plamondon
,
A.
,
,
A.
,
Larivière
,
C.
, and
Parnianpour
,
M.
,
2011
, “
,”
J. Biomech.
,
44
(
1
), pp.
84
91
.10.1016/j.jbiomech.2010.08.028
54.
Sarti
,
M. A.
,
Lison
,
J. F.
,
Monfort
,
M.
, and
Fuster
,
M. A.
,
2001
, “
Response of the Flexion-Relaxation Phenomenon Relative to the Lumbar Motion to Load and Speed
,”
Spine
,
26
(
18
), pp.
E421
E426
.10.1097/00007632-200109150-00019