A 3-D Stability-Based Dynamic Computational Model of Human Trunk Movement: Towards the Development of a Paradigm for Forensic Spinal Injury Biomechanical Analysis

Over the past few decades, biomechanical models of human movement have emerged as important tools in the investigation of possible injury pathways towards insurance claims and forensic applications. While the literature has some examples of biomechanical injury models that were successfully used to establish whether a particular injury is the consequence of an accident versus an assault [1-3], the application of biomechanical simulation and modeling techniques in forensic injury biomechanics is not yet well developed and therefore not commonly used. On the other hand, the field of forensic biomechanics in particular stands to benefit substantially from the use of time and cost-effective models, not only in terms of the flexibility in reconstructing realistic, biomechanically-correct crime and injury scenes, but also in answering a whole range of “what if ” types of questions that evaluate the credibility of a hypothesis based on solid biomechanical quantitative data.


Introduction
Over the past few decades, biomechanical models of human movement have emerged as important tools in the investigation of possible injury pathways towards insurance claims and forensic applications. While the literature has some examples of biomechanical injury models that were successfully used to establish whether a particular injury is the consequence of an accident versus an assault [1][2][3], the application of biomechanical simulation and modeling techniques in forensic injury biomechanics is not yet well developed and therefore not commonly used. On the other hand, the field of forensic biomechanics in particular stands to benefit substantially from the use of time and cost-effective models, not only in terms of the flexibility in reconstructing realistic, biomechanically-correct crime and injury scenes, but also in answering a whole range of "what if " types of questions that evaluate the credibility of a hypothesis based on solid biomechanical quantitative data.
Acute injuries in biomechanics are typically assumed to be the result of two major mechanisms: structural stability failure such as buckling, or material failure when the imposed level of stresses/strains is higher than the ultimate stresses and strains of the respective tissues. The field of occupational biomechanics has contributed to the understanding of another important mechanism of injury or musculoskeletal disorder: cumulative trauma disorder or repetitive strain injury. A low magnitude of load applied repetitively can lead to tissue fatigue, loss of stiffness, and residual strains. The inability of biological tissue to remodel/ repair the damaged tissue fast enough (i.e. the injury cycle outpacing the remolding cycle) slowly leads the micro-failures to accumulate/ propagate and allows many disorders or disabling diseases to manifest themselves leading to pain, inflammation, and to some degree of disability. Biomechanical modeling can predict the stresses and strains on the tissue once the observed or expected kinematics are input into an appropriate inverse dynamics models However, when the forces and moments are the known quantities, a forward dynamics model must be integrated to provide the kinematic results. The latter process is more complex and requires more sophisticated computational tools. When the system is kinematically and/or kinetically redundant, the modeling effort is further complicated, as the redundancies must be addressed as well. Optimal control theory is used to solve these boundary value problems; and due to the associated complexity, computational procedures are usually used instead of analytical methods based on calculus of variation.
The human body has more degrees of freedom than the minimum needed to perform various movements. Such redundancy provides the body with flexibility in coordinating these degrees of freedom and in performing tasks with efficiency [4,5]. This inherent redundancy also makes the human musculoskeletal system a suitable candidate for the use of optimal control strategies in the prediction of motion patterns, such as investigating the strategies underlying the planning and execution of human spine movements [6]. Deterministic optimization minimizes a performance index subject to a set of dynamic equations and boundary conditions over the motion interval, allowing the determination of the optimal movement profile(s). However, existing performance indices for investigating movements of the human spine are either too simplistic, resulting in nonrealistic solutions, or highly nonlinear leading to complex and computationally inefficient optimization problems [6].
A potential approach for minimizing the performance index is parameter optimization. In this approach, the nonlinear cost function of the optimization problem is transformed into a linear cost function based on the coefficients of a series expansion written for each degree of freedom [6][7][8]. Parnianpour et al. [7] parameterized one degree of freedom with a fifth-order polynomial and a Fourier expansion. Biess et al. [6] minimized a nonlinear performance index subject to different boundary conditions in the case of point-to-point and rhythmic movements of the hand. They utilized a formulation derived from a time series to satisfy the imposed boundary conditions and a series of a set of basic functions. A similar approach is utilized in this study in order to minimize three distinct nonlinear performance indices subjected to two-point boundary conditions, at the beginning and end of the trunk motion.
Spinal instability has been indicated as one of the most important risk factors in low back disorders [9][10][11][12][13][14][15][16]. Granata and Wilson [14] reported that the majority of low back injuries occur when compressive loads exceed 3400N: however, spinal instability may cause that threshold value to decrease dramatically to only 88 N. From a forensic injury point of view, spinal instability can therefore be a critical element in injury assessment and quantification.
Granata et al. investigated the influence of trunk posture [14], exertion direction [17,18], and preload [19] on spinal stability during flexion and extension exertions. Shirazi-Adl et al. [20,21] developed a multi-degree of freedom finite element model to quantify the role of passive and active tissues in spinal stability. Other studies [22,23] evaluated the effect of reflex and intrinsic mechanisms on spinal stability by utilizing nonlinear system identification methods. It was found that around 62% of the stiffness is due to the reflex mechanism [23].
However, the studies mentioned above investigated the stability of the spine only in the static postures, without consideration of the influence of acceleration and velocity. Recently, Tanaka et al. [24][25][26] evaluated the basin of stability using the Finite Time Lyapunov Exponent method. Granata and England [27] demonstrated that fast paced repetitive trunk movements diverged more quickly than slower The CNS neuromuscular intrinsic mechanism generates sufficient co-activation to stabilize the spine without any time delay associated with feedback systems [28], and can provide essential joint stiffness to partially stabilize the spine during movement. Recently, Zeinali et al. [28] assumed that stability-based optimization models can provide intrinsic stiffness to partially stabilize the spine during movement from the upright position to 60 degrees of flexion, and showed that this assumption would result in the co-activation of agonistic and antagonistic muscles. However, it was demonstrated that the intrinsic mechanism alone couldn't stabilize the spine during its motion [23]. The reflex mechanism would also increase the joint stiffness of the spine by providing appropriate feedback in the presence of disturbances [11,16,23] This paper presents a computational model of trunk performance for use in multiple applications, including forensic injury biomechanics. The model simulates discrete (point-to-point) trunk sagittal movements, allowing the quantitative assessment of joint reaction forces, considered a main causal risk factor of spinal injuries and low back disorders. The relative importance of the co-activation strategy and the particular motion syntheses strategy are also investigated. The proposed 3-D modular model is dynamic and has the capability to predict the optimal trajectories with respect to three highly prevalent and physiologically credible cost functions.

Formulation of parameter optimization
Performance indices (Table 1) typically suggested for simulating human spine movements are highly nonlinear, resulting in a corresponding optimization problem with a high computational cost and complexity. Therefore, this study adopted the strategy of parametric optimization (converting a nonlinear cost function to nonlinear functions of expansion coefficients through parameterized expression of each degree of freedom) [6]. The set of functions utilized in the expansion should form a complete set of orthogonal functions over the interval of movement. The completeness of this set is a necessary condition to be satisfied if any function, square integrable in the interval of motion, was expressed in terms of expansion coefficients. In addition, the boundary conditions imposed by the movement should be satisfied by these sets of functions [6]. The goal of nonlinear optimization is to find the expansion coefficients, while satisfying the equation of motion and the physiological boundary conditions included in the model.  Figure 1 shows the basic functions that are used in (1).
Nonlinear objective functions listed in Table 1 are transformed into nonlinear objective functions in terms of cik by substituting Equation (1) into the objective function. In this study, three objective functions were considered in order to model the optimal trajectory of trunk movement in the sagittal plane: Minimum-Peak Torque model, Minimum-Jerk model and Minimum-Energy model. The following boundary conditions were considered at the beginning and end of motion

Dynamics of the spine
A 3D inverted pendulum with a ball and socket joint at its base along with 18 muscle fascicles was utilized to model the trunk system. The dynamic equations of motion in a body coordinate system can be expressed as: Where I is the matrix of moments of inertia, W is the angular velocity vector in body coordinate system,τ is torque about the lumbosacral joint in the body coordinate system which is generated by muscles, Θ is the angular position vector in an inertial coordinate system, G is the gravitational moment vector and WW is the skew symmetric matrix based on W. Parameters utilized in this study are summarized in Table 2.

Muscle model
It was assumed that trunk motion is generated by five main trunk muscle groups: left and right Internal Oblique (IO), External Oblique (EO), Erector Spinae (ES), Latissimus Dorsi (LD) and Rectus Abdominis (RA). Some muscles were further divided into segments, which were modeled independently. For example, the Erector Spinae is divided into Longissimus Thoracic (LT) and Iliocostalis Lumborum (IL). Based on Hill's model, muscle forces can be expressed in terms of muscle length l, muscle stretch velocity , and muscle activation level  as [28][29][30]: where f is the muscle force, max f is the maximum muscle force, f l (l) is force-length relationship, f v (i) is force-velocity relationship, and f p (l) is passive force of muscles. The anatomical geometry of muscles is based on Zeinali et al. [28].

Muscle neural activation levels
IIn addition to kinematic redundancy, which allows the CNS to execute each task using different trajectories based on an inventory of infinite kinematic profile (velocity and acceleration) combinations, kinetic redundancy offers infinite muscle activation pattern possibilities to produce the same net moment profile [28]. Hence, static optimization is a powerful tool that could be effectively used to predict muscle recruitment patterns by assuming a proper cost function In this study, static optimization is used to calculate the neural activation level of muscles for a desired moment at each instant of the predicted motion. By minimizing the norm of the muscles activation level, the cost function is represented as: This function is subject to the following physiological constraints on muscle activation and the dynamic equations of motion: Where B is transformation matrix between body coordinate system and inertial coordinate system, θ ∂ ∂L is the matrix of muscle moment arms, and F is the vector of muscle forces. Values of activation level, α, can range between 0 and 1 [28][29][30].

Stability constraint
In order to ensure the stability of the spine during trunk motion, another condition must be included in the static optimization algorithm. Based on the 'equilibrium point' hypothesis [31,32], it was assumed that the spine moves by following a trajectory comprised of a sequence of equilibrium points. Therefore, the trajectory between the beginning and end of a spinal movement is stable provided that the points in this trajectory demonstrate stable equilibrium behavior. The potential energy of the spine is considered as the sum of potential energy of each muscle as well as gravitational energy:: The moment about the lumbosacral joint is calculated using the gradient of potential energy with respect to the angular position: The terms of the joint stiffness matrix can be calculated as the derivative of the moment with respect to joint angle Θ: Where K is the matrix of joint stiffness about the lumbosacral joint. This equation can be expressed in terms of muscle force and muscle stiffness as follows: is a second order term, which can be neglected.
Hence, joint stiffness can be expressed as: 18 18 18  (11) Where k i is the stiffness of the i th muscle.
In many studies on spinal stability, it was assumed that muscle stiffness changes linearly with muscle force (see, for example, [9]) That is: Where q is a constant. If Θ e corresponds to the local minimum of the lumbosacral joint total potential energy, the variation of total potential energy of the system at a stable equilibrium point must be positive [32]. In particular, an equilibrium point is stable if the matrix of joint stiffness is negative definite. This additional condition incorporated with static optimization results in muscle co-activation prediction.

Spindle model
The role of the reflex mechanism in spinal stability has been investigated by many researchers [16,22]. These studies assumed that the spindle has a significant role in controlling human trunk posture and movement. Gielen et al. [33] developed the following nonlinear model to simulate the role of spindles: Where R is the firing rate of spindles, G s is the gain of spindles, l and i are the actual length and velocity of muscles in each instance, respectively, and l d is the muscle length computed from the predicted optimal trajectories. Feedback neural activation of muscles can be computed based on a linear mapping from firing rates and typically ranges between 0 and 0.35 [30,34].

Joint reaction force
The reaction forces at the various joints, including the lumbosacral joint, are generally considered critical in injury and/or low back considerations [28][29][30]. Hence, the determination of such forces is one of the main goals of this paper. Reaction forces can be expressed in terms of their components: compression force, anterior shear force, and lateral shear force. Compression and shear forces are considered as the main risk factors leading to low back pain [30].
The joint reaction force is calculated as: Where F r is the reaction force vector at the L5/S1 joint, m is the concentrated mass at the center of mass, F muscle is the summation of individual muscle force vectors, F g is the force vector arising from gravity, and a is the linear acceleration of the concentrated mass which is determined from: In above equation, L is the position vector of the center of mass.

Numerical simulations
Five different simulation cases were run to investigate the effect of stability constraints on the neural activation of muscles (Table 3). In the first simulation (Case 1), the distribution of moment between muscles was studied in the absence of the stability constraint, while Case 2 included these constraints. Case 3 analyzed the effect of changing the value of the parameter q on the joint stiffness and stability requirements, where the joint stiffness is compared with the minimum required stiffness due to the stability constraint. Case 4 explored the effect of different spindle gains G s (0, 100, 150, 200, and 600), as well as the transmission time delay in the presence of an external disturbance of -30 N.m at the time of 0.15s for a 0.07s duration. This was done to demonstrate the role of spindles in spinal stability in the scenario where the stability constraint is not included in the simulation of trunk movement. Another important issue that affects spinal stability is the presence of time delay in the reflex mechanism. Franklin et al. reported that an increase in time delay would reduce spinal stability. Therefore, Case 5 analyzed the effect of transmission delays (0, 15, and 30 ms), while G s was fixed at 150, applying the same perturbation as in Case 4.
The compression and shear joint reaction forces at L5/S1 were computed for Cases 1 and 2 and presented in a spine body coordinate system to evaluate the effects of the stability constraint and cost functions on the internal loading of the spine. All simulations were coded and implemented in MATLAB 7 (The Mathworks, TM Massachusetts).

Results
The three objective functions considered in this study resulted in distinct optimal trajectories of trunk motion in the sagittal plane ( Figure 2). The summary statistics of the kinematic and kinetic measures (angular velocity, acceleration, net muscular torque and work) are presented in Table 4. Muscle activation and net joint reaction forces are reported in Table 5. Values of the expansion coefficients and polynomial coefficients determined from the optimization are listed in Tables 6 and 7, respectively. In the Minimum Jerk Model, the velocity profile of the optimal trajectory is bell shaped and the peak value occurs in the middle of the time interval, 500 ms. On the other hand, in the Minimum-Peak Torque and Minimum-Energy models, the peak values of the velocity profiles occurred earlier, at 320 ms, and later, at 740ms, respectively. This is in agreement with literature [7].
In Case 1 (no stability constraint), the neural activation profiles of the right flexors and right extensors for the three cost functions are shown in Figure 3. The recruitment pattern of the left back muscles is  similar to those found for the right ones since the anatomical muscle models used are symmetric about the mid-sagittal plane [35][36][37][38][39][40][41][42][43][44][45][46]. As expected, the trunk flexors are activated to launch and accelerate trunk motion until the net muscular torque becomes zero (Figure 2(D)). Among the trunk flexor muscles, RA had the largest contribution during this phase of motion. In the following phase, the trunk flexors cease activity, while the trunk extensors are recruited to decelerate the trunk motion and maintain it at 60 degrees of flexion from the upright position. In other words, IL and LT are the extensors most activated in order to slow down the trunk motion.
Activation profiles determined for the three models differ with respect to the time at which the flexors terminate their action. For instance, the termination is found to occur earlier in the case of the Minimum-Peak Torque model. Another important difference is the magnitude of the net muscular torque ( Figure 2D). The maximum value computed for the Minimum-Energy model was 184N.m, which is much larger than the 140N.m and 120N.m values predicted by the Minimum-Jerk and Minimum-Peak Torque models, respectively (approximately 24% and 35% differences). This implies that in the Minimum-Energy model, the, extensors are activated more in the second phase of motion to decelerate from peak acceleration until LT has reached its maximum level of activation at the time 0.8s. As the neural activation of muscles is biologically limited between 0 and 1, other extensors are activated more when LT reaches its highest level of activation. The evolution of muscle forces for the different models is shown in Figure 4. Force patterns of IO2 and EO1 are similar in spite of the higher level of neural activation of IO2 (Figure 3). Figures 3D through 3F depict the co-activation of agonistic and antagonistic muscles over parts of the spinal motion when the stability constraint, expressed by Equation (11), is included in the optimization routine (Case 2)..
The sensitivity of joint stiffness to the q parameter (Equation (10)) is illustrated by Figure 5 (Case 3). Increasing q results in a decrease in stiffness of the lumbosacral joint. Therefore, if q 4, co-activation is required in order to maintain the joint stiffness higher than the minimum stiffness associated with spinal stability. Figures 5B, 5D and 5F show the minimum joint stiffness as well as the joint stiffness computed for q=4 in the presence of stability constraints for each of the three models, respectively. The time and duration of co-activation during the simulated flexion movement can be determined from the stability requirements. Figures 6A and 6B, respectively, demonstrate the effect of an external perturbation on the actual angular position and the corresponding deviation of velocity profile from the optimal trajectory in the Minimum-Jerk model. The plots are shown assuming a constant time delay of 15ms in the reflex mechanism (Case 4). If the model does not include spindles, these deviations decrease to 0.05 rad and 0 rad/s at the end of the time interval, respectively. Surprisingly, higher spindle gains (>150) resulted in larger errors in the angular position and velocity profiles. Consequently, spindle gain was set equal to 150 for Case 5 of the simulation. However, the gain setting should be tuned for each time delay. Figure 7 shows the results of a parametric analysis carried out for different values of time delay. As expected, the error in the angular position profile observed at the end of the simulation time interval increased from -0.05 to -0.1 rad.
The joint reaction force profiles generated during flexion (Cases 1 and 2) are shown in Figure 8 for the three cost functions. A higher compressive force was predicted due to the co-activation when the stability constraint was included. The effect of the cost function on the internal loading was more significant than the presence of the stability constraint (Table 5 and Figure 8). Higher muscle activities and joint reaction forces can lead to muscle fatigue and hence a higher risk of intervertebral disc degeneration.

Discussion and Conclusions
The human trunk is a complex system characterized by numerous degrees of freedom. Various computational models at different levels of complexity were developed to investigate the kinematic and dynamic behaviors of the trunk. Multiple degrees of freedom models [13,35] were utilized to analyze spinal stability and stress in static postures. These models, however, have limited ability for addressing the recruitment pattern(s) of muscles [36]. Adequate rigid-body models on the other hand are capable of representing both the stability of the spinal column as well as the complex muscular coordination and recruitment patterns. Hence, a rigid body inverted pendulum model was utilized in this study to investigate the neuromuscular characteristics of the spine.
Three distinct objective functions with different time to peak velocity values were considered in order to simulate the optimal trajectory of the spine from upright position to 60 degrees of flexion. The three criteria (Minimum Jerk, Minimum Peak Torque and Minimum Energy) were selected to realistically emulate the functional biomechanics of the spine while maintaining adequate computational efficiency. Simulation results show that the Minimum Energy model yields the largest peak compressive force ( Table 5). The mean compressive joint reaction forces were respectively 2066.3N, 1847.1N, and 1637.1N as predicted by the Minimum Peak Torque, Jerk and Energy models, respectively, without imposing the stability constraints (Table 5). Parnianpour et al. considered the influence of various cost functions on the trunk motion including jerk, work, impulse, energy and peak torque [7][8]. They used these cost functions to evaluate the effects of strength impairment on the performance of the trunk. In the present study, the optimal trajectories were computed based on a complete set of basic functions, and were in good agreement with optimal profiles obtained by summing over a fifth order time polynomial and a linear combination of Fourier terms [7][8].
The Minimum Jerk model used a kinematic based cost function based on the literature [31], and hence predicts the optimal trajectory regardless of the effect of spinal dynamics and external forces. Consequently, it always has a bell-shaped velocity profile. Conversely, the Minimum Peak torque and Minimum Energy models are dynamic models based on the minimum commanded model utilized by Nakano et al. [37]. It can be seen from Figure 2 that the time to peak velocity changed in these models according to the net muscular torque profiles.
In the Minimum Energy model, for example, peak velocity occurred after the middle of the time interval as a result of the constraint to reduce the integral of the square of torque over the motion time interval. On the other hand, in the Minimum Peak torque model, the maximum velocity occurred before the middle of the time interval, most likely in order to decrease the maximum value of the net muscular torque during the second phase of motion (Figure 2(D)).
The literature has few experimental studies investigating the strategies underlying the execution of trunk flexion and patterns of muscular activity towards quantitative validation of the results in the current study. Oddsson et al. [38] asked volunteers to perform 3 series of fast flexion movements of the trunk in the sagittal plane with successively increasing the movement amplitude. They measured the EMG activity of two groups of trunk muscles: RA and ES. They reported an initiating flexor (agonist) activity during the accelerating part of the trunk flexion followed by a braking burst in the extensor (antagonist) activity, with flexor activity terminating during the decelerating part of the movement. This behavior is consistent with the numerical simulation results obtained in this study. The predicted mono-modal velocity profiles of the trunk are also similar to those observed by Ross et al. [39] in their investigation of trunk extension against various isoinertial resistance levels.
The analysis of force profiles at L5/S1 in the presence of a stability constraint (Table 5 and Figure 8) shows that the average and peak values of the compression force were 2128.3N and 2962.6N, respectively, during the Minimum Peak torque model, 1954.5N and 3313.7N, for the Minimum Jerk model, and 1792.1N and 4308.3N in the Minimum Energy model. In the Minimum Energy model, the extensors were activated more than in the other two models until LT reached the highest level of neural activation ( Figure 3F). However, the average muscle activation during motion was less than the other models, which may result in less muscle fatigue. The general trend for the reaction forces found in this study is consistent with the results by McGill et al. [40]. However, the acceleration profiles assumed in this study yield larger compression forces than those found by that group [40].
Many studies have attempted to evaluate the role of the intrinsic mechanism in the stabilization of human posture and movement. Some studies considered stability-based optimization as an appropriate tool [14,28,[41][42] and suggested that the eigenvalues of the state matrix in the spinal motion linearized equation must be always negative to maintain stability. The stability constraint implemented in this study is a conservative approximation, which underestimates the actual stability of the spinal system. In reality, impedance is related not only to the stiffness of the passive structures and muscles, but also to the damping effects of the viscous time-dependent response of the passive structures as well as the velocity-tension relationship of the active muscle.
Granata et al. [14] predicted co-contraction of agonistic and antagonistic muscles by setting the eigenvalues of the Hessian matrix representing the spinal musculoskeletal potential energy to positive. Min-Peak Torque Min-Jerk Min-Energy

Case 4
Min-Jerk _ + + Variable spindle gains G s ; transmission time delay of 15ms.

Case 5
Min-Jerk _ + + Variable transmission time delays; G s =150.  Howarth et al. [42] compared three methods to evaluate spinal stability based on the Hessian matrix, and observed that while all the methods have similar interpretation from the point of view of spinal stability, the predicted muscle recruitment patterns differ between models.
The abovementioned studies considered the static stability of the spine. Recently, few studies investigated the dynamic stability using different approaches such as Lyapunov stability criteria [43][44][45], finite time Lyapunov exponents [25][26][27]46], and stability diffusion analysis [47]. These studies, however, only considered dynamic spinal stability in the close vicinity of the upright position or during repeating motion [27]. The present study assumes that the trunk moves in a trajectory created by static equilibrium points by the neuromusculoskeletal system. Hence, each equilibrium point is stable if the joint stiffness is more negative than the rate of change of the gravitational torque.

Min-Peak Torque Model
Min-Jerk Model Mean-Energy Model

Movement Pattern
Velocity(rad/s) Mean(SD) -   Many investigations have attempted to evaluate the role of trunk muscles during various tasks. In the present work, the contribution of each muscle to the L5/S1 joint stiffness depends on its force level, as well as its squared moment arm (Equation (9)). Hence, the required stiffness is mainly provided by muscles with longer moment arms. In the case of the Minimum Jerk model, Figure 9 shows that EO1 and IO1 have the largest moment arms among all muscles during spinal motion, when the joint stiffness stability criterion is not enforced. Furthermore, it can be seen that moment arms of the muscles depend on the current joint angle, i.e. the orientation of the muscles determines their contribution to spinal stability [48]. Cholewicki et al. [49] considered the relative contribution of trunk muscles to spinal stability during various isometric exertions. They found that the elimination of any muscle could disturb the spinal stability by more than 30%, and concluded that all trunk muscles, in fact, provide an important contribution to the stability during various types of exertions, such as flexion, extension, and lateral bending. They also showed that the elimination of the External Oblique would decrease the stability index more than the other trunk muscles during the flexion exertion. The results of the present study are in good agreement with these findings.
While intrinsic mechanisms can increase the spinal joint stiffness in a feed forward manner without time delay to ensure spinal stability in the presence of small perturbations, previous studies [23,50] showed that the reflex system is necessary to guarantee that stability. In particular, Moorhouse et al. [23] reported that 62% of the required joint stiffness for spinal stability is produced by the reflex mechanism. In this study, spindles provided the required feedback to control spinal motion. The results showed that spinal stability is increased with higher spindle gain values ( Figure 7) and decreased due to higher time delays in the reflexive system ( Figure 8). These findings are consistent with the results of Franklin et al. [50], who showed that the maximum tolerable time delay values decrease with high reflex gains when stability is assumed. It should be mentioned here that the main advantage of the reflex mechanism is a decrease in the energy expenditure in comparison to other scenarios using the intrinsic mechanism alone [50].
The present study did not address certain aspects of the spinal neuromuscular system due to the associated complexity. First, the     excitation-activation muscle contraction dynamics were neglected, as they would have required complex ordinary or partial differential equations. Due to the stability-based static optimization methodology adopted, the model developed here has limited capacity for exerting differential constraints on the neural activation levels of muscles. This issue should therefore be investigated in future work, perhaps by adopting more complicated optimization algorithms such as direct collocation [51].
Since static optimization was performed, the second-order term in Equation (10) was neglected in this study. The effect of this term should be investigated in a future dynamic stability assessment of the spine.
A crude stability constraint was included in order to incorporate the stability of the spine in the model during a point-to-point movement. It is theoretically significant to consider the validity of this restrictive requirement. In other words, one may question why they should the neuromusculoskeletal system spend energy and increase wear and tear on the passive structures throughout the movement in order to meet the stability criteria? Control engineers would rather design a controller that moves the plant to a target and design the control such that the target is a point attractor. This issue must also be addressed in future work, as it is not known a priori is which part of the movement is instability encountered due to impending perturbations. It can be speculated that if previous experiences indicate a safe perturbation-free environment, the stability constraint may be relaxed, allowing more focus on the vicinity of the target so that the target itself is not overshot. However, if the environment is considered unsafe and an external perturbation is expected to be encountered, a minimum level of joint stiffness, depending on the strength of the prior perturbing forces, is likely to be invoked. In a probabilistic sense, the present simulations are closer to the latter case since a minimum joint stiffness is defined. This was confirmed during static tasks by Rashedi et al. [52].
The 'equilibrium point' hypothesis [31,53,54] is used to describe static posture or dynamic movement. However, the concept of stiffness must be expanded in order to include the orientation of the spine, distribution of its inertia, and dynamics toward the control of spinal impedance. Specifically, impedance depends on stiffness, viscosity  and inertia. The relative contributions of these terms to the overall mechanical impedance depend on the frequency of motion in addition to their relative magnitude [53][54][55]. Hence, developing a more general model of stability constraints will require investigation of impedance modulation rather than joint stiffness modulation during spinal motion. The direction of velocity is another important factor, which should be investigated in future studies, since the velocity-dependent behavior of the musculoskeletal system in response to a perturbation depends on both the direction and magnitude of the spinal velocity vector in equilibrium [54][55][56]. That would require many additional elements currently missing, such as the inclusion of all degrees of freedom in the lumbar spine, the spine-pelvic rhythm, and many more muscles with wrapping elements and via points to account for the curvature of muscles as the trunk changes its orientation.
Although the spinal stability was modeled by including a rather simple additional constraint in the static optimization routine to provide the required intrinsic joint stiffness at the L5/S1 joint, the results of the present study are in a good agreement with those of Zeinali et al. [57] in terms of predicting muscles' recruitment patterns. In this study, 60 degrees of flexion from the upright position was modeled with a stability-based optimization. The results showed that for q=4, the most significant contribution to the spinal stability is provided by EO1.
The muscular stiffness term in Equation (12) was modeled using a linear relationship, in accordance with the study by Bergmark [9]. Many researchers utilized this relationship with an average value of q =10 for to quantify stability [28]. Brown and McGill [48] showed that assuming a linear relationship of the muscle force-stiffness limits the stability consideration. They compared linear and nonlinear relationships for the muscle force-stiffness and concluded that stability would be continuously increased by increasing the muscle force. Furthermore, they found that it reaches its peak value at a specific muscle force in a nonlinear relationship. Therefore, the effect of nonlinearity of the muscle force-stiffness relationship should also be investigated in the future.
In conclusion, a 3-D dynamic optimization computational model that includes 18 anatomically oriented muscles was developed in this study for simulating the motion of the human trunk. The nonlinear optimization problem was formulated using three physiologically plausible performance indices and a stability constraint to predict muscle recruitment patterns in a fully-dynamic task. The current simulation could be used in forensic injury biomechanics applications towards establishing probabilities of injury causation. This is an important steadily growing area in biomechanics that provides biomechanical quantitative data in cases of personal injury, product and premise liability, wrongful death, and criminal cases.