Simulation of Movement in Three-Dimensional Musculoskeletal Human Lumbar Spine Using Directional Encoding-Based Neurocontrollers

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.

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

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

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:

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

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.

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:

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.

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

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.

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.

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

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

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.

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.

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

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

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.

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

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.

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

The three-burst activation pattern (agonist-antagonist-agonist)² [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].

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.

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.

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.

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.

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.