Muscle forces during running predicted by gradient-based and random search static optimisation algorithms

Muscle forces during locomotion are often predicted using static optimisation and SQP. SQP has been criticised for over-estimating force magnitudes and under-estimating co-contraction. These problems may be related to SQP's difficulty in locating the global minimum to complex optimisation problems. Algorithms designed to locate the global minimum may be useful in addressing these problems. Muscle forces for 18 flexors and extensors of the lower extremity were predicted for 10 subjects during the stance phase of running. Static optimisation using SQP and two random search (RS) algorithms (a genetic algorithm and simulated annealing) estimated muscle forces by minimising the sum of cubed muscle stresses. The RS algorithms predicted smaller peak forces (42% smaller on average) and smaller muscle impulses (46% smaller on average) than SQP, and located solutions with smaller cost function scores. Results suggest that RS may be a more effective tool than SQP for minimising the sum of cubed muscle stresses in static optimisation.     

Sensitivity of predicted muscle forces to parameters of the optimization-based human leg model revealed by analytical and numerical analyses

There are different opinions in the literature on whether the cost functions: the sum of muscle stresses squared and the sum of muscle stresses cubed, can reasonably predict muscle forces in humans. One potential reason for the discrepancy in the results could be that different authors use different sets of model parameters which could substantially affect forces predicted by optimization-based models. In this study, the sensitivity of the optimal solution obtained by minimizing the above cost functions for a planar three degrees-of-freedom (DOF) model of the leg with nine muscles was investigated analytically for the quadratic function and numerically for the cubic function. Analytical results revealed that, generally, the non-zero optimal force of each muscle depends in a very complex non-linear way on moments at all three joints and moment arms and physiological cross-sectional areas (PCSAs) of all muscles. Deviations of the model parameters (moment arms and PCSAs) from their nominal values within a physiologically feasible range affected not only the magnitude of the forces predicted by both criteria, but also the number of non-zero forces in the optimal solution and the combination of muscles with non-zero predicted forces. Muscle force magnitudes calculated by both criteria were similar. They could change several times as model parameters changed, whereas patterns of muscle forces were typically not as sensitive. It is concluded that different opinions in the literature about the behavior of optimization-based models can be potentially explained by differences in employed model parameters.     

Lumbar spinal muscle activation synergies predicted by multi-criteria cost function

The hypothesis that control of lumbar spinal muscle synergies is biomechanically optimized was studied by comparing EMG data with an analytical model with a multi-component cost function that could include (1) trunk displacements, (2) intervertebral displacements, (3) intervertebral forces; (4) sum of cubed muscle stresses, and (5) eigenvalues for the first two spinal buckling modes. The model's independent variables were 180 muscle forces. The 36 displacements of 6 vertebrae were calculated from muscle forces and the spinal stiffness. Calculated muscle activation was compared with EMG data from 14 healthy human subjects who performed isometric voluntary ramped maximum efforts at angles of 0 degrees, 45 degrees, 90 degrees, 135 degrees and 180 degrees to the right from the anterior direction. Muscle activation at each angle was quantified as the linear regression slope of the RMS EMG versus external force relationship, normalized by the maximum observed EMG.There was good agreement between the analytical model and EMG data for the dorsal muscles when the model included either minimization of intervertebral displacements or minimization of intervertebral forces in its cost function, but the model did not predict a realistic level of abdominal muscles activation. Agreement with EMG data was improved with the sum of the cubed muscle stresses added to the cost function. Addition of a cost function component to maximize the trunk stability produced higher levels of antagonistic muscle activation at low efforts than at greater efforts. It was concluded that the muscle activation strategy efficiently limits intervertebral forces and displacements, and that costs of higher muscle stresses are taken into account, but stability does not appear to be maximized. Trunk muscles are apparently not controlled solely to optimize any one of the biomechanical costs considered here.     

Evaluation of a particle swarm algorithm for biomechanical optimization

Optimization is frequently employed in biomechanics research to solve system identification problems, predict human movement, or estimate muscle or other internal forces that cannot be measured directly. Unfortunately, biomechanical optimization problems often possess multiple local minima, making it difficult to find the best solution. Furthermore, convergence in gradient-based algorithms can be affected by scaling to account for design variables with different length scales or units. In this study we evaluate a recently-developed version of the particle swarm optimization (PSO) algorithm to address these problems. The algorithm's global search capabilities were investigated using a suite of difficult analytical test problems, while its scale-independent nature was proven mathematically and verified using a biomechanical test problem. For comparison, all test problems were also solved with three off-the-shelf optimization algorithms--a global genetic algorithm (GA) and multistart gradient-based sequential quadratic programming (SQP) and quasi-Newton (BFGS) algorithms. For the analytical test problems, only the PSO algorithm was successful on the majority of the problems. When compared to previously published results for the same problems, PSO was more robust than a global simulated annealing algorithm but less robust than a different, more complex genetic algorithm. For the biomechanical test problem, only the PSO algorithm was insensitive to design variable scaling, with the GA algorithm being mildly sensitive and the SQP and BFGS algorithms being highly sensitive. The proposed PSO algorithm provides a new off-the-shelf global optimization option for difficult biomechanical problems, especially those utilizing design variables with different length scales or units.     

A comparison of optimisation methods and knee joint degrees of freedom on muscle force predictions during single-leg hop landings

The aim of this paper was to compare the effect of different optimisation methods and different knee joint degrees of freedom (DOF) on muscle force predictions during a single legged hop. Nineteen subjects performed single-legged hopping manoeuvres and subject-specific musculoskeletal models were developed to predict muscle forces during the movement. Muscle forces were predicted using static optimisation (SO) and computed muscle control (CMC) methods using either 1 or 3 DOF knee joint models. All sagittal and transverse plane joint angles calculated using inverse kinematics or CMC in a 1 DOF or 3 DOF knee were well-matched (RMS error<3°). Biarticular muscles (hamstrings, rectus femoris and gastrocnemius) showed more differences in muscle force profiles when comparing between the different muscle prediction approaches where these muscles showed larger time delays for many of the comparisons. The muscle force magnitudes of vasti, gluteus maximus and gluteus medius were not greatly influenced by the choice of muscle force prediction method with low normalised root mean squared errors (<48%) observed in most comparisons. We conclude that SO and CMC can be used to predict lower-limb muscle co-contraction during hopping movements. However, care must be taken in interpreting the magnitude of force predicted in the biarticular muscles and the soleus, especially when using a 1 DOF knee. Despite this limitation, given that SO is a more robust and computationally efficient method for predicting muscle forces than CMC, we suggest that SO can be used in conjunction with musculoskeletal models that have a 1 or 3 DOF knee joint to study the relative differences and the role of muscles during hopping activities in future studies.     

A calibrated EMG-informed neuromusculoskeletal model can appropriately account for muscle co-contraction in the estimation of hip joint contact forces in people with hip osteoarthritis

Abnormal hip joint contact forces (HJCF) are considered a primary mechanical contributor to the progression of hip osteoarthritis (OA). Compared to healthy controls, people with hip OA often present with altered muscle activation patterns and greater muscle co-contraction, both of which can influence HJCF. Neuromusculoskeletal (NMS) modelling is non-invasive approach to estimating HJCF, whereby different neural control solutions can be used to estimate muscle forces. Static optimisation, available within the popular NMS modelling software OpenSim, is a commonly used neural control solution, but may not account for an individual’s unique muscle activation patterns and/or co-contraction that are often evident in pathological population. Alternatively, electromyography (EMG)-assisted neural control solutions, available within CEINMS software, have been shown to account for individual activation patterns in healthy people. Nonetheless, their application in people with hip OA, with conceivably greater levels of co-contraction, is yet to be explored. The aim of this study was to compare HJCF estimations using static optimisation (in OpenSim) and EMG-assisted (in CEINMS) neural control solutions during walking in people with hip OA. EMG-assisted neural control solution was more consistent with both EMG and joint moment data than static optimisation, and also predicted significantly higher HJCF peaks (p < 0.001). The EMG-assisted neural control solution also accounted for more muscle co-contraction than static optimisation (p = 0.03), which probably contributed to these higher HJCF peaks. Findings suggest that the EMG-assisted neural control solution may estimate more physiologically plausible HJCF than static optimisation in a population with high levels of co-contraction, such as hip OA.     

A task-specific validation of homogeneous non-linear optimisation approaches

In biomechanics, musculoskeletal models are typically redundant. This situation is referred to as the distribution problem. Often, static, non-linear optimisation methods of the form “min: φ(f) subject to mechanical and muscular constraints” have been used to extract a unique set of muscle forces. Here, we present a method for validating this class of non-linear optimisation approaches where the homogeneous cost function, φ(f), is used to solve the distribution problem. We show that the predicted muscle forces for different loading conditions are scaled versions of each other if the joint loading conditions are just scaled versions. Therefore, we can calculate the theoretical muscle forces for different experimental conditions based on the measured muscle forces and joint loadings taken from one experimental condition and assuming that all input into the optimisation (e.g., moment arms, muscle attachment sites, size, fibre type distribution) and the optimisation approach are perfectly correct. Thus predictions of muscle force for other experimental conditions are accurate if the optimisation approach is appropriate, independent of the musculoskeletal geometry and other input required for the optimisation procedure. By comparing the muscle forces predicted in this way to the actual muscle forces obtained experimentally, we conclude that convex homogeneous non-linear optimisation approaches cannot predict individual muscle forces properly, as force-sharing among synergistic muscles obtained experimentally are not just scaled versions of joint loading, not even in a first approximation.     

A hybrid static optimisation method to estimate muscle forces during muscle co-activation

The general static optimisation (GSO) process is one of various muscle force estimation methods due to its low computational requirements. However, it can show biased muscle force estimation under muscle co-contraction. In the present study, we introduced a novel hybrid static optimisation (HSO) method to estimate reasonable muscle forces during muscle co-activation movements using more specific equality constraints, i.e. agonist and antagonist muscle moments predicted from a new correlation coefficient approach. The new method was evaluated for heel-rise movements. We found that the proposed method improved the potential of antagonist muscle force estimation in comparison to the GSO solutions. The proposed HSO method could be applied in biomechanics and rehabilitation, for example.     

Coupled objective function to study the role of abdominal muscle forces in lifting using the kinematics-driven model

To circumvent the existing shortcoming of optimisation algorithms in trunk biomechanical models, both agonist and antagonist trunk muscle stresses to different powers are introduced in a novel objective function to evaluate the role of abdominal muscles in trunk stability and spine compression. This coupled objective function is introduced in our kinematics-driven finite element model to estimate muscle forces and to identify the role of abdominal muscles in upright standing while lifting symmetrically a weight at different heights. Predictive equations for the compression and buckling forces are developed. Results are also compared with the conventional objective function that neglects abdominal muscle forces. An overall optimal solution involving smaller spinal compression but higher trunk stability is automatically attained when choosing muscle stress powers at and around 3. Results highlight the internal oblique muscle as the most efficient abdominal muscle during the tasks investigated. The estimated relative forces in abdominal muscles are found to be in good agreement with the respective ratios of recorded electromyography activities.     

A comparative study of muscle force estimates using Huxley's and Hill's muscle model

Determination of muscle forces in individual muscles is often essential to assess optimal performance of human motion. Inverse dynamic methods based on the kinematics of the given motion and on the use of optimisation approach are the most widely used for muscle force estimation. The aim of this study was to estimate how the choice of muscle model influences predicted muscle forces. Huxley's (1957, Prog Biophys Biop Chem. 7: 255–318) and Hill's (1938, Proc R Soc B. 126: 136–195) muscle models were used for determination of muscle forces of two antagonistic muscles of the lower extremity during cycling. Huxley's model is a complex model that couples biochemical and physical processes with the microstructure of the muscle whereas the Hill's model is a phenomenological model. Muscle forces predicted by both models are within the same range. Huxley's model predicts more realistic patterns of muscle activation but it is computationally more demanding. Therefore, if the overall muscle forces are to be assessed, it is reasonable to use a simpler implementation based on Hill's model.     

Hierarchical genetic algorithm versus static optimization-investigation of elbow flexion and extension movements

The applicability of static optimization (and, respectively, frequently used objective functions) for prediction of individual muscle forces for dynamic conditions has often been discussed. Some of the problems are whether time-independent objective functions are suitable, and how to incorporate muscle physiology in models. The present paper deals with a twofold task: (1) implementation of hierarchical genetic algorithm (HGA) based on the properties of the motor units (MUs) twitches, and using multi-objective, time-dependent optimization functions; and (2) comparison of the results of the HGA application with those obtained through static optimization with a criterion "minimum of a weighted sum of the muscle forces raised to the power of n". HGA and its software implementation are presented. The moments of neural stimulation of all MUs are design variables coding the problem in the terms of HGA. The main idea is in using genetic operations to find these moments, so that the sum of MUs twitches satisfies the imposed goals (required joint moments, minimal sum of muscle forces, etc.). Elbow flexion and extension movements with different velocities are considered as proper illustration. It is supposed that they are performed by two extensor muscles and three flexor muscles. The results show that HGA is a suitable means for precise investigation of motor control. Many experimentally observed phenomena (such as antagonistic co-contraction, three-phasic behavior of the muscles during fast movements) can find their explanation by the properties of the MUs twitches. Static optimization is also able to predict three-phasic behavior and could be used as practicable and computationally inexpensive method for total estimation of the muscle forces.     

Constraining spine stability levels in an optimization model leads to the prediction of trunk muscle cocontraction and improved spine compression force estimates

A major limitation of optimization models of the spine has been the inability to accurately predict trunk muscle co-activity. Antagonist muscle activity is thought to be necessary to maintain adequate levels of spine stability but, in turn, creates increased loading on the spine. It is thus hypothesized that the CNS attempts to optimize the relationship between spine loading and spine stability in determining muscular activation patterns. This study presents an optimization model of the spine in which stability was constrained to target levels predicted from regression equations of independent loading variables. Objective functions were set to either minimize the sum of the cubed muscle forces or minimize the sum of the squared intervertebral forces at the L4-L5 disc level. Results demonstrate that the inclusion of stability constraints in optimization simulations produced realistic predictions of antagonist muscle activity and predictions of spine compression levels that agree more closely with EMG-based estimates, compared to simulations in which stability was unconstrained. It was concluded that spinal stability is a vital consideration for the CNS when dictating trunk muscle recruitment patterns.     

Static and dynamic optimization solutions for gait are practically equivalent

The proposition that dynamic optimization provides better estimates of muscle forces during gait than static optimization is examined by comparing a dynamic solution with two static solutions. A 23-degree-of-freedom musculoskeletal model actuated by 54 Hill-type musculotendon units was used to simulate one cycle of normal gait. The dynamic problem was to find the muscle excitations which minimized metabolic energy per unit distance traveled, and which produced a repeatable gait cycle. In the dynamic problem, activation dynamics was described by a first-order differential equation. The joint moments predicted by the dynamic solution were used as input to the static problems. In each static problem, the problem was to find the muscle activations which minimized the sum of muscle activations squared, and which generated the joint moments input from the dynamic solution. In the first static problem, muscles were treated as ideal force generators; in the second, they were constrained by their force-length-velocity properties; and in both, activation dynamics was neglected. In terms of predicted muscle forces and joint contact forces, the dynamic and static solutions were remarkably similar. Also, activation dynamics and the force-length-velocity properties of muscle had little influence on the static solutions. Thus, for normal gait, if one can accurately solve the inverse dynamics problem and if one seeks only to estimate muscle forces, the use of dynamic optimization rather than static optimization is currently not justified. Scenarios in which the use of dynamic optimization is justified are suggested.     

Passive and active mechanical properties of the superficial and deep digital flexor muscles in the forelimbs of anesthetized Thoroughbred horses

The superficial (SDF) and deep digital flexor (DDF) muscles are critical for equine forelimb locomotion. Knowledge of their mechanical properties will enhance our understanding of limb biomechanics. Muscle contractile properties derived from architectural-based algorithms may overestimate real forces and underestimate shortening capacity because of simplistic assumptions regarding muscle architecture. Therefore, passive and active (=total - passive) force-length properties of the SDF and DDF muscles were measured directly in vivo. Muscles from the right forelimbs of four Thoroughbred horses were evaluated during general anesthesia. Limbs were fixed to an external frame with the muscle attached to a linear actuator and load cell. Each muscle was stretched from an unloaded state to a range of prefixed lengths, then stimulated while held at that length. The total force did not exceed 4000 N, the limit for the clamping device. The SDF and DDF muscles produced 716+/-192 and 1577+/-203 N maximum active isometric force (F(max)), had ascending force-length ranges (R(asc)) of 5.1+/-0.2 and 9.1+/-0.4 cm, and had passive stiffnesses of 1186+/-104 and 1132+/-51 N/cm, respectively. The values measured for F(max) were much smaller than predicted based on conservative estimates of muscle specific tension and muscle physiological cross-sectional area. R(asc) were much larger than predicted based on muscle fiber length estimates. These data suggest that accurate prediction of the active mechanical behavior of architecturally complex muscles such as the equine DDF and SDF requires more sophisticated algorithms.     

A physiology-based inverse dynamic analysis of human gait using sequential convex programming: a comparative study

This paper presents an enhanced version of the previously proposed physiological inverse approach (PIA) to calculate musculotendon (MT) forces and evaluates the proposed methodology in a comparative study. PIA combines an inverse dynamic analysis with an optimisation approach that imposes muscle physiology and optimises performance over the entire motion. To solve the resulting large-scale, nonlinear optimisation problem, we neglected muscle fibre contraction speed and an approximate quadratic optimisation problem (PIA-QP) was formulated. Conversely, the enhanced version of PIA proposed in this paper takes into account muscle fibre contraction speed. The optimisation problem is solved using a sequential convex programing procedure (PIA-SCP). The comparative study includes PIA-SCP, PIA-QP and two commonly used approaches from the literature: static optimisation (SO) and computed muscle control (CMC). SO and CMC make simplifying assumptions to limit the computational time. Both methods minimise an instantaneous performance criterion. Furthermore, SO does not impose muscle physiology. All methods are applied to a gait cycle of six control subjects. The relative root mean square error averaged over all subjects, ε(RMS), between the joint torques simulated from the optimised activations and the joint torques obtained from the inverse dynamic analysis was about twice as large for SO (ε(RMS) = 86) as compared with CMC (ε(RMS) = 39) and PIA-SCP (ε(RMS) = 50). ε(RMS) was at least twice as large for PIA-QP (ε(RMS) = 197) than for all other methods. As compared with CMC, muscle activation patterns predicted by PIA-SCP better agree with experimental electromyography (EMG). This study shows that imposing muscle physiology as well as globally optimising performance is important to accurately calculate MT forces underlying gait.     

Effect of optimization criterion on spinal force estimates during asymmetric lifting

Optimization-based muscle force prediction models of the lumbar region are used in research and ergonomic practice, usually as a subroutine of a job analysis software package. Various optimization criteria have been put forward for use in rationally selecting a set of muscle forces to satisfy moment equilibrium, including the sum of cubed muscle contraction intensities and a double linear programming procedure for minimizing the spinal compression force and maximum muscle contraction intensity. A laboratory study was conducted to determine whether these two model formulations produce significantly different estimates of spinal forces for a dynamic asymmetric lift. Although statistically significant differences were found between the predictions of the two models, the difference in peak spinal compression force was only 1.1%.     

Muscle compensation strategies to maintain glenohumeral joint stability with increased rotator cuff tear severity: A simulation study

Rotator cuff tear (RCT) in older adults may cause decreased muscle forces and disrupt the force balance at the glenohumeral joint, compromising joint stability. Our objective was to identify how increased RCT severity affects glenohumeral joint loading and muscle activation patterns using a computational model. Muscle volume measurements were used to scale a nominal upper limb model’s peak isometric muscle forces to represent force-generating characteristics of an average older adult male. Increased RCT severity was represented by systematically decreasing peak isometric muscle forces of supraspinatus, infraspinatus, and subscapularis. Five static postures in both scapular and frontal planes were evaluated. Results revealed that in both scapular and frontal planes, the peak glenohumeral joint contact force magnitude remained relatively consistent across increased RCT severity (average 1.5% and −4.2% change, respectively), and a relative balance of the transverse force couple is maintained even in massive RCT models. Predicted muscle activations of intact muscles, like teres minor, increased (average 5–30% and 4–17% in scapular and frontal planes, respectively) with greater RCT severity. This suggests that the system is prioritizing glenohumeral joint stability, even with severe RCT, and that unaffected muscles play a compensatory role to help stabilize the joint.     

Computational biodynamics of human knee joint in gait: from muscle forces to cartilage stresses

Using a validated finite element model of the intact knee joint we aim to compute muscle forces and joint response in the stance phase of gait. The model is driven by reported in vivo kinematics-kinetics data and ground reaction forces in asymptomatic subjects. Cartilage layers and menisci are simulated as depth-dependent tissues with collagen fibril networks. A simplified model with less refined mesh and isotropic depth-independent cartilage is also considered to investigate the effect of model accuracy on results. Muscle forces and joint detailed response are computed following an iterative procedure yielding results that satisfy kinematics/kinetics constraints while accounting at deformed configurations for muscle forces and passive properties. Predictions confirm that muscle forces and joint response alter substantially during the stance phase and that a simplified joint model may accurately be used to estimate muscle forces but not necessarily contact forces/areas, tissue stresses/strains, and ligament forces. Predictions are in general agreement with results of earlier studies. Performing the analyses at 6 periods from beginning to the end (0%, 5%, 25%, 50%, 75% and 100%), hamstrings forces peaked at 5%, quadriceps forces at 25% whereas gastrocnemius forces at 75%. ACL Force reached its maximum of 343 N at 25% and decreased thereafter. Contact forces reached maximum at 5%, 25% and 75% periods with the medial compartment carrying a major portion of load and experiencing larger relative movements and cartilage strains. Much smaller contact stresses were computed at the patellofemoral joint. This novel iterative kinematics-driven model is promising for the joint analysis in altered conditions.     

Development of a novel MATLAB-based framework for implementing mechanical joint stability constraints within OpenSim musculoskeletal models

The Static Optimization (SO) solver in OpenSim estimates muscle activations and forces that only equilibrate applied moments. In this study, SO was enhanced through an open-access MATLAB interface, where calculated muscle activations can additionally satisfy crucial mechanical stability requirements. This Stability-Constrained SO (SCSO) is applicable to many OpenSim models and can potentially produce more biofidelic results than SO alone, especially when antagonistic muscle co-contraction is required to stabilize body joints. This hypothesis was tested using existing models and experimental data in the literature. Muscle activations were calculated by SO and SCSO for a spine model during two series of static trials (i.e. simulation 1 and 2), and also for a lower limb model (supplementary material 2). In simulation 1, symmetric and asymmetric flexion postures were compared, while in simulation 2, various external load heights were compared, where increases in load height did not change the external lumbar flexion moment, but necessitated higher EMG activations. During the tasks in simulation 1, the predicted muscle activations by SCSO demonstrated less average deviation from the EMG data (6.8% −7.5%) compared to those from SO (10.2%). In simulation 2, SO predicts constant muscle activations and forces, while SCSO predicts increases in the average activations of back and abdominal muscles that better match experimental data. Although the SCSO results are sensitive to some parameters (e.g. musculotendon stiffness), when considering the strategy of the central nervous system in distributing muscle forces and in activating antagonistic muscles, the assigned activations by SCSO are more biofidelic than SO.