Sensitivity Analysis of Gait Models for Realistic Variation of Muscle-tendon Parameters

Team Members

• Katherine Poggensee
• Jon Stingel

Project Video

Project Files

Files for this project can be found on SimTK at https://simtk.org/projects/sensitivity.

Background

Humans are all unique, and as a result, the way in which we move is all different. Going from one place to another we each have the option of how to get there; we can run, we can walk, or we can shuffle or get there using another method. Through all the differences in sizes and speeds of humans and their motion, studies have shown that we may all follow different optimization patterns in the selection of how we actually plan and execute gait (Ackermann, et al., 2010). In order to properly design assistive devices and surgical interventions, an understanding of what these subconscious optimizations are is required. Previous studies by Anderson et al. (2001) and Ackermann et al. (2010) have shown that there are several different types of optimization that we may be doing while in normal gait. One method of optimization looks at minimizing energy consumption while another method is the minimization of muscle fatigue. Ackermann et al. concluded that an optimization of muscle fatigue seemed to produce simulation results that matched typical walking gait kinematics, however previous studies by Anderson and Pandy (2001) showed kinematic simulations that matched normal gait using an energy consumption optimization was used.

In addition to the variety of optimization frameworks, there is significant variability in the individual parameters of the model. There is a significant body of literature quantifying average physiological parameters using both cadaveric and in vivo measurements (Rajagopal, et al., 2016). These studies compile values from different sources, likely resulting in non-negligible measurement variability. Subject-specific differences are also important in individualized modeling. Certain parameters are scaled as a proportion of kinematic measurements from motion capture systems. These values have uncertainty in how the markers are placed as well as in the relative ratios of parameters such as muscle vs. tendon lengths. Other parameters, i.e., the maximum isometric force of a muscle, can be directly measured but are also subject to large degrees of variation. Determining maximum isometric force by maximum voluntary contraction can be done in several different ways for a single muscle, with different results depending on the configuration of the joint and other unmeasurable factors in the nervous system (e.g., motivation, fatigue). Regardless of the mode, these measurements can also be insufficient to normalize electromyography, as certain tasks will result in muscle activations above the "maximum." Some upper-limb muscles can show activations in tasks like boxing at values over double what is seen during the maximum voluntary contraction trials (Halaki & Ginn, 2012).

Because of the large degree of uncertainty in model parameter measurements, there is an underlying need for an analysis on the stability of these optimization methods. Previous studies have observed how well different types of optimization methods will produce movement that is actually seen experimentally. Myers et al. (2014) characterized the propagation in uncertainty in measurement error and parameter estimation. While movement artifact uncertainty propagated universally in all stages of simulation, uncertainty in muscle parameters caused 1.7 times greater impact on the uncertainty in the muscle force output. Muscle force was most sensitive to uncertainty in maximum isometric force and tendon slack length, but this was not universally true across all muscles. Because these parameters are difficult to measure, simulations that are robust to these uncertainties are necessary to produce meaningful results.

The goal of this study is to observe how muscle-tendon parameter uncertainty affects optimization and if those errors fall within experimentally-observed ranges.

Research Question

Does modifying the muscle-tendon parameters of key muscles used in walking alter their activation and the activation of surrounding muscles for fixed kinematics?

Methods

In order to complete the project, we adapted a workflow similar to the one specified in De Groot et et al. (2016), and adapted by our project mentors, Nick Bianco and Chris Dembia. The general process is outlined below, beginning from model adaptation, and ending with the muscle-redundancy solver which provides us with the muscle activations and excitations that we are looking to observe. 

1. Adapting and scaling a model
2. Inverse kinematics and dynamics
3. Modifying muscle-tendon parameters
4. Muscle redundancy solver and direct collocation

Adapting and Scaling a model

In order to apply kinematic data to a model, first, the model must be adapted to match the subject for which the study is being performed. In our study, the model first began as the generalized musculoskeletal muscle model created by Apoorva Rajagopal et al. (2016). This model has 80 muscle-tendon units in the lower limbs, 17 torque actuators, and 37 degrees of freedom. This model, however, was simplified to a model with only 29 degrees of freedom, 17 torque actuators, and 18 muscle-tendon units in the lower limbs. This model was developed by Carmichael Ong et al. (2017) for similar simulation studies. The simpler model still incorporates a level of detail and numbe of muscles so that meaningful results can be obtained, while greatly reducing the complexity from the full model. The difference in models can be seen in Figure 1. 

The simplified model was then scaled in order to match the dimensions of the subject for which kinematic data was captured. The scaling process involves tuning the parameters of the generalized model to match specific heights, weights, and lengths of the subject. In this way, the model becomes customized to the individual, though no scaling is exactly perfect, otherwise our study would not be needed. We would like to thank our mentors (Nick Bianco and Chris Dembia) and members of the teaching team (Apoorva Rajagopal and Carmichael Ong) for performing the necessary model adaptations and scaling. We were lucky to be provided with a model that had already been scaled to the subject data. 

Inverse Kinematics and Dynamics

The next step in the analysis takes the kinematic data that was captured for the subject, and imports it into OpenSim. Once all of the data has been loaded with the musculoskeletal model, inverse kinematics and dynamics can be performed using the respective tools in the OpenSim software. In inverse kinematics, the position marker data captured is used to calculate the positions, velocities, and accelerations of body segments throughout a motion's duration. At this point we not only have positions and orientations of body segments, but we can also determine the lengths and velocities of the individual muscles in the model. This information is very important for later resolving the muscle redundancy problem. Additionally, from inverse kinematics, we know the moment arms of each of the muscles during the prescribed motion. 

Inverse dynamics utilizes the positions, velocities, and accelerations of the joint angles to determine the needed joint-torques that produce the motion. Including the ground reaction forces from the subject data increases the accuracy of the computations. Finally, using the moment arms for each of the muscles and joints, we can determine the resulting forces at the joints. However, since the human body is overactuated, in that there are more muscle actuators than degrees of freedom, the system is overdetermined. Therefore, we must pass the needed torques, moment arms, and muscle properties to a muscle redundancy solver. This will allow us to determine a solution based on an optimization. The whole process that we followed throughout the project is summarized in Figure 2.

Modifying muscle-tendon parameters

To perform the sensitivity analysis, we leveraged the data gathered from several sources and reported in Rajagopal et al. (2016) in order to establish a range to which muscle-tendon parameters could be changed. Values listed are the baseline parameter value with one standard deviation in parentheses. This data was gathered using various methods, including cadaver studies.

In order to perform our sensitivity analysis, we needed to vary parameters in the model by some regulated amount. For our study, we focused on the parameters that were consolidated in Rajagopal et al. 2016:

• Maximum isometric force
• Optimal fiber length
• Tendon slack length
• Pennation angle

We varied each of these parameters in the model by plus or minus one standard deviation. Parameters were varied individually three muscles, the medial gastrocnemius, soleus, and tibialis anterior. Shown below is an example section from a model file, highlighting the parameters that we altered throughout the project. Once this model was altered, we solved for activations using the muscle redundancy solver. In this way, we performed a fixed kinematic and dynamic optimization, but introduced variation in the muscle-tendon parameters individually to observe the changes in the activations just as a result of the optimization. 

Example model construction for one muscle (medial gastrocnemius):

<Millard2012EquilibriumMuscle name="med_gas_r">
<!--The set of points defining the path of the muscle.-->
<GeometryPath> ... </GeometryPath>
<!--Maximum isometric force that the fibers can generate-->
<max_isometric_force>4935.74134</max_isometric_force>
<!--Optimal length of the muscle fibers-->
<optimal_fiber_length>0.0551036008753214</optimal_fiber_length>
<!--Resting length of the tendon-->
<tendon_slack_length>0.4307981501</tendon_slack_length>
<!--Angle between tendon and fibers at optimal fiber length expressed in radians-->
<pennation_angle_at_optimal>0.16568155</pennation_angle_at_optimal>
<!--Compute muscle dynamics ignoring tendon compliance. Tendon is assumed to be rigid.-->
<ignore_tendon_compliance>false</ignore_tendon_compliance>
<!--The linear damping of the fiber.-->
<fiber_damping>0.01</fiber_damping>
<!--Assumed initial activation level if none is assigned.-->
<default_activation>0.01</default_activation>
<!--Activation lower bound.-->
<minimum_activation>0.01</minimum_activation>
<!--Active-force-length curve.--> ... </ActiveForceLengthCurve>
<!--Force-velocity curve.-->
<ForceVelocityCurve> ... </ForceVelocityCurve>
<!--Passive-force-length curve.-->
<FiberForceLengthCurve> ... </FiberForceLengthCurve>
<!--Tendon-force-length curve.-->
<TendonForceLengthCurve> ... </TendonForceLengthCurve>
</Millard2012EquilibriumMuscle>

Muscle redundancy solver and direct collocation

For each of the different variations of the project, the model, parameters, kinematics, and dynamics were passed into a direct collocation pipeline. These methods are based on the direct collocation optimal control framework for solving the muscle redundancy published by Friedl De Groote (2016). In this method, the optimizer minimizes a specified objective function, and computes the muscle activations and excitations for that solution. For our project, we focused on a single objective function of muscle activations squared, as shown below, in which w1 is a heavy weight applied to the activations term. 

This objective function is minimized while meeting the constraints for the problem, including the dynamics for the system which are passed as constraints with a heavy weight. We initially selected activations squared as the objective since it is an accepted function of how humans naturally optimize during tasks such as walking. 

Limitations

Though our project has taken what we believe is a good step in determining just how sensitive simulations are to the underlying muscle tendon parameters that we specify in our models, there are certain limitations. First, we would like to address that the muscle-tendon parameters from person to person vary greatly. As a result, no model is going to fit any one person perfectly. The generalized parameters that were used to create the model come from an average of values obtained in a combination of different studies using various measurement techniques. This combination can introduce greater variability in the values that are obtained for specific parameters. Additionally, tendon slack lengths and optimal fiber forces can't be directly measured in vivo, and the approximations that we have for them could then carry error across all subjects. Therefore, the amount by which we vary each of the muscle-tendon parameters may not be the best representation of how much those parameters can vary in humans. It can be extremely difficult to create a well-scaled model of a subject. Motion capture markers can move, creating distortions in the measurements, and each individual is shaped differently; there is not a perfect way to scale the model to match each dimension and parameter. As a result, the model that we used will not be a perfect representative of the actual human that data was captured for. 

Finally, we recognize that our study is limited in its processing of the muscle activations. Though we have found a qualitative difference in the activations patterns of the muscles, a detailed statistical analysis is needed in order to quantify just how different the activations are for each of the different sets of conditions. Muscle-tendon interactions, as well as the ways in which they work together to create movement, is complicated, and the effect of changing different parameters is not always clear. This analysis may lead to a better understanding of what parameters are the most sensitive and could aid in the development of more accurate models. 

Future Work

In working on this project, we needed to compile a large number of resources in the form of model data, software packages, and scripts developed by a number of different people. This became somewhat of a challenge, as it became difficult to navigate the structure of several different code-bases, as well as ensure that the adapted version of these different methods was working properly for the model and data that we were using. However, with all of these tools and resources available to us, there are many possibilities for different aspects of human movement to study. As a result, it can be difficult to pinpoint a specific research question.

Once we identified the sensitivity analysis of muscle-tendon parameters on muscle activation optimization as the research topic, there were several avenues in which we would have liked to travel down further. In our study, we varied four parameters in three separate muscles individually. In the future, it would be useful to not only expand to other muscles but to examine the effects of multiple parameters changed within the model, as there are likely interdependencies. This could give more insight as to the ways in which muscles coordinate to perform tasks. Similarly, performing these analyses on different kinds of motion or different tasks would also be a great direction for future work, as the patterns and ranges of activations are task dependent.

Similarly, in the future, this work should be expanded to different objective functions as well. In our exploration, we used activation squared as the optimization function, as minimizing the sum of muscle activations squared has been found in previous studies (Ackermann et al. 2010) to produce motion very similar to the ways in which our bodies naturally select. It would be interesting through to expand to objectives such as the minimization of metabolic cost over the course of the motion. The sensitivity to these parameters with respect to different objectives is one question that we originally wanted to explore, though were unable to address in this study.

All of these directions for future work will help to build and expand the understanding of how sensitive the result of complex optimizations are as a result of the underlying muscle-tendon parameters that are used in model creation. Our work has taken a good first step in understanding what parameters are the most influential to the resulting muscle-activations obtained from an optimization using direct collocation methods. Understanding these sensitivities allows us to know if the results we obtain are accurate, and also indicates which parameters can be the most influential as to why a movement is happening the way it does in the body. We can then create more accurate and detailed models, and form a better understanding biomechanics as a whole. 

Acknowledgments

We would like to thank the ME485 teaching team, class, and members of the Neuromuscular Biomechanics Lab at Stanford for their assistance and guidance throughout this project, as well as for their aid in providing subject data, and helpful scripts needed to complete the project tasks. 

References

Ackermann, M., Bogert, A.J.V.D. (2010). Optimality principles for model-based prediction of human gait. Journal of Biomechanics 43, 1055–1060. doi:10.1016/j.jbiomech.2009.12.012

Alexander, R.McN., Bennet-Clark, H.C. "Storage of elastic strain energy in muscle and other tissues." Nature 265.5590 (1977): 114. doi:10.1038/265114a0.

Anderson, F.C., Pandy, M.G. (2001). Dynamic Optimization of Human Walking. Journal of Biomechanical Engineering 123, 381. doi:10.1115/1.1392310

De Groote, F., Kinney, A. L., Rao, A. V., & Fregly, B. J. (2016). Evaluation of Direct Collocation Optimal Control Problem Formulations for Solving the Muscle Redundancy Problem. Annals of Biomedical Engineering, 44(10), 2922–2936. doi:10.1007/s10439-016-1591-9

Halaki, M., Ginn, K. Normalization of EMG signals: To normalize or not to normalize and what to normalize to?. Computational intelligence in electromyography analysis-a perspective on current applications and future challenges. InTech, 2012. doi: 10.5772/49957

Millard, M., Uchida, T., Seth, A., & Delp, S. L. (2013). Flexing Computational Muscle: Modeling and Simulation of Musculotendon Dynamics. Journal of Biomechanical Engineering, 135(2), 021005. doi:10.1115/1.4023390

Myers, C. A., Laz, P. J., Shelburne, K. B., & Davidson, B. S. (2015). A probabilistic approach to quantify the impact of uncertainty propagation in musculoskeletal simulations. Annals of biomedical engineering, 43(5), 1098-1111. doi:10.1007/s10439-014-1181-7

Ong,C. F., Geijtenbeek T., Hicks J. L., & Delp S. L. (2017). Predictive Simulations of Human Walking Produce Realistic Cost of Transport at a Range of Speeds. XVI International Symposium on Computer Simulation in Biomechanics July 20th – 22nd 2017. Gold Coast, Australia.

Rajagopal, A., Dembia, C. L., DeMers, M. S., Delp, D. D., Hicks, J. L., & Delp, S. L. (2016). Full-Body Musculoskeletal Model for Muscle-Driven Simulation of Human Gait. IEEE Transactions on Biomedical Engineering, 63(10), 2068–2079. doi:10.1109/TBME.2016.2586891