1 EMG BASED MUSCLE FORCE ESTIMATION By JONATHAN PAUL WALTER A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNI VERSITY OF FLORIDA 2012
2 2012 Jonathan P. Walter
3 T o my colleagues and future bio mechanics m ay this k nowledge help y ou achieve a g reater s tate of e nlightenment and fore sight.
4 ACKNO WLEDGMENTS Thank you to my family, friends, and colleagues for the support they provided t hroughout my life I am especially thankful for the support and love provided by my wife, Dr. Nicole Kahhan. Thanks to all those people who have helped Nicole and I accompl ish our goals together and for h ousing us in times of need. Acknowled gments should be provided to the Computational Biomechanics Laboratory and Orthopaedic Biomechanics Laboratory for their help and contributions to this project I am especially thankful to Dr. Allison Kinney and to Dr. Benjamin Fregly for their g uidance and technical contributions T hank s to my dissertation committee members especially Dr. Scott Banks, for their ass istance and approval of this work Thanks to the m any collaborators w hom paved the way for this project and were fundamental in collecting the special data required for this proj ect. A cknowledgements should also be provided to the research funding programs o f NSF, NI H t he Shiley Center for Orthopaedic Research, and the University of Florida Alumni Fellowship program for their financial support o f the current research
5 TABLE OF CONTENTS page ACKNOWLEDGMENTS ................................ ................................ ................................ .. 4 LIST OF TABLES ................................ ................................ ................................ ............ 7 LIST OF FIGURES ................................ ................................ ................................ .......... 8 LIST OF ABBREVIATIONS ................................ ................................ ........................... 10 ABSTRACT ................................ ................................ ................................ ................... 13 CHAPTER 1 INTRODUCTION ................................ ................................ ................................ .... 15 2 EVALUATION OF CONTACT FORCES WITH MEDIAL COMPARTMENT KNEE OSTEOARTHRITIS ................................ ................................ ...................... 17 Background ................................ ................................ ................................ ............. 17 Methods ................................ ................................ ................................ .................. 18 Results ................................ ................................ ................................ .................... 22 Discussion ................................ ................................ ................................ .............. 23 3 EVALUATION OF MUSCULOTENDON MODELS WITH EMG BASED MUSCLE FORCE ESTIMATION ................................ ................................ ............ 33 Background ................................ ................................ ................................ ............. 33 Methods ................................ ................................ ................................ .................. 35 Dynamic Musculoskeletal Model Creation ................................ .............................. 36 Kinematics ................................ ................................ ................................ .............. 37 Joint Loads ................................ ................................ ................................ ............. 38 EMG to Force Model ................................ ................................ .............................. 40 EMG Based Optimizations ................................ ................................ ...................... 43 Initial Guess ................................ ................................ ................................ ............ 47 Parameter Calibration ................................ ................................ ............................. 47 Model Testing ................................ ................................ ................................ ......... 48 Results ................................ ................................ ................................ .................... 49 Discussion ................................ ................................ ................................ .............. 49 4 EVALUATION OF NEURAL SYNERGIES FOR MUSCLE FORCE ESTIMATIONS ................................ ................................ ................................ ....... 62 Background ................................ ................................ ................................ ............. 62 Methods ................................ ................................ ................................ .................. 64 Dynamic Musculoskeletal Model Creation ................................ .............................. 65
6 Kinematics ................................ ................................ ................................ .............. 66 Joint Loads ................................ ................................ ................................ ............. 67 Neural Synergy Factorization ................................ ................................ .................. 68 EMG to Force Model ................................ ................................ .............................. 68 EMG Based Optimizations ................................ ................................ ...................... 70 Initial Guess ................................ ................................ ................................ ............ 74 Parameter Calibration ................................ ................................ ............................. 74 Model Testing ................................ ................................ ................................ ......... 75 Results ................................ ................................ ................................ .................... 76 Discussion ................................ ................................ ................................ .............. 77 5 CONCLUSIONS ................................ ................................ ................................ ..... 94 APPENDIX A DISCRETIZATION OF ACTIVATION AND CONTRACTION DYNAMICS .............. 96 Background ................................ ................................ ................................ ............. 96 Methods ................................ ................................ ................................ .................. 97 Activation Dynamics ................................ ................................ ................................ 97 Activation Dynamics Numerical Solution ................................ ................................ 98 Contraction Dynamics ................................ ................................ ............................. 99 Contraction Dynamics Numerical Solution ................................ ............................ 107 Results and Discussion ................................ ................................ ......................... 110 B EMG TO FORCE MODELLING ................................ ................................ ............ 111 Motivation ................................ ................................ ................................ ............. 111 Optimization and Calibration ................................ ................................ ................. 112 LIST OF REFERENCES ................................ ................................ ............................. 116 BIOGRAPHICAL SKETCH ................................ ................................ .......................... 122
7 LIST OF TABLE S Table page 2 1 Percent reductions in peak and impulse values relative to normal gait. ............. 26 2 2 Prediction R 2 values between tested variables. ................................ .................. 27 3 1 Literature of EMG based models for any joint and any task. .............................. 53 3 2 Literature de scribing EMG based models for the full leg during gait. ................. 54 3 3 Net Load R 2 values for test trial. ................................ ................................ ......... 55 3 4 EMG R 2 values for test trial. ................................ ................................ .............. 56
8 LIST OF FIGURES Figure page 2 1 Trial means of primary variables. ................................ ................................ ........ 28 2 2 Trial m eans of secondary variables. ................................ ................................ ... 29 2 3 Changes in primary peak and impulse values. ................................ ................... 30 2 4 Changes in secondary peak values. ................................ ................................ ... 31 2 5 R 2 value sensitivity to rotation of knee adduction axes. ................................ ...... 32 3 1 Exponent selection method for experimental tracking errors. ............................. 57 3 2 Exponent selection method for boundary errors. ................................ ................ 57 3 3 Average Net Load matching R 2 values for calibration and prediction.. ............... 58 3 4 Average EMG matching R 2 values for calibration and test trials. ........................ 59 3 5 Average muscle forces differences to full model for calibration trial. .................. 60 3 6 Comparison of muscle forces between models for calibration trial. .................... 61 3 7 Optimized knee anterior posterior translation. ................................ .................... 61 4 1 Net load tracking quantities. ................................ ................................ ............... 83 4 2 Net loads modeled with desired/actual net load for Match cases. ...................... 84 4 3 Net loads modeled with desired/actual net load for PreCalPred cases. ............. 85 4 4 Net loads modeled with desired/actual net loads for Pred cases. ....................... 86 4 5 EMG tracking quantities at optimal time shifts. ................................ ................... 87 4 6 EMG tracking normalized RMS errors at optimal time shifts. .............................. 88 4 7 EMG tracking R 2 values at optimal time shifts. ................................ ................... 89 4 8 Contact force quantities for medial, lateral and total forces. ............................... 90 4 9 Contact forces modeled and desired/actual for Match cases. ............................ 91 4 10 Contact forces modeled and desired/actual for PreCalPred cases. .................... 91 4 11 Contact forces modeled and desired/actual for Pred cases. .............................. 92
9 4 12 Muscle forces with large differences and large contributions to contact forces for Match cases. ................................ ................................ ................................ 92 4 13 Muscle forces with large differences and large contributions to contact forces for PreCalPred cases. ................................ ................................ ........................ 93 4 14 Muscle forces with lar ge differences and large contributions to contact forces for Pred cases. ................................ ................................ ................................ ... 93
10 LIST OF ABBREVIATION S A DD B REV Adductor brevis A DD L ONG Adductor longus A DD M AG Adductor magnus AA Adduction Abduction AP Anterior Posterior body centered direction B IFEM SH Short head biceps femoris B IFEM LH Long head biceps femoris B SP B spline parameterization C AL Calibration of model parameters C T Computerized t omography D ELAY O NLY P ure electro mechanical time delay with a pure force generator mus cle model D OF Degrees of f reedom EDL Extensor digitorum longus EHL Extensor hallucis longus EMG Electromyography FDL Flexor digitorum longus FE Flexion extension axis FHL Flexor hallucis longus F ULL P ure electrical time delay with activation dynamics and compliant tendon contraction dynamics model G EM Gemellus G L M AX Gluteus maximus G L M ED Gluteus medius G L M IN Gluteus minimus
11 G RAC Gracilis GRF Ground reaction forces (and moments) I LIAC Iliacus I MP Impulse IE Internal external K AM Knee a dduction m om ent K FM Knee f lexion m oment L AT G AS Lateral gastrocnemius M ATCH Matching of contact forces with parameter calibration M ED G AS Medial gastrocnemius ML Medial lateral M RI Magnetic r esonance i mage M CF Medial c ontact f orce N Newt ons N O A CT D YNAMS P ure electrical t ime delay with complia nt tendon contraction dynamics model ODE Ordinary differential equation P ECT Pectineus P ER B REV Peroneus brevis P ER L ONG Peroneus longus P ER T ERT Peroneus tertius P IRI Piriformis P K 1 Peak 1 corresponding to ~25% stance phase P K 2 Pea k 1 corresponding to ~75% stance phase P RE C AL P RED Prediction of contact forces with pre calibrated parameters P RED Prediction of contact forces with calibration of parameters
12 P SO Psoas P URE F ORCES P ure electrical time delay with activation dynamics and a pure force generator muscle model Q UAD F EM Quadratus femoris RF Rectus femoris RMS Root mean square S ART Sartorius S EMI M EM Semimembranosus S EMI T EN Semitendinosus SI Superior inferior S OL Soleus S TIFF T ENDON P ure electrical time delay with activation dyn amics with stiff tendon contraction dynamics model S YN Neural synergy parameterization T CS Tibial f ixed c oordinate s ystem T IB A NT Tibialis anterior T IB P OST Tibialis posterior TFL Tensor fascia lata VAF Variance accounted for V AS I NT Vastus intermedius V A S L AT Vastus lateralis V AS M ED Vastus medialis
13 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy EMG BASED MUSCLE FORCE ESTIMATION By Jonathan P. Walter December 2012 Chair: Benjamin Fregly Major: Mechanical Engineering E xcessive contact forces in the medial compartment of the knee during walking are believed to progress articular cartilage damage in o steoarthrit is A bn ormally large adduction moments of the knee from gait analyses have been linked to p rogression and severity of this disease The results of the first study presented here indicate that reducing the knee adduction moment does not necessarily guarantee a red uction in corresponding contact forces Therefore, a more accurate measure ment of the loads on the knee, especially the balance between contact and muscle forces, is desired for a subject specific biomechanical analysis of subjects with neuromusculoskeleta l disorders. However, t he inability to easily measure muscle forces in vivo along with the inability to calculate unique muscle forces directly has greatly hindered achievement of this goal. The use of electromyography (EMG) to constrain shapes of muscle excitations has been helpful to produce muscle force estimates but require EMG to force models with unknown model parameters and EMG collection is limited to a subset of leg muscles The concept of neural synergies is developed here to parameterize the m uscle excitations of all muscles based on a small set of time varying neural commands. The second and third studies presented here investigate d how the complexity of EMG to
14 force models and the inclusion of EMG tracking or neural synergy parameterization affect estimates of muscle and contact forces in the leg during walking using optimization methods The investigation s were performed using experimental motion capture, ground reaction, EMG and knee contact loads collected from a subject implanted with a force measuring tibial prosthesis. A small trend for better contact force predictions was found with the use of synergies to constrain muscle excitations and also with the removal of EMG tracking. It is believed that EMG based constraint me thods suffer fr om the difficulties of accurate EMG collection model parameter uncertainty, and the lack of experimental excitation magnitude constraints. Further research should be performed to determine how well the complexity of EMG to force models affects contact for ce estimates and to determine the best model parameter calibration techniques for EMG based muscle force estimation s
15 CHAPTER 1 INTRODUCTION The ability to walk is one of the most important tasks that humans perform. Ambulation is commonly used as an indic ator for quality of life in a wide range of clinical studies. Problems that arise from neurological and biomechanical disorders in the human body can lead to pain and the eventual inability to walk. Biomechanical evaluation and treatment of these disorders on a subject specific level are essential to the advancement of the knowledge concerning these disorders Treatments for neuro biomechanical disorders can range from physical therapy sessions that are intended to train muscles and motions, to surgical int ervention s that are intended to correct major deformities. Treatment outcomes are tested by study ing walking patterns range of motions, and skeletal geometry ( using imaging techniques such as magnetic resonance imaging ( MRI ) and computerized tomography ( C T ) ) For example, a subject diagnosed with osteoarthritis of the knee who suffers from serious knee pain during walking, may undergo surgery to replace the articulating bone surfaces with implants. S urgeons and physical therapists can then evaluate the bio mechanical loading conditions that are present in the knee before an d after surgery. Typical analysi s in a gait laboratory includes measurement of marker location s attached to body segment s and the magnitude of loads applied to the ground by the feet Pati ent specific modeling and a dynamic analys is of these measurements enable clinicians to evaluate the external or overall loading conditions at any joint of the body. For a more thorough analysis it is desirable to ascertain the magnitude of the internal l oads at the knee, including loads from muscles applied to the bone and contact loads between articulating surfaces. However, due to the inability to directly measure the
16 muscle and contact forces without extraneous strain measuring implant devices these m easurements are not commonly available. A biomechanical analysis provides insight into the timing and magnitudes of these loads by utilizing musculoskeletal models of the subject along with a combination of complex numerical analysis techniques and dynamic analyses. E xcessive contact forces of the medial compartment of the knee during walking are believed to progress articular cartilage damage in o steoarthrit is A bnormally large adduction abduction ( AA ) moments of the knee from gait analyses have been linke d to p rogression and severity of this disease An accurate estimation of knee contact forces would be a powerful aid for osteoarthritis treatment. The first study presented here examined the relationship between contact forces and knee AA moments (Walter e t al., 2010) to determine how well the AA moment could predict changes in contact fo rces. The second and third studies were performed to examine the ability of novel muscle forces methods to predict unique muscle and knee contact force estimates. Since mu scle forces cannot be easily measured or calculated, numerical optimization methods and electromyography (EMG) have been used to constrain patterns and magnitudes of muscle excitations in order to produce muscle force estimates. The use of neural synergies was developed as a novel method to constrain the muscle activity of all muscles based on a small set of time varying neural commands and a small set of synergistic pathways. T he complexity of EMG to force models as well the inclusion of EMG tracking and n eural synergy constraints were varied to determine their e ffect on estimated muscle forces in the leg and predicted contact forces in the knee during walking.
17 CHAPTER 2 EVALUATION OF CONTAC T FORCES WITH MEDIAL COMPARTMENT KNEE OSTEOARTHRITIS Background Med ial compa rtment knee osteoarthritis plagues a lar ge portion of the population ( Helmick et al., 2008 ) The disease process is often initiated by changes in motions and loads at the knee joint such as those following anterior cruciate ligament or meniscal i njury. Since overloading of the medial compartment in particular has been hypothesized to play a critical role, many treatment approaches have made reduction of medial contact force ( MCF ) during gait their primary target. Unfortunately, no accurate method currently exists for determining in vivo medial and lateral knee contact forces in individual patients. External measurements and dynamic analyses can determine the net forces and moments experienced by a joint, but knowledge of muscle and ligament forces is needed to determine contact forces accurately. Since contact force predictions from musculoskeletal computer models have yet to be thoroughly validated, researchers have proposed the knee adduction moment ( KAM ) as a surrogate measure for medial compartm ent load ( Andriacchi 1994 ) I t has be e n observed that patients with high AA moments require increased antagonistic muscle loads to maintain the balance between the medial and lateral compartments ( Schipplein and Andriacchi 1991 ) This moment normally exh ibits two peaks during stance phase, with the first one typically being the largest. The value of the first peak during stance phase has been shown to be highly correlated with disease severity ( Sharma et al., 1998 ) and the r ate of disease progression ( Miy azaki et al., 2002 ) The angular impulse of the AA moment has been shown to be well correlated with radiographic based disease severity for mild and moderate patients ( Thorp et al., 2006 )
18 as well as symptomatic knee pain in mild patients ( Thorp et al., 200 7 ). T he instantaneous value of this moment during stance phase has been shown to be highly correlated with MCF measurements provided by an instrumented knee implant ( Zhao et al., 2007 ) However, no study has shown that changes in the peak KAM are an accura te indicator of changes in peak MCF during stance phase. Proper expression of the KAM on an anterior posterior ( AP ) axis is critical to determine the load split between medial and lateral compartments. Expressing the moment on a different axis can result in significant differences of the corresponding moment ( Schache et al., 2008 ) This can lead to varying conclusions about how the moment changes between different gait patterns as well as how this quantity is related to other kinetic and kinematic quantiti es of gait. An improperly aligned AP axis could lead to a misrepresentation of the load balance between the medial and lateral knee compartments. The current study uses gait data collected from a subject with a force measuring knee replacement to determin e whether changes in both peaks of the KAM curve are accurate indicators of corresponding changes in both peaks of the MCF curve. In order to account for internal external ( IE ) rotation of the shank, the AP axis was aligned with a line connecting the heel and toe markers. The results challenge the current hypothesis that reduction of the largest peak in the KAM curve necessarily means a reduction in the largest peak of the MCF curve. Methods One patient with a force measuring knee replacement (male, right k nee, age 83, mass 68 kg, height 1.7 m, body mass index 23.5, neutral leg alignment) performed overground gait with simultaneous collection of internal knee contact force, external
19 video motion, and external ground reaction data ( Fregly et al., 2009 ) Insti tutional review board approval and informed consent were obtained. The implant utilized a custom tibial prosthesis instrumented with four uniaxial force transducers, a micro transmitter, and an antenna, where the transducers measured compressive force at t he four corners of the tibial tray ( D'Lima et al., 2005 ) The subject performed five trials each for three different pole gait. These gait patterns have been shown t o reduce both peaks of the knee AA ( Fregly et al. 2007 ) The subject performed all three gait patterns at his self selected walking speed of 1.23 m/s. The subject also performed joint range of motion trials at the ankle, hip, and knee along with a static standing pose. Using data from ground reaction load s and three dimensional video marker motion locations, the i nverse dynamic knee loads were calculated using a patient specific dynamic gait model. The equations of motion were derived with Autolev using K method. A full body model was used for parameter calibration which consisted of 27 degrees of freedom ( DOF ) Six DOFs were used to place and orient the pelvis segment relative to the fixed laboratory coordinate system. The remaining DOFs included thr ee rotations at the hip, one at the knee, two at the ankle, three at the back, two at the shoulder, and one at the elbow. Nominal model parameters were found from mo tion capture using the standing pose. Segment masses, mass centers, and moments of inertia were estimated using regression relationships They were then calibrated with gait capture data to reduce any residual loads on the pelvis via non linear least squares optimization techniques. All o rientations of joint axes and locations of joint centers i n segmental frames were calibrated by optimization to minimize the error between
20 simulated and experimental marker positions during joint range of motion trials and one gait trial The nominal knee flexion extension (FE) axis was defined as the line conne cting the femoral epicondyle markers. This functional axis was used as a common axis between the thigh and shank segments. The knee FE angle was therefore calculated about this axis. The knee joint center was defined to be equidistant from knee epicondyle markers. The knee joint center served as the origin of a tibial fixed coordinate system ( T CS) system The superior inferior (SI) axis was defined as the line between the knee and ankle joint centers. The medial lateral ( ML ) axis was defined as the cross pr oduct of the proximal/distal with a line connecting the heel and index toe markers. The AP axis was then defined as the cross product of the other two axes and therefore remained in a plane connecting the shank SI axis and the foot AP axis. This allowed an AP axis which would account for IE rotation of the shank. A separate shank foot model was created using the calibrated full body gait model's parameters. This allowed the lower leg to be more closely aligned to the experimental data and therefore provide d better estimates of joint moments and forces. The full body model's knee center and TCS served as the first 6 DOFs for both sides. This re aligned TCS was then used to define the f lexion/ e xtension and AA joint moments at the knee. The KAM was defined as positive when it is a moment which would adduct the lower leg, and similarly with the knee extension moment. The flexion/extension moment was also calculated about the functional FE axis (the FE axis which the knee was rotated about for kinematics) for fur ther analyses.
21 The MCF was determined from manipulation of data from four load cells of an instrumented knee implant. A previously validated regression equation for this patient's implant was used to determine the load split between the medial and lateral femoral condyles ( Zhao et al., 2007 ) The peak values of KAM and MCF during stance phase were calculated for each gait pattern. This was done by finding local maxima that were closest to 25% and 75% of the stance phase for Peak 1 and Peak 2 respectively. Impulse values for the MCF and the KAM were calculated by integrating the curve with respect to time throughout stance phase. Peak values from the knee flexion moment, knee flexion angle, and hip internal rotation angle were also calculated since they are known to change significantly for medial thrust gait (Fregly et al., 2007) Since no consistent peaks could be found for the hip internal rotation angle, a window average between 15 25% and 65 75% stance phase were used for surrogate Peak 1 and Peak 2 meas urements, respectively. Changes in peak values for medial thrust and walking pole gait relative to normal gait were analyzed statistically using a two tailed Mann Whitney U test, which is the non parametric equivalent to a two tailed t test. The level of s tatistical significance was set at p 0.05. The ability of the KAM measures to predict the measures of the MCF was investigated using linear regression analysis and the corresponding R 2 values. In order to determine the sensitivity of results to the dire ction of the AP tibial axis in which the KAM was expressed th is axis was rotated about the shank IE rotation axis. This rotation was done at 1 degree increments from 25 to 25 degrees. This rotation represent s the addition of varying amounts of knee flexi on moment into the KAM and acts as a sensitivity study of errors in the AP axis direction.
22 Results The peak KAM was substantially reduced in both modified gait patterns (Fig ure 2 1). These changes did not correspond to reductions in peak MCF (Table 2 1). During early stance phase, both modified gait patterns had statistically significant reductions in the pea k of the KAM while n either gait pattern had a statistically significant change in the peak of the MCF (Fig ure 2 2 ). In contrast, during the late r part of stance phase, only the walking pole gait had a statistically significant reduction in peak KAM while both gait patterns had statistically significant reductions in MCF (Figure 2 3 ). The impulse and angular impulse all showed significant reductions for each of the modified gait patterns. Changes consistent with those desired for medial thrust were found in the knee flexion angle, knee extension moment, and hip internal rotation angle (Fig ure 2 2 ). First peak changes of the extension moment and flexion angle were consistent with the first peak changes of the KAM but not the MCF (Figure s 2 3, 2 4). Throughout stance phase the knee FE moment of both walking pole and medial thrust increased. The knee flexion angle of walking pole saw an increase in the firs t half of stance phase similar to medial thrust but remained close to normal in the second half. The hip IE rotation angle remained close to normal except for one significant difference during the first half of stance of medial thrust gait. Rotation of th e AP axis created large variations in prediction ability (Figure 2 5). At zero rotation, r egressi o n analysis revealed that the first peak of the KAM curve was a poor predictor of the first peak in MCF (Table 2 2 ). T he second peak to second peak ( R 2 = 0 .68) and angular impulse to impulse ( R 2 = 0 .56) predictions both had relatively high R 2 values. The first peak to first peak predictions showed dramatic improvement at an internal rotation of the AA axis of 20 degrees. However, t his rotation was found to have
23 the lowest R 2 value when u sing the first peak of the KAM to predict the second peak and the impulse of the MCF. No amount of rotation allowed the first peak of the MCF to be well predicted by either the second peak or impulse of the KAM. A rotation of 4 d egrees was an optimal rotation to predict the second peak of the MCF. Optimal rotations using the angular impulse occurred at 9 degrees for predicting the second peak MCF ( R 2 = 0 .60) and 15 degrees for predicting the impulse of the MCF ( R 2 = 0 .74). The lar gest prediction R 2 value of all combinations occurred when using the second peak of the KAM to predict the impulse of the MCF at a rotation of 5 degrees ( R 2 = 0 .79). Discussion This study evaluated the relationship between changes in both KAM measures and i n the corresponding MCF measures The availability of instrumented knee implant data collected simultaneously with the gait data provided a unique opportunity to investigate this relationship. The results indicate that reducing the peak KAM does not necess arily guarantee a reduction in corresponding peak MCF Because of the lack of expected changes in the first peak of the MCF in this patient, other factors must play an important role in determining changes produced by a clinical intervention such as gait modification. Identification of these other factors is difficult given the current set of gait data used in this study. Since the first peak of the MCF curve was well fitted by a constant coefficient in the linear regression analysis, adding other covariat es to the regression analysis produced little improvement in the fit. Looking at important changes between the walking patterns, i t seems possible that the increased knee flexion during the first half of stance phase of medial thrust and walking pole gait necessitated increased peak knee extension moment. This moment increase may have been the result of increased force generation by knee extensor muscles,
24 which would also ten d to increase joint compression. An increase in joint compression would offset the reduction produced by knee medialization or load transfer through the walking poles. Even when walking poles were used, which reduce d the ground reactions applied at the feet, the first peak of the total contact force in the knee remained statistically the same as those from normal gait. This further implies that an increase in knee extension moment must have increased joint compression. In order to ensure that changes of the knee extension were accurate, this moment was also calculated about the functional FE axis which resulted in similar statistical differences of peak values between gait trials Despite the poor predictions found using the first peak for this patient, both the second peak and impulse variations between gait patterns seemed to be well p redicted. Angular impulse has been shown in the past to be well tied to clinical measures ( Thorp et al., 2006 ) and was found in this study to be a good predictor of medial force impulse. Previou s studies have demonstrated high correlation of the peak KAM to clinical measures as well ( Sharma et al., 1998 ; Miyazaki et al., 2002 ) While this patient was performing modified gait patterns, abnormal muscle activations may have occurred. Someone performing their natural walking motion at different times during os teoarthritis progression could have more normal muscle activations and the KAM changes could be more predictive of medial force changes. For example, if no gait retraining treatments are performed, subjects may not increase the knee FE angle of moment and thus KAM changes may be more significantly correlated to changes in MCF. Considering this, changes in the first peak of the KAM should still be an important variable to consider for gait analysis.
25 The modification of the AP axis created significant differe nces in the KAM and the predictions made from it. Internal rotation of the AP tibial axis by 20 degrees resulted in better prediction of the MCF peaks. This occurs due to a combination of the non rotated AA and FE moments that creates a lack of first peak changes between trial types of the KAM. The differences in the second peak between gait types were also sensitive to the rotation of the AP axis. These differences in prediction abilities support the idea that expression of the KAM is important to consider for gait analysis. While these findings show limitations of the ability of the KAM to predict changes in the MCF they are limited to only one patient studied. This patient has significantly altered knee architecture including a prosthetic knee implant, no anterior cruciate ligament and a corrected neutral axial alignment (no varus and valgus rotation). Despite these differences, the patient exhibited high functionality and had a normal appearing gait. Because of variability of knee architecture as well a s gait modifications at an individual level, the ability of the KAM to correctly predict compartmental load balance could also be variable. In order to properly asse s s changes in gait modifications, the calculation of simultaneous contact, muscle, and liga ment forces should be performed to understand the internal components that make up the overall loads a joint experiences. Gait modifications that aim at reducing the KAM may still offload the medial compartment of the knee. To be more certain of effectiven ess, the use of validated musculoskeletal models that include contact modelling are suggested for a more thorough and accurate gait analysis.
26 Table 2 1 Percent reductions (mean std) in peak and impulse values relative to normal gait. Measure Medial t h rust Walking p ole Knee a dduction m oment f irst p eak 32.4 (0.8)* 33.1 (0.9)* Knee a dduction m oment s econd p eak 14.9 (1.4) 47.2 (1.0)* Knee a dduction m oment a ngular i mpulse 37.8 (0.2)* 43.9 (0.2)* Medial c ontact f orce f irst p eak 4.5 (0.3) 4.7 (0.3) Media l c ontact f orce s econd p eak 19.9 (0.9)* 41.8 (0.7)* Medial c ontact force i mpulse 11.9 (0.2)* 28.0 (0.1)* D enotes statistically significant difference from normal trials.
27 Table 2 2. Prediction R 2 values between tested variables Adduction m oment f irst p eak Adduction m oment s econd p eak Adduction m oment a ngular i mpulse Medial f orce f irst p eak 0.30 0.00 0.03 Medial f orce s econd p eak 0.34 0.68 0.57 Medial f orce i mp ulse 0.38 0.75 0.56
28 Fig ure 2 1 Trial means of primary variables. A ) K nee ad duction mo ment B ) M edial contact force
29 Fig ure 2 2. Trial means of secondary variables. A ) K nee flexion moment B ) K nee flexion angle C ) H ip internal rotation angle.
30 Figure 2 3. Changes in primary peak and impulse values A) K nee adduction moment. B) M edia l contact force. ( D enotes statistically significant difference from normal trials )
31 Fig ure 2 4 Changes in secondary peak values A ) K nee flexion moment. B ) K nee flexion angle. C ) H ip internal rotation angle. ( D enotes statistically significant differ ence from normal trials )
32 Fig ure 2 5 R 2 value sensitivity to rotation of knee adduction axes.
33 CHAPTER 3 EVALUATION OF MUSCUL OTENDON MODELS WITH EMG BASED MUSCLE FORCE ESTIMATION Background Accurate assessment of human muscle forces during walking c ould significantly aid in the analysi s and treatment of common neuromusculo skeletal disorders such as osteoarthritis (Fregly et al., 2007), stroke (Higginson et al., 200 6), and cerebral palsy (Arnold and Delp, 2005). However, the inability to easily measur e muscle forces in vivo along with the inability to calculate unique muscle forces directly has greatly hindered the achievement of this goal. The neuromuscular system is redundant with a large number of muscles controlling the joints in the body and there fore, muscle force prediction through optimization is difficult. This muscle force indeterminacy problem has motivated the development of a variety of muscle force estimation methods. These methods typically seek to include as many constraints as possible to create a more unique set of estimated muscle forces. Constraints can come from a variety of sources including kinematic (marker motion, fluoroscopy), kinetic (ground reaction, implant loa ds), and EMG measurements. EMG driven models have sought to const rain muscle excitations directly by utilizing processed EMG signals as the excitation signal (Table 3 1). These models require an EMG to force model to account for the discrepancies found between experimental measurements of muscle electrical activity and measured net joint moments. However, the model components used to define the EMG to force relationship, and combinations of these components, varies greatly between studies. For example, one of the original studies that used a series of models to relate EM G data to muscle force (Hof and Van den Berg, 1981), modeled similar EMG to force
34 relationships that are modeled in more recent papers (2010; Kumar et al., 2012), but utilized very different equations to do so. The main model components that can be varied include parameter s controlling the excitation time delay, activation dynamics and nonlinearity, and muscle contraction dynamics. The excitation time delay is usually modeled as a pure time shift and can represent the small amount of time observed between t he onset of EMG activity and the onset of a measurable change in load. The time delay has also been modeled with a low pass filter that is not zero phase shifted and thus effectively applies a time delay. Activation dynamics, which represent the slower ris e and fall of muscle force compared to EMG activity, is typically modeled using either a first order (He et al., 1991; van den Bogert et al., 2011) or second order (Lloyd and Besier, 2003) set of ordinary differential equations ( ODE ) To represent the obse rved nonlinearity between processed EMG signals and measured joint moments, a model can be used that either uses one equation for the entire range of excitation or that separates the range into two parts, one which is nonlinear for low excitation and one w hich is linear for higher excitations (Table 3 1) The relationship between muscle activation and muscle force is known to be both length and velocity dependent and these are typically modeled with contraction dynamic equations based on the Hill type muscl e model (Zajac, 1989). A tendon element, in series with the muscle model, can also be included and is modeled as either a compliant or infinitely stiff spring (Table 3 1) There are only a small number of reported EMG based models that simulate all joints of the leg and are applied to walking tasks (Table 3 2). To our knowledge, these studies have not rigorously tested different model components and have not done so
35 with experimental validation measures The most recent and highest fidelit y models, (Barrett et al., 2007 ), have modeled the hip, knee and ankle joints simult aneously using EMG signals directly as muscle excitation. Minimization of the number of model parameters and determination of which model components significantly affect muscle force predict ions would improve the development of simple and computationally efficient EMG based models (Heine et al., 2003 ; Anderson and Pandy, 2001b ). This study investigates how the complexity of EMG to force models affects estimated muscle forces in the leg durin g walking. The investigation was performed using experimental motion capture, ground reaction, and EMG data collected from a subject implanted with a force measuring tibial prosthesis. Muscle forces were estimated for a complete walking cycle using a subje ct specific lower extremity musculoskeletal model constructed in OpenSim from full leg CT data available from the subject. The complexity of custom EMG to force models was varied by including or omitting activation dynamics, tendon compliance, and muscle f orce length velocity properties. For each EMG to force model investigated, static optimization was used to knee, and ankle), contact loads (knee), and processed EMG signals (hip, knee, and ankle) simultaneously across the gait cycle. Muscle force predictions generated by the different models were compared to assess how reduction of model complexity affected the estimated muscle forces. M ethods Experimental data (Fregly et a l., 2012) were collected from a single subject implanted with a force measuring knee replacement (female, left knee, age: 69 yrs, mass: 78 kg, height: 167 cm). Institutional review board approval and subject informed
36 consent were obtained prior to testing. The subject performed normal overground walking trials with simultaneous collection of GRF marker trajectories, EMG, and knee implant loads. All gait trials were performed at an average self selected walking speed of 1.3 m/s. The subject also performed j oint range of motion trials at the ankle, hip, and knee and a static standing pose to be used for calibration of model parameters of a rigid body, inverse dynamics model EMG was collected for the biceps femoris long head (BifemLH) semimembranosus (SemiM em) medial gastrocnemius (MedGas) lateral gastrocnemius (LatGas) tensor fascia lata (TFL) vastus lateralis (VasLat) rectus femoris (RF) sartorius (Sart) gracilis (Grac) soleus (Sol) tibialis anterior (TIbAnt) peroneus longus (PerLong) adductor ma gnus (AddMag) and gluteus medius (GlMed) Knee implant loads contained three orthogonal forces and three orthogonal moments measured at the center of the tibial tray (Kirking et al., 2006). Dynamic Musculoskeletal Model Creation A rigid body dynamics mod el was created of the subject using GRF and marker motion data. Joint centers, joint axes, and inertial parameters were optimized to reduce marker error and pelvis residuals (Reinbolt et al., 2005). Bone geometries of the e created from 0.67 mm slices of CT scans Contact model geometry was created from laser scans of the subjects knee implants and aligned to the bone geometries from CT. Visually estimated joint center locations at the hip (middle of femoral head) knee ( middle of the distal femoral end), and ankle ( middle of talus), as well as bony landmark marker locations ( anterior and posterior superior iliac spine, femoral and tibial condyles, second distal phalange of the foot, and calcaneal tuberosity ), were alig ned to the calibrated joint centers and static marker data to align the bony geometries to the rigid body dynamics model. A twelve DOF contact
37 model was used at the knee replacing the standard pin joint used for joint center and knee FE axis calibration (Z hao et al., 2007; Lin et al., 2010). A musculoskeletal model was created in OpenSim combining the knee contact model, bone geometries, and calibrated rigid body model. Muscle attachment and wrapping surface parameters were taken from a nominal OpenSim mod el based on human cadaver data (Ward et al., 2009; Arnold et al., 2010). Transformations found from scaling and aligning the bones between the nominal and subject's model s were used to adjust the muscle attachments and wrappings surfaces from the nominal m odel. Adjustments to these locations were needed to account for differences between the standard and subject specific geometries especially for the anatomical differences between male and female pelvic bone dim ensions. An additional locked joint with six D OF s was included at the tibial origin to calculate net joint moments that are mutually orthogonal and calculated about the center of the tibial tray where knee implant loads were measured. Muscles in the model include the BifemLH, BifemSH, LatGas, MedGas, SemiMem, SemiTen, TFL, VasMed, VasInt, VasLat, RF, Sart, Grac, TibAnt, Sol, AddBrev, AddLong, four portions of the AddMag, Gem, three portions of the GlMax, three portions of GlMed, three portions of GlMin, Iliac, Pect, Piri, Pso, QuadFem, EDL, EHL, FDL, F HL, PerBrev, PerLong, PerTert, TibPost (see Abbreviations for full muscle names). Kinematics Tibio femoral knee kinematics for each normal gait trial were developed using the twelve DOF contact model and static optimization. The S I force (Fy) ML force (Fz) and AA moment (Fx) measured from the instrumented kne e implant were matched with an elastic foundation based contact model (Zhao et al., 2008). The FE angle was
38 prescribed from inverse kinematics of the marker locations. Th e IE rotation angle was prescribed to be constant and determined from the average during stance from fluoroscopy based kinematics (Fregly et al., 2012). AP translation was also prescribed and based on fluoroscopy. The p atello femoral kinematics were calcul ated with static optimization such that all loads were balanced on the patella from contact, patellar ligament springs, and nominal guesses of the quadriceps muscle forces (Lin et al., 2010). With tibio femoral and patello femoral kinematics held constant an inverse kinematics analysis was performed in OpenSim (De lp et al., 2007) to determine three pelvis translations (tx, ty, tz) three pelvis rotations (rx, ry, rz) three hip rotations for each leg (rx, ry, rz) one contralateral knee rotation (rz) and two ankle rotations for each leg ( rx, rz ). A fifth order Woltring filter (Woltring, 1985) with a low pass cutoff frequency of 6 Hz was used to smooth the generalized positions and to determine the generalized velocities and accelerations for each degree of freedom of the model. Joint Loads A custom built inverse dynamics analysis was created in OpenSim (Delp et al., 2007) using prescribed translations, velocities, and accelerations along with inclusion of applied GRF and knee implant loads. GRF and knee co ntact loads were filtered with a zero lag fourth order Butterworth filter with a cutoff frequency of 6Hz to be consistent with both the kinematics and EMG data. Loads calculated with this approach represent the external inverse dynamic loads minus the cont ributions from internally applied implant loads. This net load should be balanced by the remaining internal structures which only include muscles for this study. The contributions of knee ligaments to the SI force and AA moment were inferred from measured contact loads during a knee laxity
39 trial test (Fregly et al., 2012) determined to be minimal relative to loads that occur during gait and were thus omitted from the current model. The effect of each muscle at the joint level is needed to determine the c ontributions that should balance out the net inverse dynamic loads. The Moment Arm Analysis in OpenSim (Delp et al., 2007) was used to calculate the amount of load applied to each degree of freedom of every joint produced by one unit of force from each mus cle individually. During simulation, each muscle's force production was multiplied by these unit contributions to determine the entire net load produced by the muscles for an individual gait trial. To account for differences between the subject specific mu sculoskeletal geometry and the generic model's muscle attachments, via points and wrapping surfaces, adjustments to these moment arms were parameterized by B spline curves with five nodes. The B spline coefficients for each moment arm for each muscle were optimized with a maximum adjustment of 1 cm. Passive joint moments were also applied at the hip and ankle based on equations from Davy and Audu ( 1987) with parameters adjusted to match the best fit of experimental data (Silder et al., 2007). Passive moment s were not included at the knee because of the lack of an ACL in the subject, which may contribute significantly to the passive flexion moments. The effect on moment arms and net inverse dynamic loads from modifying the AP translation of the knee was also quantified using offsets between 6 mm and 6 mm at 2 mm increments at each time frame. Offsets to the moment arms and net loads were fit with a straight line relative to the amount of AP translation so that the AP translation could be optimized without tim e intensive re calculation of moment arms and inverse dynamic net loads
40 EMG to Force Model An EMG to force model was created that described the relationship for each muscle between muscle excitation (which is linked to EMG patterns) and muscle force. EM G data was high pass filtered (zero lag fourth order Butterworth at 30 Hz), rectified, and low pass filtered (zero lag fourth order Butterworth at 6 Hz) to create a desired set of EMG/excitation profiles. Each muscle signal was normalized by the maximum pr ocessed EMG signal achieved throughout all normal gait trials maximum voluntary contraction trials and other maximum force evoking trials (Fregly et al., 2012). Excitation curves were parameterized by a set of 25 modifiable B spline coefficients. Each e xcitation curve was then normalized so that it ranged from zero to one, and a modifiable excitation gain parameter was used to scale each curve. The excitation gain parameter was included to efficiently control how much a muscle is excited relative to its maximum muscle force without needing to modify the individual B spline coefficients. The modification of the gain parameter was constrained such that the scaled excitation curve was similar in magnitude to the scaled EMG curve. This constraint especially r estricted these gains such that the scaled model excitations were not larger than the scaled EMG curve since this would be akin to decreasing the maximum EMG value measured. The excitation gain parameters were grouped between similar muscles to maintain a similar level of m uscle force and to thus avoid individual muscles from being turned off completely. The excitation scale parameters were grouped according to similar muscle function and/or innervation (such as the three Vastus muscles). The modification of gain parameters for muscles with EMG was introduced based on the assumption that EMG signals during gait are not perfectly
41 scaled when normalized by the maximum EMG magnitude during maximum voluntary contraction, and to allow for all muscles to be treat ed similarly. EMG to activation was modeled using time delays, activation dynamics and nonlinearities. A time delay was modeled by using either 1) a very small pure time delay to account for the electrical t ransmission time (at most 20 ms; Corcos et al., 1 992) or by 2) using a pure time delay model which allows the excitation to be moved forward in time by up to 150 ms (Winters and Stark, 1988; Jia et al., 2011) and therefore modelling both an electrical and mechanical delay. The larger electro mechanical time delay model was used only when neither activation dynamics nor contraction dynamics were included. A discretized form of activation dynamics was used which included a time constant for both activation and deactivation, ( Appendix A, He et al., 1991; va n den Bogert et al., 2011). Each time constant and time delay could be modified on a per muscle basis. Along with the need for a time delay and activation dynamics, the relationship between measured muscle EMG and force has been shown to be nonlinear. A tw o part model was used here based on experimental data which allowed for excitations from 0 to 0.3 to be modeled with a nonlinear exponential equation and for excitations above 0.3 to be modeled linearly (Manal and Buchanan, 2003). A single parameter for th is model ranged from 0 to 0.12 (Manal and Buchanan, 2003) and was modified on a per muscle basis. The passive and active musculotendon properties were modeled with two sets of musculotendon models based on Hill type muscle equations ( Appendix A, Zajac, 198 9). Each model included force length, force velocity, and pennation angle relationships using an active element in parallel with a passive element, both in series with a tendon
42 element. The nonlinear differential equations for the calculation of the applie d tendon force for these models were re expressed in a discretized form for computational efficiency and to avoid optimization of an initial force guess. A compliant tendon model included a series elastic tendon and thus required the use of a Newton Raphso n root finding algorithm and numerical approximations of time derivatives to calculate the unknown tendon force, and the muscle and tendon lengths. An infinitely stiff spring was used for the stiff tendon model. This model avoided the interaction between t endon length and muscle length and thus enabled a direct calculation of the applied tendon force to the body. Alternatively, when no contraction dynamics were considered, a pure force generator model was used that scaled activations by the muscle strength parameter to calculate muscle force. Musculotendon parameters, including tendon slack lengths and optimal muscle fiber lengths were restricted during optimization such that the difference between modeled and the scaled literature values was minimized Tend on compliance/stiffness gains were also allowed to vary during optimization in order to account for differences between muscles and variations away from the standard relationship for tendon stiffness as observed in previous s tudies for the Achilles tendon (Muramatsu et al., 2001 ). This allowed the parameter calibration process to control the amount of effective time delay introduced from a compliant tendon model The resulting normalized muscle fiber lengths were constrained during optimization such that th e minimum was at most 0.8 and the maximum was at least 0.7 (Garner and Pandy, 2003). The muscle strength parameters were initially grouped together such that they were a scaled version of the literature values using a single gain value. After one round of parameter calibrations, muscle strength values were allowed
43 to ch ange independently but deviation were minimized to be within +/ 15 % of the the muscle volumes, PCS A values, and the specific tension values found in literature (Holzbaur et al., 2007; Ward et al., 2009; Arnold et al., 2010). EMG Based Optimizations Optimizations were performed with design criteria of balancing all available net joint loads with muscl e forces matching muscle excitations to EMG signals for muscles with recorded EMG. A nonlinear least squares optimization routine was used in MATLAB The optimization cost function can be summarized as follows : The terms in this cost function represent the experimental tracking errors, E T the boundary errors, E B the parameter errors, E P and the muscular effort, E T The minimization of muscular effort was induced by including the excitation signals (after excitation gains were applied) with no weight and a summation of values squared: The cost function term used for tracking experimental net loads and excitations, including the square exponent performed by the optimizer, was:
44 N et load matching errors were quantified by R 2 values to match both th e net load shapes and magnitudes closely: Weights were chosen to scale this error so that the cost function error w as equal to one at the desired cutoff: Errors for EMG tracking were also represented by R2 values between the experimental EMG signals and the modeled excitations. The weights for this error were chosen such that the break po int would be at 0.95 for muscles with high EMG signal confidence, 0.90 for muscles with medium confidence and 0.85 for muscles with low confidence. Confidence level was based on the depth and size of muscle bellies as well as possibility of cross talk cont amination. VasInt and VasMed excitations were tracked to the VasLat EMG with a desired error of 0.75 since no EMG data was available for these muscles and because it was assumed that they would be similar but not exactly equivalent SemiTen excitations wer e tracked to the SemiMem EMG signal. To increase the disparity between errors above and below the desired R Square error values, an exponent was used which raised the scaled error to an even exponent
45 greater than two. If errors are scaled such that the sca led error is equal to one when the error is equal to the desired error, then, when raised to an exponent gre ater than one, the resulting cost function values for errors above the desired error will be increased and the values for errors below the desired e rror will be decreased. The exponent of this cost function term was selected by finding the even exponent that increased the final cost function value by a multiple of 3 when the error was greater th an the desired error by 5 % and similarly decreased the co st function value by a multiple of 3 when the error was below the desired error by 5 %. This rationale corresponded to the exponent of 20. For various calculated variables, including muscle activations and errors between modeled and reported muscle strength s, strict boundaries were desired to ensure that variables remained within realistic physiological limits. Several of these variables such as activations required both upper and lower bounds A monomial function with a large exponent was used for cost fun ction terms t o provide a single and continuous ly differentiable function that would greatly penalize values that are closely approaching and outside bounds, but would not (significantly) penalize values inside bounds First, a variable was offset by the av erage value of the desired range, for example, 0.5 was subtracted from activations which centered the error at zero. Then, the offset variable was scaled by one over half of the desired range so that the error was equal to one at the desired boundary value on either side F or example, offset activations were multiplied by such that activations of zero and one would correspond to errors of 1 and 1. Finally, the offset and scaled variable was raised to a large exponent to greatly increase the disparity between errors above and below the cutoff. This exponent was selected by comparing various options and choosing the even exponent which
46 increased the final error by a multiple of 3 when the variable exceeded the boundary by 0. 5 %, an d similarly decreased the error by a multiple of 3 when the variable was within the desired boundary by 0.5 % This rationale corresponded to an exponent of 200 The cost function term for maintaining boundaries is represented as follows: This equation represents, respectively, the excitation and activation bounds, deviations of muscle strengths from the scaled literature values, moment arm adjustment bounds, maximum bound for the minimum normalized muscle lengths, minimum boun d for the maximum normalized muscle lengths, and normalized force velocity coefficient bounds. The terms included in the cost function to minimize model parameter changes was: This equation represents, respectively, the design variable modification error, tendon slack and optimal muscle fiber length deviation from scaled literature values, and activation/deactivation time constant deviation from scaled literature values. Weights for
47 these terms were chosen such that calculated expression inside the parentheses would equal one at the desired boundary levels, similar to the experimental tracking scaling method. These boundary levels were set to 10% for the design variable modification weight, w 1 15% for the muscle tendon length w eights, w 2 and w 3 and 100% for activation time constant weights, w 4 and w 5 Initial Guess Initial excitation guesses for muscles with EMG were set to the processed EMG signal. For muscles without EMG, the excitation gains were set to a constant value of 0 .5 For muscles with EMG, these were set to the value needed to scale the excitation to be equal in magnitude to the scaled EMG signal. Excitation patterns were set to match EMG from similar muscles or from EMG found in the literature. Time delays were set to 20 ms for both time delay models. Time constants were set to those values found in literature (Winters and Stark, 1988 ). Nonlinear parameters were set to zero. Nominal pennation angle and muscle strengths were set to values found in the literature (Arn old et al., 2010). Optimal muscle fiber length and tendon slack length were both set to literature values (Arnold et al., 2010) that were modified to have normalized muscle fiber lengths that ranged from 0.8 to 1.2 throughout a single gait cycle. Paramete r Calibration As a first step in parameter calibration, sequential modification of individual parameter types was performed to improve convergence rates. Since interaction between the different parameters is also important, simultaneous optimization of eve ry combination of two design parameters was performed after individual calibration. Cost function errors for the calibrations also included a term that represented changes in the design variables from the initial guess to better guide each optimization pr oblem during
48 parameter calibration. Optimization of all parameters together, and then of each calibration was repeated until the overall cost function error was reduced by le ss than 1.0% between the averages of the last five and the previous five optimizations. Model Testing The main sets of model combinations were set up to determine the amount of model complexity required for realistic and reliable EMG to force modelling F ive cases of EMG to force models were tested, 1) pure electrical time delay with activation dynamics and compliant tendon contraction dynamics (Full), 2) pure electrical time delay with compliant tendon contraction dynamics (NoActDynams), 3) pure electrica l time delay with activation dynamics with stiff tendon contraction dynamics (StiffTendon), 4) pure electrical time delay with activation dynamics and a pure force generator muscle model (PureForce), and 5) pure electro mechanical time delay with a pure fo rce generator muscle model (DelayOnly). For each case, all model parameters were calibrated individually using one gait trial. Each case also included optimization of a nonlinear activation model, a moment arm adjustment model, and an AP knee translation a djustment model. The ability for each case to track EMG signals and to balance the inverse dynamics loads was tested by using an average error of all EMG to excitation R 2 errors and all net load R 2 errors. Muscle forces during calibration for each model we re compared to those from the fully complex (Full) model by taking the root mean square ( RMS ) difference for each muscle force and normalizing it by the maximum value of that muscle force from the full model After calibration of each model case with a sin gle gait trial, another gait trial was used to determine how well each model with calibrated model parameters could fit multiple data sets. Only excitation
49 patterns and gains were allowed to vary for these optimizations. The ability for each case to track EMG signals and to balance the inverse dynamics loads were also compared for these trials. R esults All model calibration cases reached a solution where overall average Net Load R 2 values were at least 0.98 with an overall standard deviation of 0.004 ( Figu re 3 3 ), and overall average EMG R 2 values were at least 0.87 with an overall standard deviation of 0.045 ( Figure 3 4 ). The testing trial cases achieved an overall average Net Load R 2 value of 0.98 with an overall standard deviation of 0.009 ( Figure 3 3 ) a nd an overall average EMG R 2 value of 0.83 with an overall standard deviation of 0.08 ( Figure 3 4 ). The case with the lowest tracking errors for the test trial was the NoDelay model with an average Net Load R 2 value of 0.98, standard deviation of 0.006, an d an average EMG R 2 value of 0.87, standard deviation of 0.05. The case with the largest tracking errors for the test trial was the PureForce model with an average Net Load R 2 value of 0.97, standard deviation of 0.009, and an average EMG R 2 value of 0.81, standard deviation of 0.099. The StiffTendon case had the highest overall muscle force differences from the Full model case for calibration with 25% normalized RMS for muscles without EMG and 12% normalized RMS for muscles with EMG. The DelayOnly model ha d the highest muscle force RMS difference for muscles with EMG, at 14%, and the lowest muscle force RMS difference for muscles without EMG, at 12%, ( Figure 3 5 ). D iscussion This study successfully estimated muscle forces that matched a large amount of exp erimental gait measurements for two gait trials using a musculoskeletal model of the
50 lower limb All five levels of complexity for the EMG to force model were able to match the net loads from an inverse dynamic analysis, and the EMG signals collected for a portion of the modeled muscles T he simplest model, DelayOnly, performed the best in that it resulted in lower net and EMG errors, and required the least amount of time to run. However the tracking errors for net loads and EMG signals were not different e nough between the worst and best model cases to conclude that any case represented a more superior model. Also, l arge differences between models were found in the resulting muscle forces for both calibration and testing trial s (Figure 3 6 ) These differenc es between models, along with the lack of variation in measurement tracking errors, does not support the hypothesis that the matching of all available experimental measurements creates a unique estimate of muscle force production even with the matching of knee contact loads A large amount of constraint was applied to the model, not only from net loads and EMG signals, but also from the use of the minimization of muscle excitations. Initial results showed significant differences in the estimated muscle fo rces for muscles without EMG signals when minimization of excitations in the cost function was not included For the case when excitations were not minimized, the smaller muscles without EMG signals tended to be left on at their maximum production througho ut the gait cycle. It is believed that the optimization algorithm would use these muscles as an easy path to balancing the net loads. These highly constrained muscle force predictions are unique due to experimentally known knee contact loads. Because thes e loads were matched at the knee, high fidelity simulations were created for model testing capabilities. However
51 these contact loads at the knee are not typical measurements available for subjects. Without these loads known, EMG based models may not accura tely estimate the correct amount of co contraction at the knee. Accurate knee kinematics may be important for muscle force estimation. Knowledge of the contact loads was also used to determine three of these kinematics, but AP translation was allowed to va ry during muscle force optimization. The resulting knee anterior posterior translation was similar to fluoroscopy and was found to be moderately helpful in reducing net load and EMG tracking errors (Figure 3 7). While the methods here utilized and matched the known contact loads, the modelling method presented is also capable of match ing the reduced set of net loads that is typically available in the leg. While the results of this study may have differed significantly if these loads were not matched, the in creased level of net load constraint created more unique muscle force estimates and thus a more difficult test for each EMG to force model. Further constraint, and muscle force uniqueness, may also come from a rigorously tested ligament model which was not available at the time of this study. Regression equations were used to model passive forces in this model. While these equations were based on experimental data, they may not be an accurate measurement of all net ligament loads present at the hip, knee, a nd ankle. The methods tested represent a range of model complexity supported by recent literature, and it was assumed that the most complex model (Full) was the most physiological model. Because no model was able to match the experimental measurements lar gely better than the others, and because each model case produced
52 different muscle forces than the most physiological model, it is suggested that the most physiological model should be used for future musculoskeletal modelling No other study has been abl e to model human gait with as high a level of accuracy to date. All tested EMG to force models were able to reproduce the experimental measurements but with a large variance in estimated muscle forces. Further research should be performed to determine the capability of these models to accurately estimate muscle forces.
53 Table 3 1. Literature of EMG based models for any joint and any task. Paper Time d elay Activation d ynamics Nonlinear Contraction d ynamics Amarantini and Martin, 2004 One Part Barett et al., 2007 Pure 2ndOrder One Part Compliant tendon Lloyd and Besier 2003 Pure 2ndOrder One Part Compliant tendon Bogey et al., 2005 2ndOrder Compliant tendon Buchanan et al., 2005 Pure 2ndOrder Two Part Compliant tendon Buchanan et al., 1998 Compliant tendon Besier et al., 2009 Pure 2ndOrder Two Part Compliant tendon Doorenbosch and Harlaar, 2003 Filter Gerus et al., 2010 Filter Stiff tendon Granata and Marras, 1995 Stiff tendon Hayashibe et al., 2009 Pure 2ndOrder One Part Compliant tendon Heine et al., 2003 Pure 2ndOrder Two Part Compliant tendon Heintz and Gutierrez Farewik, 2007 2ndOrder Compliant tendon Hof et al., 1987 Compliant tendon Hof and Van den Berg, 1981 Filter Exponent Compliant tendon Ja cobs et al., 1996 1stOrder Compliant tendon Jonkers et al., 2002 1stOrder Compliant tendon Kumar et al., 2012 Pure 2ndOrder One Part Compliant tendon McGill, 1992 Stiff tendon Nussbaum and Chaffin, 1998 Pure One Part Stiff tendon Piazz a and Delp, 1996 Pure 1stOrder Compliant tendon Shao et al., 2009 Pure 2ndOrder Two Part Compliant tendon White and Winter, 1992 Filter Stiff tendon
54 Table 3 2. Literature describing EMG based models for the full leg during gait. Paper Time d elay Activation d ynamics Nonlinear Contraction d ynamics Barrett et al., 2007 Pure 2ndOrder One Part Compliant tendon Heintz and Gutierrez Farewik, 2007 2ndOrder Compliant tendon Jonkers et al., 2002 1stOrder Compliant tendon Neptune et al., 2004 Pure 2ndOrder Compliant tendon Piazza and Delp, 1996 Pure 1stOrder Compliant tendon White and Winter, 1992 Filter Stiff tendon
55 Table 3 3. Net Load R 2 values for test trial Net Load Full DelayOnly NoActDynams StiffTendon NoMuscDynams Knee Fy 0.962 0.973 0.963 0.962 0.970 KneeTx 0.965 0.975 0.965 0.966 0.973 KneeTz 0.962 0.973 0.964 0.962 0.971 HipTx 0.980 0.986 0.985 0.984 0.983 HipTy 0.983 0.985 0.977 0.981 0.985 HipTz 0.980 0.985 0.981 0.982 0.982 AnkleTx 0.982 0.985 0.966 0.982 0.9 86 AnkleTz 0.984 0.987 0.984 0.983 0.984
56 Table 3 4. EMG R 2 values for test trial. Muscles are listed in order of increasing overall error. Muscle Full DelayOnly NoActDynams StiffTendon NoMuscDynams TibA nt 0.900 0.939 0.915 0.911 0.930 A dd M ag 0.895 0.918 0.904 0.901 0.913 S ol 0.897 0.918 0.886 0.926 0.900 PerL ong 0.860 0.920 0.876 0.873 0.902 LatG as 0.859 0.901 0.885 0.862 0.882 BifemLH 0.859 0.897 0.894 0.846 0.877 VasL at 0.848 0.892 0.857 0.864 0.878 MedG as 0.848 0.854 0.876 0.875 0.845 S emi M em 0.846 0.857 0.880 0.819 0.851 G l M ed 0.805 0.863 0.831 0.829 0.836 RF 0.788 0.857 0.818 0.756 0.841 Sart 0.683 0.799 0.745 0.676 0.745 TFL 0.681 0.790 0.731 0.669 0.749 Grac 0.667 0.790 0.716 0.664 0.731
57 Figure 3 1. Exponent selection method f or experimental tracking errors. Figure 3 2 Exponent selection method for boundary errors.
58 Figure 3 3 Average Net Load matching R 2 values for calibration and prediction. Error bars indicate one standard deviation throughout 8 matched net loads.
59 Figure 3 4 Average EMG matching R 2 values for calibration and test trials. Error bars indicate one standard deviation throughout 14 matched EMG signals.
60 Figure 3 5 Average muscle forces differences to full model for calibration trial. RMS differen ce normalized to maximum force from all cases for each muscle. Error bars indicate one standard deviation.
61 Figure 3 6 Comparison of muscle forces between models for calibration trial. Muscles chosen based on largest moments for muscles with EMG and wi thout. Figure 3 7 Optimized knee anterior posterior translation.
62 CHAPTER 4 EVALUATION OF NEURAL SYNERGIES FOR MUSCLE FORCE E STIMATIONS Background Biomechanical models of the human body are an essential tool for the analysis of dynamic movement. These models have been used in a variety of applications including sports training, film production, injury and crash analysis, and the scientific study of several health problems such as osteoarthritis (Fregly et al., 2007), stroke (Higginson et al., 2006), and cerebral palsy (Arnold and Delp, 2005). Gait related disorders are studied by modelling the dynamic interactions between the muscular, skeletal, and nervous systems of the human body. The outputs of these biomechanical analyses help guide physical therapi sts, orthopaedic surgeons and developers of surgical implants to better design treatments for patients. Furthermore, as these models become more sophisticated, treatments in the future can become more subject specific and thus increase the likelihood of tr eatment success and longevity. Research using these neuromusculo skeletal models has been limited however by the inability to easily measure in vivo muscle forces and the lack of validation of current techniques (Erdemir et al., 2007). The human body has more actuators than degrees of freedom at most joints and is thus a redundant mechanical system. This disables the use of simple dynamic equations to calculate the muscle forces needed to actuate a measured motion. The most common method to solve this prob lem is to utilize optimization methods to determine an optimal set of muscle forces based on muscular effort or energy (Piazza and Delp, 1996; Neptune and Hull, 1998; Anderson and Pandy, 2001 a ). Another common method is to constrain these techniques using EMG, which essentially measures the relative amount of neural input to muscles during dynamic
63 activities (Lloyd and Besier, 2003; Buchanan et al., 2005). While EMG driven models have become the forefront of experimentally guided muscle force estimations, t hey are limited by the inability to easily measure excitations for all of the muscles and by the inaccuracies introduced from electrical crosstalk and noise (De Luca, 1997; Disselhorst Klug et al., 2009). It has been proposed that the neurological system in the body attempts to solve the muscle redundancy problem by sending out a limited number of time varying commands from the central nervous system (Dimitrijevic et al., 1998; MacKay Lyons, 2002; Ivanenko et al., 2004; Ting and Macpherson, 2005; Kutch and Valero Cuevas, 2012). This concept of neural synergies has been proposed to constrain the muscle activity of all muscles based on a small set of time varying neural commands and a small set of s ynergistic pathways or synergy vectors (Taga, 1995; Byadarhal y et al., 2012; Kutch and Valero Cuevas, 2012). The neural commands may be a set of feed forward signals sent from a central pattern generator which are combined in synergistic pathways to branch out into multiple muscles. E ach muscle would have one set of synergy vectors that directs relative portions of each neural command to that muscle. If synergy commands and vectors are calculated from experimentally measured EMG data from a subject being analyzed, this may allow for excitations from all muscles to be constrained with experimentally derived neural constraints. Previous work has shown that gait can be simulated by five to six time varying signals that are combined with experimentally derived synergy vectors for specific muscles. H owever these studies d id not allow excitations of all muscles to be parameterized with neural synergies and did not have means to validate estimated muscle and contact forces (Neptune et al., 2009;
64 Allen and Neptune, 2012). It is hypothesized for the current study, that this no vel addition to the muscle force estimation process will create a muscle force estimation technique capable of accurately re producing overall joint loads at the knee. The goal of this study was to evaluate whether parameterization of muscle excitations wi th experimentally derived synergy relationships and/or tracking shapes of EMG signals produces unique muscle force and accurate contact force estimates. The investigation was performed using experimental motion capture, ground reaction, and EMG data collec ted from a subject implanted with a force measuring tibial prosthesis. Muscle forces were estimated for a complete walking cycle using a subject specific lower extremity musculoskeletal model constructed in OpenSim from full leg CT data available from the subject. The level of constraint was varied by including or omitting the parameterization of muscle excitations with neural synergies, and the constraint of matching the shapes of muscle excitations to EMG signals. For each level of constraint investigated static optimization with minimized muscle excitation was used to evaluate them predicting them with pre calibrated parameters, and predicting them without pre calibrated parameters M etho ds Experimental data (Fregly et al., 2012) were collected from a single subject implanted with a force measuring knee replacement (female, left knee, age: 69 yrs, mass: 78 kg, height: 167 cm). Institutional review board approval and subject informed consen t were obtained prior to testing. The subject performed normal overground walking trials with simultaneous collection of GRF marker trajectories, EMG, and knee implant loads. All gait trials were performed at an average self selected walking speed of 1.3 m/s. The subject also performed joint range of motion trials at the ankle, hip, and
65 knee and a static standing pose for calibration of a rigid body inverse dynamics model EMG was collected for the BifemLH SemiMem MedGas LatGas TFL VasLat RF Sart Grac Sol TibAnt PerLong AddMag and GlMed Knee implant loads contained three orthogonal forces and three orthogonal moments measured at the center of the tibial tray (Kirking et al., 2006). EMG data were high pass filtered (zero lag fourth order Butte rworth at 30 Hz), rectified, and low pass filtered (zero lag fourth order Butterworth at 6 Hz) to create a desired set of EMG/excitation profiles. Ground reaction and knee contact loads were filtered with a zero lag fourth order Butterworth filter with a c utoff frequency of 6 Hz to be consistent with both the kinematics and EMG data. Dynamic Musculoskeletal Model Creation A rigid body dynamics model was created of the subject using GRF and marker motion data. Joint centers, joint axes, and inertial paramete rs were optimized to reduce marker error and pelvis residuals throughout joint range of motion trials, a static trial, and one gait trial pelvis were created from 0.67 mm slices of CT scans Contact model geometry of the knee was created from laser scans of the subjects knee implants and aligned to the bone geometries from CT. Visually estimated joint center locations at the hip, knee and ankle as well as bony landmark marker locations were aligned to the calibrated joint centers and static marker data to align the bony geometries to the rigid body dynamics model. A twelve DOF contact model was used at the knee replacing the standard pin joint used for joint center and knee flexion axis calibration (Zhao et al., 2007; Lin et al., 2010). A musculoskeletal model was created in OpenSim combining the knee contact model, bone geometries, and calibrated rigid body model. Muscle attachment and
66 wrapping surface parameters were taken from a nomin al OpenSim model based on human cadaver data (Ward et al., 2009; Arnold et al., 2010). Transformations found from scaling and aligning the bones between the nominal and subject's model were used to adjust the muscle attachments and wrappings surfaces from the nominal model. An additional locked joint with six DOFs was included at the tibial origin to calculate net joint moments that are mutually orthogonal and calculated about the center of tibial tray where knee implant loads were measured. Muscles in the model include d the BifemLH, BifemSH, LatGas, MedGas SemiMem SemiTen, TFL, VasMed VasInt VasLat R F Sart Grac TibAnt Sol AddBrev AddLong, four portions of the AddMag, Gem three portions of the GlMax three portions of GlMed three portions of GlM in Iliac Pect Piri Pso QuadFem EDL EHL FDL FHL PerBrev PerLong PerTert TibPost Kinematics Tibio femoral knee kinematics for each normal gait trial were developed using the twelve DOF contact model and static optimization. The SI force, ML for ce and AA moment measured from the instrumented kn ee implant were matched with an elastic foundation based contact model (Zhao et al., 2008). The FE angle was prescribed from inverse kinematics of the marker locations. The IE rotation angle was prescribed to be constant and determined from the average during stance from fluoroscopy based kinematics (Fregly et al., 2012). AP translation was initially prescribed based on the same fluoroscopy but adjusted during muscle force optimizations The p atello femoral kinematics were calculated with static optimization such that all six net loads were balanced on the patella from contact, patellar ligament springs, and nominal guesses of the quadriceps muscle forces (Lin et al., 2010).
67 With tibio femoral and patello f emoral kinematics held constant, an inverse kinematics analysis was performed in OpenSim (Delp et al., 2007) to determine three pelvis translations (tx, ty, and tz) three pelvis rotations (rx, ry, and rz) three hip rotations for each leg ( rx, ry, and rz ) one contralateral knee rotation (rz) and two ankle rotations for each leg ( rx, and rz ). A fifth order Woltring filter (Woltring, 1985) with a low pass cutoff frequency of 6 Hz was used to smooth the generalized positions and to determine the generalized velocities and accelerations for each degree of freedom of the model. Joint Loads An inverse dynamics analysis was used in the OpenSim environment using prescribed translations, velocities, and accelerations along with inclusion of applied GRF and knee i mplant loads. Loads calculated with this approach represent the external inverse dynamic loads minus the contributions from internally applied implant loads. This net load should be balanced by the remaining internal structures which only include musc les f or this study. The MomentArm Analysis in OpenSim (Delp et al., 2007) was used to calculate the amount of load applied to each DOF of every joint produced by one unit of force from each muscle individually. Medial, lateral and total contact forces at the kn ee were calculated using regression equations derived from experimental data using the modeled subjects instrumented knee implant loads and an elastic foun dation based contact model Contact forces were calculated for the model by subtracting the contribut ions of muscle forces from the net loads from inverse dynamic analysis and applying the regression equation that transformed three sensitive knee loads into medial, lateral and total contact force. Using these regression equations
68 and the net load contribu tions from each muscle, the contribution of each knee muscle to medial, lateral, and total contact forces were also calculated. Neural Synergy Factorization Neural synergy commands and vectors were calculated using a non negative matrix factorization (NMF) (Lee and Seung, 1999) for the one gait cycle that was analyzed To allow equal weightings between all muscles, each muscle signal was normalized by the maximum processed EMG signal achieved throughout the gait cycle analyzed. Synergy commands and vectors were initially set to random values and an alternating least squares algorithm was performed. The number of synergies was increased until average variance accounted for (VAF) was at least 0.95 and the addition of a synergy did not increase the average VAF by more than 3% (Cappellini et al., 2006). Six time varying commands were found from this analysis and were used for excitation parameterization during simulations. EMG to Force Model Excitation curves were parameterized by either a set of 25 modifiable B spline coefficients or by a combination of time varying neural commands controlled by synergy vectors. Synergy vectors for all muscles were treated as unknown parameters to be optimized. Synergy commands could also be optimized but were restricted to be ve ry similar (R 2 > 0.95) to the commands calculated from NMF using the experimental EMG signals of the gait trial used Each excitation curve was then normalized so that it ranged from zero to one, and a modifiable excitation gain parameter was used to scale each curve. The excitation scale parameters were grouped according to similar muscle function and/or innervation (such as the three Vastus muscles).
69 EMG to activation was modeled using time delays, activation dynamics and nonlinearities. A time delay was modeled by either using a very small pure time delay to account for the electrical t ransmission time (at most 20 ms; Corcos et al., 1992) The delay was optimized for each muscle but minimized to be close to the average between all muscles. A discretized f orm of activation dynamics w as used which included a time constant for both activation and deactivation, ( Appendix A ). Each time constant and time delay was optimized for each but minimized to be close to the average between all muscles. A two part nonline arity model was used here based on experimental data which allowed for excitations from 0 to 0.3 to be modeled with a nonlinear exponential equation and for excitations above 0.3 to be modeled linearly (Manal and Buchanan, 2003). A single parameter for thi s model that ranged from 0 to 0.12 ,(Manal and Buchanan, 2003), was optimized for each muscle but minimized to be close to the average between all muscles. The passive and active musculotendon properties were modeled based on Hill type muscle equations ( App endix A, Zajac, 1989). The contractile element included an active element with force length, force velocity, and pennation angle relationships in parallel with a passive element. The contractile element (muscle) was in series with a compliant linearly elas tic element ( tendon ) Musculotendon parameters, including tendon slack lengths and optimal muscle fiber lengths were bounded and constrained during optimization to vary between +/ 15% of scaled literature values. Tendon compliance/stiffness gains were also allowed to vary during optimization and minimized to be close to the average between all muscles. Muscle strength values were allowed to change independently within +/ 10% of scaled literature values.
70 EMG Based Optimizations Optimizations were performed with design criteria of balancing all available net joint loads with muscle forces matching muscle excitations to EMG signals for muscles with recorded EMG. A nonlinear least squares optimization routine was used in MATLAB The optimization cost function can be summarized as follows : The terms in this cost function represent the experimental tracking errors, E T the boundary errors, E B the parameter errors, E P and the muscular effort, E T The minimization of muscular effo rt was induced by including the excitation signals (after excitation gains were applied) with no weight and a summation of values squared: The cost function term used for tracking experimental net loads and excitations, includ ing the square exponent performed by the optimizer, was: N et load matching errors were quantified by R 2 values to match both th e net load shapes and magnitudes closely:
71 Weights were chosen to scale this error so that the cost function error was equal to one at the desired cutoff: Errors for EMG tracking were also represented by R2 values between the experimental EMG signals and the modeled excitations. The weights for this error were chosen such that the break point would be at 0.95 for muscles with high EMG signal confidence, 0.90 for muscles with medium confidence and 0.85 for muscles with low confidence. Confidence level was based on the depth and size of muscle bellies as well as possibility of cross talk contamination. VasInt and VasMed excitations were tracked to the VasLat EMG with a desired error of 0.75 since no EMG data was available for these muscles and because it was assumed that they would be similar but not exactly equivalent SemiTen excitations were tracked to the SemiMem EMG signal. To increase the disparity between errors above and below the desired R Square error values, an exponent was used which raised the scaled error to an even exponent greater than two. If errors are scaled such that the scaled error is equal to one when the error is equal to the desired error, then, when raised to an exponent gre ater than one, the resulting cost function values for errors above the des ired error will be increased and the values for errors below the desired error will be decreased. The exponent of this cost function term was selected by finding the even exponent that increased the final
72 cost function value by a multiple of 3 when the err or was greater th an the desired error by 5 % and similarly decreased the cost function value by a multiple of 3 when the error was below the desired error by 5 %. This rationale corresponded to the exponent of 20. For various calculated variables, including muscle activations and errors between modeled and reported muscle strengths, strict boundaries were desired to ensure that variables remained within realistic physiological limits. Several of these variables such as activations required both upper and low er bounds A monomial function with a large exponent was used for cost function terms t o provide a single and continuous ly differentiable function that would greatly penalize values that are closely approaching and outside bounds, but would not (significan tly) penalize values inside bounds First, a variable was offset by the average value of the desired range, for example, 0.5 was subtracted from activations which centered the error at zero. Then, the offset variable was scaled by one over half of the desi red range so that the error was equal to one at the desired boundary value on either side F or example, offset activations were multiplied by such that activations of zero and one would correspond to errors of 1 and 1. Finally the offset and scaled variable was raised to a large exponent to greatly increase the disparity between errors above and below the cutoff. This exponent was selected by comparing various options and choosing the even exponent which increased the final er ror by a multiple of 3 when the variable exceeded the boundary by 0. 5 %, and similarly decreased the error by a multiple of 3 when the variable was within the desired boundary by 0.5 % This rationale corresponded to an exponent of 200 The cost function ter m for maintaining boundaries is represented as follows:
73 This equation represents, respectively, the excitation and activation bounds, deviations of muscle strengths from the scaled literature values, moment arm adjustment bou nds, maximum bound for the minimum normalized muscle lengths, minimum bound for the maximum normalized muscle lengths, normalized force velocity coefficient bounds neural command shape tracking error, and neural command bounds The terms included in the c ost function to minimize model parameter changes was: This equation represents, respectively, the design variable modification error, tendon slack and optimal muscle fiber length deviation from scaled literature values, and ac tivation/deactivation time constant deviation from scaled literature values. Weights for these terms were chosen such that calculated expression inside the parentheses would equal one at the desired boundary levels, similar to the experimental tracking sca ling
74 method. These boundary levels were set to 10% for the design variable modification weight, w 1 15% for the muscle tendon length weights, w 2 and w 3 and 100% for activation time constant weights, w 4 and w 5 Initial Guess Initial excitation guesses fo r all muscles were set to random values ranging from zero to one. Synergy vectors were set to random values ranging from zero to one. Synergy commands were set to the commands calculated from NMF. Excitation gains were set to 0.5 for muscles without EMG an d to the corresponding gain needed to scale the excitation to the level of the scaled EMG signal for muscles with EMG. Time delays were set to 20 ms. Time constants were set to values found in literature (Winters and Stark, 1988). Nonlinear parameters were set to 0.85 based on the average of a previously calibrated musculoskeletal model of the same subject. Nominal pennation angle and muscle strengths were set to values found in literature (Arnold et al., 2010). Optimal muscle fiber length and tendon slack length were both set to literature values (Arnold et al., 2010) that were modified to have normalized muscle fiber lengths that ranged from 0.8 to 1.2 throughout a single gait cycle. Parameter Calibration Model parameters and excitations were optimized to match optimization design criteria. As a first step in parameter calibration, sequential modification of individual parameter types was performed to improve convergence rates. Since interaction between the different parameters is also important, simultane ous optimization of each design parameter with excitation shapes was performed after individual calibration. Cost function errors for the calibrations also included a term that represented changes in the
7 5 design variables from the initial guess to better gu ide the solution to a better solution similar to having a slow cooling rate in a simulated annealing method (Appendix B) alone were performed last. Optimization s of a m om ent arm adjustment model, and a knee AP translation adjustment model were also included throughout the calibration process. Parameter and excitation pattern calibration was repeated until the overall cost function error was reduced by less than 1.0% betw een consecutive rounds of optimizations. Model Testing The main sets of model combinations were set up to determine the amount of constraint required for accurate contact force predictions. Four cases of constraint were tested: 1) minimize excitations usi ng B splines for excitations (Bsp), 2) minimize excitations and track EMG using B splines for excitations (BspEMG), 3) minimize excitations and track EMG using synergies to parameterize excitations (SynEMG), and 4) minimize excitations using synergies to p arameterize excitations (Syn). For each case, three subtests were performed to evaluate the necessity of well calibrated parameters : 1) all model parameters calibrated and excitations optimized while matching three knee loads with contact contributions inc luded (Match), 2) using the parameters pre calibrated in the first subtest, reset and re optimize d excitations to match only one knee load without contributions from contact (PreCalPred), and 3) all model parameters calibrated and excitations optimized whi le matching only one knee load without contributions from contact (Pred). The ability for each case to track or predict medial and lateral contact forces, to track or predict EMG signals, and to
76 balance the inverse dynamics loads was tested by using an av erage error of all contact force R 2 and RMS errors, EMG to excitation R 2 errors and all net load R 2 errors. R esults Net load tracking R 2 values ranged from 0.96 to 0.98 with an average of 0.97 and an overall standard deviation < 0.01 (Figures 4 1, 4 2, 4 3 4 4). Individual net load RMS errors normalized to the RMS values of desired net loads ranged from 8.8% to 13%, with an average of 10% and an overall standard deviation of 1.0%. For cases that tracked EMG, EMG tracking R 2 values ranged from 0.71 to 0.95 with an average of 0.86 and an overall standard deviation of 0.05 (Figures 4 5, 4 6, 4 7). For cases that tracked EMG, individual EMG RMS errors, normalized to the max values of desired EMG signals, ranged from 6.7% to 14%, with an average of 9.9% and an overall standard deviation of 1.6%. For cases that did not track EMG and used synergies for parameterization, EMG tracking R 2 values ranged from 4.5 to 0.85 with an average of 0.05 and an overall standard deviation of 1.1. For the same cases, individual EMG RMS errors, normalized to the max values of desired EMG signals, ranged from 10% to 47%, with an average of 23% and an overall standard deviation of 7.6%. For cases that did not track EMG and used B splines for parameterization, EMG tracking R 2 values ranged from 3.9 to 0.70 with an average of 0.32 and an overall standard deviation of 1.1. For the same cases, individual EMG RMS errors, normalized to the max values of desired EMG signals, ranged from 14% to 53%, with an average of 28% and an overall s tandard deviation of 9.1%. For cases where knee contact loads were matched (Match), the average and standard deviation of R 2 values for medial, lateral, and total contact force were 0.97
77 (<0.01), 0.96 (<0.01), and 0.99 (<0.01), respectively, (Figures 4 8, 4 9, 4 10, 4 11). For cases where pre calibrated parameters were used and knee contact loads were predicted (PreCalPred), the average and standard deviation of R 2 values for medial, lateral, and total contact force were 0.92 (0.01), 0.83 (0.07), and 0.90 (0.03), respectively. For cases where parameters were calibrated and knee contact loads were predicted (Pred), the average and standard deviation of R 2 values for medial, lateral, and total contact force were 0.83 (0.05), 0.76 (0.11), and 0.82 (0.06), resp ectively. The highest averages of R 2 values between medial, lateral, and total contact forces, for contact force predictions were achieved by the synergy without EMG tracking case (Syn) with a margin of <0.03 R 2 for PreCalPred cases and by Syn with a lar ger margin of 0.10 R 2 for Pred cases The median R 2 value of all of contact force estimates with synergies (Syn and SynEMG, R 2 = 0.92) was higher than the median for estimates with B splines (Bsp and BspEMG, R 2 = 0.88). The median R 2 value of all of contac t force s without EMG tracking (Syn and Bsp, R 2 = 0.91) was higher than the median for estimates with EMG tracking (SynEMG and BspEMG, R 2 = 0.89). D iscussion This study tested the ability to match and predict contact forces at the knee with a musculoskeleta l model while comparing the use of neural synergies and EMG tracking. All methods tested were able to match and predict knee contact forces well while matching net loads from an inverse dynamic analysis at the hip, knee, and ankle. The best contact force p redictions, where knee contact loads were not explicitly matched by the optimization algorithm, were produced by the method with synergy parameterization and no EMG tracking. This case was effective at predicting the first and second peak of
78 the medial, la teral, and total contact force during the stance phase of gait. These results support the hypothesis that the use of synergies to constrain all muscle excitations aids in the prediction of contact forces. However the use of synergies only decreased the RMS errors of contact forces predictions relative to the use of B splines when EMG signals were not explicitly tracked by the optimization algorithm. EMG tracking tended to create an underestimation of the second peak of contact forces at the knee during stan ce (Figure 4 11); even when model parameters were pre calibrated to match contact forces (Figure 4 10). The use of EMG tracking also decreased the overall ability to match the lateral contact forces for all cases (Figure 4 8). This underestimation of cont act forces during contact force prediction can be attributed to decreased muscle forces at the knee during the second half of stance especially when parameters were simultaneously calibrated with contact force predictions, (Figures 4 13, 4 14). Four muscle s had both relatively large contributions to the knee contact forces and relatively large differences between the main test methods. These muscles were the LatGas, MedGas, VasMed, and RF. The shape of the MedGas EMG signal was estimated poorly by both syne rgies and B spline parameterization when the EMG signal was not explicitly tracked by the optimization algorithm. This muscle was the largest contributor to the SI force and AA moment at the knee (and also the medial, lateral and total contact forces) duri ng the second half of stance. During preliminary investigations it was found that when muscles were tied to EMG signals and parameters were adjusted to match net loads, the maximum muscle fo rce produced tended to decrease. This occurs as an effort by the optimization algorithm to reduce EMG tracking errors without undesired contribution s to
79 the net load tracking errors from those muscles Thus when MedGas excitation is tracked to the EMG signal, the optimizer will tend to estimate smaller muscle forces, an d thus contact forces will also be reduced when this muscle is active (second half of stance). This relationship between tracking EMG and decreased muscle forces was not as strongly pronounced when pre calibrated parameters were used however. Thus, proper calibration of the EMG to force model parameters may be essential for accurately estimating muscle and contact forces during gait. The RF was another muscle for which the EMG signal was estimated poorly by both synergies and B spline parameterization when the EMG signal was not explicitly tracked by the optimization algorithm. This muscle also contributes a relatively large amount to the medial, lateral, and total force during the second half of stance. The second peak of this muscles force was reduced sub stantially during EMG tracked predictions without pre calibrated parameters which contributed to the underestimation of the second peak of the contact forces. This also occurred when pre calibrated parameters were used but muscles were not constrained by s ynergies. This seems to indicate that both correctly calibrated parameters of the EMG to force model and the use of synergies may be essential for estimating muscle and contact forces during gait. The use of synergies to parameterize all muscle excitation s increased the ability to estimate experimental EMG signals for a majority of muscles, (when not explicitly matched in the cost function). Certain muscles were not estimated well by either Bsp or Syn including the GlMed and MedGas. The GlMed EMG signal wa s suspected to be a bad collection of the muscles excitation due to a very low signal during gait trials. This could be due to signal reduction from the amount of tissue between the muscle and the
80 electrode as well from movement of the electrode relative to the excitation center of the muscle during gait. Overall, the general inability to reproduce EMG signals may be attributed to the lack of excitation constraint when using B splines (Bsp PreCalPred) or the lack of well calibrated parameters (BspPred and SynPred). Furthermore, because the use of synergies without EMG tracking and with pre calibrated parameters did not reproduce all EMG signals well, this may indicate that either the method of synergy parameterization used did not provide enough constraints or that the EMG signals of certain muscles were not perfectly collected from surface electrodes for this patient during gait (and thus should not be tracked) The matching of contact forces at the knee explicitly also did not always constrain the muscle f orces to be unique (Figure 4 12), supporting a general trend of not enough or improper constraint being applied to the muscle force estimation process from EMG tracking alone. Limitations include the exclusion of knee ligaments which may contribute to knee contact forces, especially during early swing phase. Experimental measurement of the contact loads during a passive knee flexion trial resulted in small contact forces except when fully flexed. This implies that during a normal gait trial when the knee is not flexed greatly, the ligament loads may not be large and may thus be ignored. Inclusion of generic ligaments at the knee resulted in applied loads only during the very end of the stance phase and early swing phase, indicating that ligament loads may no t contribute greatly to the remainder of the stance phase. Also for this subject, a model of the anterior cruciate ligament was not needed since it was removed during knee implantation. Another limitation that may affect these results is that only one gai t trial was used for parameter calibration, instead of maximum voluntary contraction trials.
81 Applicability of the calibrated model parameters to other motions would be necessary to determine the uniqueness of these parameters, however due to the long amoun ts of time needed for parameter optimizations, this test was not performed. Additional subjects may have provided more distinction between the methods tested yet because of the requirement of experimentally collected knee loads, data from only one subject was available for analysis. Since good estimates of the contact force were predicted ( best Pred average R 2 = 0.88 best Pred average RMS ~ 8% of peak force) it is believed that this is due to the use of highly patient specific musculoskeletal geometry. F urther testing is needed to determine how the contact force estimates change when a scaled generic model is used. Kinematics of the knee, including ML and SI translation, were also set as constants in the estimation process. These were set to match the exp erimental con tact forces from the gait trial analyzed. Thus, some of the knee kinematics were in effect pre calibrated. While this model estimation process optimized the AP translation of the knee, the ability of synergies and EMG tracking methods to predi ct knee kinematics correctly would also be helpful. The optimization of the six basic time varying neural synergy commands increased the time of optimizations by a large amount since a single change in one command creates the need to re calculate model out puts for all muscles. While this was found to be very helpful when matching contact forces along with tracking EMG and net loads, it may not be necessary for good contact force estimation. Finally, weighting of the net load and EMG tracking was also found to affect predictions of contact forces. Because of this, weights were chosen to produce generally good results during initial testing.
82 T his study was effective at reproducing a large amount of experimental measurements and at predict ing contact forces at the knee during gait using a musculoskeletal model The results seem to indicate that the use of synergies to constrain muscle excitations increased the ability to accurately predict contact forces. This small trend for better contact force predictions wa s found with the use of synergies to constrain muscle excitations but without the use of EMG tracking. It is believed that EMG based constraint met hods suffer from the difficulties excitation signal with EMG model parameter un certainty, and the lack of excitation magnitude constraints. Further research should be performed to determine the best model parameter calibration techniques and the effects of ligament modeling knee kinematics, and subject specific modeling
83 Figur e 4 1. Net load tracking quantities Thin lines indicate range of all values. Bars indicate range of the second and third quartiles (middle 50%) Plus sign indicates median. Normalized to the RMS value of desired net load.
84 Figure 4 2. Net loads modeled with desired/actual net load for Match cases
85 Figure 4 3 Net loads modeled with desired /actual net load for PreCalPred cases
86 Figure 4 4 Net loads modeled with desired /actual net load s for Pred cases
87 Figure 4 5 EMG tracking quantities at optimal time shifts. Thin lines indicate range of all values. Bars indicate range of the second and third quartiles (middle 50%) Plus sign indicates median. Normalized to RMS values of desired EMG.
88 Figure 4 6 EMG tracking normalized RMS errors at optimal time shifts. Normalized to RMS values of desired EMG.
89 Figure 4 7 EMG tracking R 2 values at optimal time shifts.
90 Figure 4 8. Contact force quantities for medial, lateral and total forces.
91 Figure 4 9. Contact forces modeled and desired /actual for Match cases Figure 4 10. Contact forces modeled and desired /actual for PreCalPred cases
92 Figure 4 11. Contact forces modeled and desired /actual for Pred cases Figure 4 12. Muscle forces with large differences and large contribut ions to contact forces for Match cases
93 Figure 4 13. Muscle forces with large differences and large contributions to contact forces for PreCalPred cases Figure 4 14. Muscle forces with large differences and large contributions to contact forces for Pred cases
94 CHAPTER 5 CONCLUSIONS The first study found that the standard method for gait analysis evaluation of medial knee osteoarthritis was not a good indicator of the experimental contact forces in the knee. Therefore, in order to properly asse ss changes in gait modifications, the calculation of contact and muscle forces is desired to evaluate the treatments for osteoarthritis and other musculoskeletal disorders. To be able to calculate these variables, a validated musculoskeletal model is neede d to understand, without uncertainty, the balance of loads present at joints. The musculoskeletal model and muscle force estimation methods presented here were able to simulate human gait with the highest level of accuracy to date. All tested EMG to force models were able to reproduce the experimental measurements but a large variance occurred in estimated muscle forces between methods. One goal of these models is to obtain rapid estimation of muscle forces for real time analysis of human movements. A simp le model is desired that takes less computational time and is still capable of estimating muscle and contact forces accurately. However, many researchers believe that a more physiological model, and thus more complex and time consuming, is necessary to obt ain accurate measurements, especially when models are calibrated with one type of motion and applied to other types of motion. Further research should therefore be performed to determine the capability to accurately estimate muscle forces of the simplest a nd most complex EMG to force models in a variety of motions The additional level of constraint that was applied to the EMG based optimizations process with neural synergies was found to be an effective method to simplify the
95 parameterization of muscle exc itations, reproduce the large amount of experimental measurements, and to predict knee contact forces well. However, only minor increases in the ability to predict contact forces resulted from the addition of synergy parameterization. Conversely to neural synergies, the inclusion of EMG tracking constraints, without properly calibrated parameters, was found to decrease the ability to predict contact forces, especially for the lateral force. This may indicate that either certain EMG signals which were tracke d should not have been tracked, that the full physiological model is incomplete or incorrect, or that new methods of model parameter calibration need to be created and evaluated. Another possible limitation from the EMG and neural synergy constraint is tha t neither method is able to constrain the magnitudes of excitations for muscles without EMG. For muscles with EMG, there is also the need to determine how appropriate the scaling of EMG signals to the EMG from maximum voluntary contraction trials is for ac curate estimation of muscle excitation magnitudes. Wh ile the EMG based optimization methods presented here allowed for minimized differences in both the shape and magnitude of the muscle excitation signal relative to the scaled EMG signal, the necessity to do so should also be tested in future research. These future research goals are necessary for the advancement of musculoskeletal modelling and need to be evaluated with highly subject specific musculoskeletal models, such as the one presented here. It i s also necessary that model evaluation methods include some means of validation based on experimental measurements, such as contact forces at any particular joint. To accomplish this goal, implantation of instrumented joint replacements in more subjects, e specially subjects who have a high level of ambulation, is also suggested.
96 APPENDI X A DISCRET IZATION OF ACTIVATIO N AND CONTR ACTION DYNAMICS Background T o achieve the most physiological and accurate force predictions possible w hen predicting muscle forces during human movement using inverse rather than forward dynamic optimization method s (i.e., static optimization), it may be necessary to model the full dynamic properties of musculotendon actuators. To achieve this goal, we seek to model activation and contraction dynamics in a way that is numerically accurate, computationally efficient, and easy to implement. While the first order ODEs defining activation dynamics (converts measured muscle excitations into calculated muscle activations) and contraction dynamics (converts calculated muscle activations into calculated muscle forces) could be numerically integrated (e.g., using ODE45 in MATLAB ), this solution approach is computationally costly and subject to drift and stability problems. Discretizing the relevant ODEs over all time frames via finite difference approximations provides an alt ernate solution approach that can achieve a better combination of accuracy, efficiency, and simplicity. Activation dynamics can be defined by a linear first order ODE where the first time derivative of activation can be discretized using a finite differenc e approximation. Applying the discretized ODE to all time frames produces a linear algebraic system of equations in the unknown activations, which can be found using basic linear algebra methods. In contrast, contraction dynamics can be defined by a nonlin ear first order ODE where the first time derivative of tendon force can be discretized using a finite difference approximation. In this case, applying the discretized ODE to all time frames produces a nonlinear algebraic system of equations
97 in the unknown tendon forces, which can be found using a Newton Raphson nonlinear root finding method. To achieve computational efficiency in the iterative root finding process, we derive all contraction dynamics expressions in closed form so that an analytic Jacobian (i .e., first derivative) matrix can be formulated for use by the root finding algorithm. Methods Activation Dynamics We first consider the first order ordinary differential equation that relates musculotendon excitatio n t o mus culo tendon activation where is typically calculated from an experimental electromyography signal. He et al. (1991) and van den Bogert et al ( 2011) proposed the following first order ODE to describe activat ion dynamics ( A 1 ) In this equation, is an activation rate constant, is a deactivation rate constant, and is a time varying quantity defined to simplify the subsequent derivation. Given activation and deactivation time constants and respectively, the following two recursive equations are used to find and :
98 Activation Dynamics Numerical Solution If we disc retize equation ( A 1 ) using a central difference approximation for we obtain ( A 2 ) For the first node we impose the assumption that the first point of the activation is equal to the excitation. This requires that the starting first node is at a time well ahead of the desired first time point for errors to dissipate ( A 3 ) At the l ast node, we employ the use of a ba ckward difference approximation: and apply this equation to the last node to yield or
99 ( A 4 ) Applying equations ( A 2 ) for all intermediate time points along with ( A 3 ) and ( A 4 ) or the first and last time points respectively, a linear algebraic system of equations can be formed and solved for the unknown activations Thus, these equations change a discretized muscle excitat ion time history into a discretized muscle activation time history Contraction Dynamics We next consider the first order ordinary differential equation that relates activation, muscle tendon length, a nd muscle tendon velocity to tendon force. The derivation comes directly from Zajac (1989) with additions from Anderson (2007). A variable pennation angle will be used along with an optimal muscle fiber length that is dependent on the level of activation. Optimal muscle fiber length, can be calculated using the following relationship from Lloyd and Besier (2003): where Thus, the current value of is a function of a constant value and the time varying value of activation For a simple linear model of tendon compliance, tendon force is defined as ( A 5 ) with
100 In these equations, is tendon stiffness, is tendon length, is tendon slack length, and is muscle peak isometric force at full activation. Since we do not know (and have no easy way to find) we cannot solve equation ( A 5 ) directly for tendon force Thus, we consider instead the first time derivative of tendon force, which can be found by taking th e first time derivat ive of equation ( A 5 ) : To formulate the first order ODE for tendon force, the goal therefore becomes to develop an expression for tendon velocity in terms of known quantiti es. The desired expression for can be derived by starting with the relationship between muscle tendon length muscle length and tendon length incor porating a variable pennation angle : ( A 6 ) To derive an expression for we take the first time derivative of equation ( A 6 ) noting that changes as a function of time: ( A 7 ) or
101 Since musculotendon velocity can be calculate d directly from a geometric musculoskeletal model, determination of requires deriving expressions for muscle length muscle velocity pennation angle and the first time derivative of pennation angle in terms of known quantities. Expressions for and can be derived by making the assumption that muscle volume is constant. F or this assumption, muscle width w is also constant and can be defined as a function of pennation angle when the muscle length is equal to the optimal fiber length : At any other muscle fiber length pennation angle is therefore related to muscle width w via ( A 8 ) Since w is constant, can be found from ( A 9 ) assuming muscle fiber length is known. Furthermore, differentiation of e quation ( A 8 ) with respect to time yields which can be solved for assuming and are known: ( A 10 )
102 Incorporation of equation ( A 10 ) into equation ( A 7 ) yields the following simplified expression for : ( A 11 ) Note that muscle fiber length has dropped out of this equation but is still needed to find from equation ( A 9 ) Since musculotendon velocity can be calculated directly from a geome tric musculoskeletal model, the problem of finding is reduced to that of finding and To derive an expression for we begin by making use of the forc e velocity prop erty of muscle. Anderson (2007) and Thelen ( 2003) provide the following relationship for normalized muscle force as a function of normalized muscle lengthening (rather than shortening) velocity : ( A 12 ) In this equation, and where is chosen to be positive when the muscle is lengthening rather than shortening. With this choice, muscle velocity is the first time derivative of muscle length rather than the negative of this der ivative. E quation ( A 12 ) can then be inverted analytically to define as a function of and with the use of the definition of normalized muscle velocity, ( ) the muscle velocity can be calculated However, to avoid numerical issues with the
103 inversion of this equation, a polynomial equation was fit to this inverted equation and used: ( A 13 ) To find we consider how contributes to the in t element, ( A 14 ) which expresses the relationship between contractile element force and muscle activation muscle peak isometric force normalized muscle force leng th properties and normalized muscle force velocity properties Solving equation ( A 14 ) for yields ( A 15 ) Equation ( A 15 ) necessitates finding two new quantities, and where both quantities will depend on To find we first note that tendon force is related to muscle force and penn ation angle via the following geometric relationship:
104 Furthermore, can be decomposed contractile element force and passive element force : We express (i.e., muscle passive force length properties) using the analytic expression ( A 16 ) The normalized muscle length is defined as the muscle length divided by the optimal muscle fiber length By making smooth and continuous over the statements allowed) for the derivative can be derived in the discretized final ODE w ith respect to as needed to facilitate formulation of a fast and accurate root finding solution for over all time frames. Combining equations ( A 5 ) through ( A 16 ) yields the desired expression for : To find (i.e., normalized muscle active force length properties), we use the following Gaussian like analytic function prop osed by Thelen (2003):
105 The next quantity we need to find is which is required to calculate as well as Re arranging equation ( A 6 ) above, we have ( A 17 ) where can be calculated directly from a geometric musculoskeletal model. Furthermore, can be defined as a function of based on equation ( A 5 ) : Subs tituting equations ( A 11 ) into equation ( A 17 ) yields the following expression : ( A 18 ) The final issue to resolve is that we e nd up with two equations, equations ( A 9 ) and ( A 18 ) containing and with appearing nonlinearly in both equations. To solve these two equations for and we re write both equations in the following form: Squaring both sides of both equations and adding the resulting equations together eliminates and produces the following equation for :
106 equation ( A 9 ) can then be used to solve for We are now in a position to define the precise sequence of calculations needed to formulate a nonlinear first order ODE defining how tendon force varies as a function of time t Given a , and we perform our calculations in the following order:
107 Combining these equations may yield a closed form expression for as a function of The resulting first order ODE is nonlinear in our dependent variable Thus, if we discretize it using a finite difference approximation for we will need to employ a nonlinear root finding method for multiple equations (such as Newton Raphson) to solve for the time history of However, when tendon is modeled as rigid, the nonlinear first order ODE for contraction dynamics reduces to a single algebraic equation that can be solved explicitly at ea ch time frame for tendon force. Contraction Dynamics Numerical Solution In simplified form, for a musculotendon model with a compliant tendo n, our final equation for can be written as where contains all of the nonlinear terms in Similar to our first order ODE for activation dynamics this first order ODE can be discretized using a central difference approximation for : ( A 19 ) Assuming slope continuity at the start and end of the cycle the equation for the first node becomes
108 while the equation for the last node is ( A 20 ) Note that each of the equations above is in the form of a nonlinear zero equation as r equired for root finding. Thus, we end up with a system of n nonlinear equations in the n unknowns through To solve these equations for the unknown values of across all time frames, we first define a vector of unknowns x where If is the current guess for the unknown values of across all time frames, then the new guess can be found from the Newton Raphson root finding method for a nonlinear system of equations: In this equation, the vector (a s defined in equations ( A 19 ) through ( A 20 ) ) and the square Jacobian matrix (as defined by taking the derivative of with respect to ) are formulated as follows:
109 and At the solution, we will have to within the specified tolerance. The condition number of should be checked to verify that the matrix is wel l conditioned. To calculate the Jacobian matrix, symbolic manipulation needs to be performed to calculate the analytic derivatives of the nonlinear function of tendon force (derivative of with respect to
110 Results and Discussion To analyze the output of this model, a comparison of the discretized form of the contraction dynamics equations to the forward integrated form was performed for forty muscles of the lower limbs. Guessed excitations along with cal culated musculotendon lengths and velocities from an inverse kinematics solution using a musculotendon model of the leg during a normal walking cycle were used as inputs. Parameter values for both the activation dynamics and musculotendon contraction dynam ics were taken from literature. The average time taken for calculation was 3.2 s for the integration method and 0.02 s for the discretized method. This time reduction of 99.4% is a significant step towards creating a more computationally efficient muscul otendon model. The calculated muscle forces were very similar to each other with an average R 2 value of 0.98 and an average RMS difference of 7 N (relative to a maximum force of ~1000 N).
111 APPENDIX B EMG TO FORCE MODELLI NG Motivation Constraining muscle activations or excitations based on experimental measurement may or may not be appropriate for obtaining unique muscle force estimates. While many methods seek a solution that minimizes the sum of muscle activations or stresses raised to some power at any time instance, this optimality assumption may not be suitable for subjects with impaired neural control systems. EMG driven models have sought to constrain muscle excitations directly by utilizing processed EMG signals. These models require an EMG to force model to account for the discrepancies found between experimental measurements of muscle electrical activity and measured net joint moments. However, due to the low fidelity of some EMG signals and the uncertainty of the reliability of the EMG to force re lationship, a greater level of freedom for the muscle excitation patterns may be needed. The need to adjust both the maximum excitation value and the muscle strength is especially necessary with the inclusion of the contraction dynamics models. These model s create a nonlinear time varying affect between changes in activation and changes in maximum muscle force. A scaling of the muscle activation by a certain amount will not result in the same forces as a scaling of the muscle strength by the same amount. T o maintain dynamic consistency between a musculoskeletal model and its experimental inputs, researchers typically compare muscle and other model contributions to the net loads at any joint with the net loads c alculated from inverse dynamics some model that describes different net loads and contributions. A total of six
112 net inverse dynamic loads can be calculated at any joint but it is common practice to consider only degrees of freedom that undergo large motions. While many studies aim to match only the sag ittal plane net loads at the hip, knee, and ankle, some also include the IE rotation and AA moments at the hip and the inversion/eversion moment at the ankle. The FE moment is typically the only load considered at the knee since FE is the primary motion at the knee and muscles are assumed to be the primary contributors to this moment. Other loads, such as contact and ligament loads, are assumed to contribute significantly to the other five net loads at the knee joint. If contact loads are known, then the ne t SI force along with the net AA moment at the knee may also be used as constraints to increase the uniqueness of muscle force estimates. Optimization and Calibration The trust region reflective algorithm for nonlinear least squares root finding in MATLAB was used and consistently produced the best and quickest results when compared to other optimization techniques including sequential quadratic programming, interior point, and particle swarm methods. This optimization routine allows for a large number of design variables and cost function errors as well as lower and upper bounds on the design variables. Weights are placed on each error term and are picked based on the desired magnitude of error for each error term. Weights are picked such that magnitudes o f errors above the cutoff will be greater than one and errors below the cutoff will be less than one. The weighted error is th e n raised to a power greater than two depending on the desired level of constraint for the particular error. For errors that are d esired to be constrained at a certain level, the power that the weighted error is r a ised to can be increased. For example, to constrain excitations to be below one, a weight of one is placed on the modeled excitations and this error is raised to a large
113 po wer in order to increase the amount of error for any excitations above one and to decrease the error on any excitations below one. This control of the level of constraint can be applied to each error taking into account that each error is also squared in t he non linear least squares optimization algorithm. Numerical gradients were calculated using a user defined central finite difference algorithm. For faster calculations of these gradients, model calculations that are known not to change with the modificat ion of a single design variable were not re calculated and instead the cost function would use variables from the nominal cost function evaluation for that specific calculation. For example, if muscle strength is modified for one muscle, then all calculati ons for the other muscles are not necessary along with any calculations that precede the equations that involve the muscle strength. This can reduce the time for gradient calculation by 90%. Parallel processing of the design gradient calculation was also u sed to speed up these calculations by splitting up multiple design variable modifications into multiple local servers. To avoid issues with gradient based optimization rout ines being held at local minima, optimization s were restarted once for each step of the overall calibration routine This helped to ensure that the solution was not held in a local minima based on numerical gradients. Calibration of all model parameters and muscle excitations together created ill conditioned numerical gradient matrices a nd greatly increased the time taken per optimization step. To avoid this issue several modifications to the optimization scheme were made. The first step was to scale the design variables to ensure consistency of gradient magnitude between parameters and t o be certain that all design variables were between zero and one. P arameter calibration involved
114 optimization of a large number of unknown design parameters, and initial optimization test s with all model parameters included as design variables simultaneous ly failed to converge well This was avoided by careful manipulation of the initial guesses for model parameters and muscle excitations, along with an optimization scheme that gradually changes the design variables Initial optimizations for each model tes ting optimized all design variables simultaneously, but also included error terms representing changes in the design variables from the initial guess. After a limited amount of optimization iterations, the optimization was stopped and a calibration and adj ustment process was begun. The calibration and adjustment process began by calibrating each parameter separately with optimization The minimization of changes in design variables was also performed for these individual parameter optimizations to avoid lar ge changes in parameters which could be deleterious to the overall convergence of the optimization problem by setting the solution in a bad local minimum This is similar to the o ptimization method of simulated annealing with slow cooling r ates. After indi vidual parameter calibrations, the first study performed optimizations of each combination of two model parameters, while in the second study this was replaced by every combination of one model parameter with the parameters to adjust excitation shapes. Aft er these combined parameter calibrations a sequence of adjustments through optimization was performed twice sequentially which i ncluded, in order, the modification of 1) all model parameters simultaneously 2) muscle contributions to net loads and 3) mus cle excitation patterns and gains. This entire calibration and adjustment process was then repeated until iterations did not decrease the overall net loads by more than one percent. When parameters did not calibration, and only excitation patterns were
115 des ired to be modified, a sequence of optimization was performed including modification of, in order, 1) excitation gains, 2) excitation patterns, and 3) excitation gains and patterns. Initial tests were performed by optimization s that did not track EMG signa ls, and did not include any EMG to for ce models. This solution created a best fit set of muscle forces that are desired to be output from the EMG to force model. It also serves as a testing platform to ensure that the musculoskeletal model was properly bui lt and scaled to match the experimental data. The resulting muscle activations were similar to the EMG signals and thus provided good initial guesses of muscle excitation patterns for muscles without EMG. With the best guesses of muscle forces for EMG to force modelling muscle excitations were then reset to the experimental EMG signals and the time delay model between excitations and activations can be used. Initial guesses for time delays and time constants were chosen based on the time between the peaks of the optimized act ivations and the measured EMG. Tracking of the EMG was then introduced to the cost function terms, and the complexity of the model was built using the delays and/or activation dynamics model nonlinear model and either a rigid tendon or compliant tendon muscle model.
116 LIST OF REFERENCES Allen, J.L., Neptune, R.R., 2012. Three dimensional modular control of human walking. Journal of Biomechanics 45, 2157 2163. Amarantini, D., Martin, L., 2004. A method to combine numerical optimizatio n and EMG data for the estimation of joint moments under dynamic conditions. Journal of Biomechanics 37, 1393 1404. Anderson, F.C. Pandy, M.G. 2001 a Dynamic optimization of human walking. Journal of Biomech anical Eng ineering 123, 381 390. Anderson, F.C. Pandy, M.G. 2001 b Static and dynamic optimization solutions for gait are practically equivalent. Journal of Biomech ics 34, 153 161. Anderson, F. C. 2007. BioE215 Physics based Simulation of Biological Structures: Equations for Modeling the Forces Gener ated by Muscles and Tendons. Web based document. Andriacchi, T.P. 1994. Dynamics of knee malalignment. The Ortho paedic Clin ics of N orth A merica 25, 395 403. Arnold, A.S., Delp, S.L., 2005. Computer modeling of gait abnormalities in cerebral palsy: applica tion to treatment planning. Theoretical Issues in Ergonomics Science 6, 305 312. Arnold, E.M., Ward, S.R., Lieber, R.L., Delp, S.L., 2010. A model of the lower limb for analysis of human movement. Ann als of Biomedical Engineering 38, 269 279. Barrett, R.S. Besier, T.F., Lloyd, D.G., 2007. Individual muscle contributions to the swing phase of gait: An EMG based forward dynamics modelling approach. Simulation Modelling Practice and Theory 15, 1146 1155. Besier, T.F., Fredericson, M., Gold, G.E., Beaupr, G.S ., Delp, S.L., 2009. Knee Muscle Forces during Walking and Running in Patellofemoral Pain Patients and Pain Free Controls. Journal of Biomechanics 42, 898 905. Bogey, R.A., Perry, J., Gitter, A.J., 2005. An EMG to force processing approach for determining ankle muscle forces during normal human gait. IEEE Trans actions on Neural Syst ems and Rehabil itation Engineering 13, 302 310. Buchanan, T.S., Delp, S.L., Solbeck, J.A., 1998. Muscular resistance to varus and valgus loads at the elbow. Journal of Biomechani cal Engineering 120, 634 639. Buchanan, T.S., Lloyd, D.G., Manal, K., Besier, T.F., 2005. Estimation of muscle forces and joint moments using a forward inverse dynamics model. Med icine and Sci ence in Sports and Exerc ise 37, 1911 1916. Byadarhaly, K.V., Per door, M.C., Minai, A.A., 2012. A modular neural model of motor synergies. Neural Netw orks 32, 96 108.
117 Cappellini, G., Ivanenko, Y.P., Poppele, R.E., Lacquaniti, F., 2006. Motor patterns in human walking and running. Journal of Neurophysiology 95, 3426 3437 Corcos, D.M., Gottlieb, G.L., Latash, M.L., Almeida, G.L., Agarwal, G.C., 1992. Electromechanical delay: An experimental artifact. Journal of Electromyography and Kinesiology 2, 59 68. Davy, D.T., Audu, M.L., 1987. A dynamic optimization technique for pr edicting muscle forces in the swing phase of gait. Journal of Biomechanics 20, 187 201. Delp, S.L., Anderson, F.C., Arnold, A.S., Loan, P., Habib, A., John, C.T., Guendelman, E., Thelen, D.G., 2007. OpenSim: open source software to create and analyze dynam ic simulations of movement. IEEE Transactions on Biomedical Engineering 54, 1940 1950. De Luca, C.J., 1997. The use of surface electromyography in biomechanics. Journal of A pplied Biomechanics 13, 135 163. Dimitrijevic, M.R., Gerasimenko, Y., Pinter, M.M., 1998. Evidence for a spinal central pattern generator in humans. Ann als of the New York Academy of Sciences 860, 360 376. Disselhorst Klug, C., Schmitz Rode, T., Rau, G., 2009. Surface electromyography and muscle force: limits in sEMG force relationship a nd new approaches for applications. Clin ical Biomech anics (Bristol, Avon) 24, 225 235. D'Lima, D.D., Townsend, C.P., Arms, S.W., Morris, B.A., Colwell, C.W. Jr. 2005 An implantable telemetry device to measure intra articular tibial forces. Journal of Bio mech anics 38, 299 304. Doorenbosch, C.A.M., Harlaar, J., 2003. A clinically applicable EMG force model to quantify active stabilization of the knee after a lesion of the anterior cruciate ligament. Clin ical Biomech anics (Bristol, Avon) 18, 142 149. Erdemir A., McLean, S., Herzog, W., van den Bogert, A.J., 2007. Model based estimation of muscle forces exerted during movements. Clin ical Biomech anics (Bristol, Avon) 22, 131 154. Lima, D.D., 2012. Grand challenge competition to predict in vivo knee loads. Journal of Orthopaedic Research 30, 503 513. Fregly B.J., D'Lima D.D., Colwell C.W. Jr. 2009 Effective gait patterns for offloading the medial compartment of the knee. Journ al of Orthopaedic Res earch 27, 1016 1021. Fregly, B.J., Reinbolt, J.A., Rooney, K.L., Mitchell, K.H., Chmielewski, T.L., 2007. Design of patient specific gait modifications for knee osteoarthritis rehabilitation. IEEE Trans actions on Biomedical Engineering 54, 1687 1695.
118 Garner, B.A., Pandy, M.G., 2003. Estimation of musculotendon properties in the human upper limb. Annals of Biomedical Engineering 31, 207 220. Gerus, P., Rao, G., Buchanan, T.S., Berton, E., 2010. A clinically applicable model to estimate t he opposing muscle groups contributions to isometric and dynamic tasks. Ann als of Biomedical Engineering 38, 2406 2417. Granata, K.P., Marras, W.S., 1995. An EMG assisted model of trunk loading during free dynamic lifting. Journal of Biomechanics 28, 1309 1317. Hayashibe, M., Guiraud, D., Poignet, P., 2009. EMG based neuromuscular modeling with full physiological dynamics and its comparison with modified Hill model. Engineering in Medicine and Biology Society Annual International Conf erence of the IEEE 200 9 6530 6533. He, J., Levine, W.S., Loeb, G.E., 1991. Feedback gains for correcting small perturbations to standing posture. IEEE Transactions on Automatic Control 36, 322 332. Heine, R., Manal, K., Buchanan, T.S., 2003. Using Hill Type Muscle Models and E MG Data in a Forward Dynamic Analysis of Joint Moment. Journal of Mechanics in Medicine and Biology 03, 169 186. Heintz, S., Gutierrez Farewik, E.M., 2007. Static optimization of muscle forces during gait in comparison to EMG to force processing approach. Gait and Posture 26, 279 288. Helmick, C.G., Felson, D.T., Lawrence, R.C. 2008. Estimates of the prevalence of arthritis and other rheumatic conditi ons in the United States: Part 2. Arthritis and Rheumatism 58 15 25. Higginson, J.S., Zajac, F.E., Neptune R.R., Kautz, S.A., Delp, S.L., 2006. Muscle contributions to support during gait in an individual with post stroke hemiparesis. Journal of Biomechanics 39, 1769 1777. Hof, A.L., Pronk, C.N., van Best, J.A., 1987. Comparison between EMG to force processin g and kinetic analysis for the calf muscle moment in walking and stepping. Journal of Biomechanics 20, 167 178. Hof, A.L., Van den Berg, J., 1981. EMG to force processing I: An electrical analogue of the Hill muscle model. Journal of Biomechanics 14, 747 7 58. Holzbaur, K.R. Murray, W.M., Gold, G.E., Delp, S.L., 2007. Upper limb muscle volumes in adult subjects. Journal of Biomechanics 40, 742 749. Ivanenko, Y.P., Poppele, R.E., Lacquaniti, F., 2004. Five basic muscle activation patterns account for muscle activity during human locomotion. The Journal of Physiol ogy 556, 267 282.
119 Jacobs, R., Bobbert, M.F., van Ingen Schenau, G.J., 1996. Mechanical output from individual muscles during explosive leg extensions: the role of biarticular muscles. Journal of Biome chanics 29, 513 523. Jia, B., Kim, S., Nussbaum, M.A., 2011. An EMG based model to estimate lumbar muscle forces and spinal loads during complex, high effort tasks: Development and application to residential construction using prefabricated walls. Internat ional Journal of Industrial Ergonomics 41, 437 446. Jonkers, I., Spaepen, A., Papaioannou, G., Stewart, C., 2002. An EMG based, muscle driven forward simulation of single support phase of gait. Journal of Biomechanics 35, 609 619. Kirking, B., Krevolin, J. multiaxial force sensing implantable tibial prosthesis. Journal of Biomechanics 39, 1744 1751. Kumar, D., Rudolph, K.S., Manal, K.T., 2012. EMG driven modeling approach to muscle force and joint load estimations: case study in knee osteoarthritis. Journal Orthopaedic Research 30, 377 383. Kutch, J.J., Valero Cuevas, F.J., 2012. Challenges and new approaches to proving the existence of muscle synergies of neural origin. PLoS Comput. Biol. 8, e1002434. Lee, D.D., Seung, H.S., 1999. Learning the parts of objects by non negative matrix factorization. Nature 401, 788 791. Lin, Y.C., Walter, J.P., Banks, S.A., Pandy, M.G., Fregly, B.J., 2010. Simultaneous prediction of muscle and contact forces in the knee d uring gait. Journal of Biomechanics 43, 945 952. Lloyd, D.G., Besier, T.F. 2003 An EMG driven musculoskeletal model to estimate muscle forces and knee joint moments in vivo. Journal of Biomech ics 36, 765 776. MacKay Lyons, M., 2002. Central pattern gener ation of locomotion: a review of the evidence. Phys ical Ther apy 82, 69 83. Manal, K., Buchanan, T.S., 2003. A one parameter neural activation to muscle activation model: estimating isometric joint moments from electromyograms. Journal of Biomechanics 36, 1 197 1202. McGill, S.M., 1992. A myoelectrically based dynamic three dimensional model to predict loads on lumbar spine tissues during lateral bending. Journal of Biomechanics 25, 395 414. Miyazaki, T., Wada, M., Kawahara, H. 2002 Dynamic load at baseline can predict radiographic disease progression in medial compartment knee osteoarthritis. Annals of Rheumatic Diseases 61, 617 622.
120 Muramatsu, T., Muraoka, T., Takeshita, D., Kawakami, Y., Hirano, Y., Fukunaga, T., 2001. Mechanical properties of tendon and aponeurosis of human gastrocnemius muscle in vivo. Journal of Applied Physiology 90, 1671 1678. Neptune, R.R., Clark, D.J., Kautz, S.A., 2009. Modular Control of Human Walking: A Simulation Study. Journal of Biomechanics 42, 1282 1287. Neptune, R.R., Hull, M.L., 1998. Evaluation of performance criteria for simulation of submaximal steady state cycling using a forward dynamic model. Journal of Biomechanical Engineering 120, 334 341. Neptune, R.R., Zajac, F.E., Kautz, S.A., 2004. Muscle mechanical work requir ements of mass is significant. Journal of Biomechanics 37, 817 825. Nussbaum, M.A., Chaffin, D.B., 1998. Lumbar muscle force estimation using a subject invariant 5 parameter EMG based m odel. Journal of Biomechanics 31, 667 672. Piazza, S.J., Delp, S.L., 1996. The influence of muscles on knee flexion during the swing phase of gait. Journal of Biomechanics 29, 723 733. Reinbolt, J.A., Schutte, J.F., Fregly, B.J., Koh, B.I., Haftka, R.T., G eorge, A.D., Mitchell, K.H., 2005. Determination of patient specific multi joint kinematic models through two level optimization. Journal of Biomechanics 38, 621 626. Schache, A.G., Fregly, B.J., Crossley, K.M., Hinman, R.S., Pandy, M.G. 2008 The effect of gait modification on the external knee adduction moment is reference frame dependent. Clinical Biomechanics (Bristol, Avon) 23, 601 608. Schipplein, O.D., Andriacchi, T.P. 1991 Interaction Between Active and Passive Knee Stabilizers. Journal of Orthop aedic Res earch 113 119. Shao, Q., Bassett, D.N., Manal, K., Buchanan, T.S., 2009. An EMG driven model to estimate muscle forces and joint moments in stroke patients. Computers in biology and medicine 39, 1083 1088 Sharma, L., Hurwitz, D.E., Thonar, E.J. 1998 Knee adduction moment, serum hyaluronan level, and disease severity in medial tibiofemoral osteoarthritis. Arthritis and Rheumatism 41, 1233 1240. Silder, A., Whittington, B., Heiderscheit, B., Thelen, D.G., 2007. Identification of Passive Elastic Joint Moment Angle Relationships in the Lower Extremity. Journal of Biomechanics 40, 2628 2635. Taga, G., 1995. A model of the neuro musculo skeletal system for human locomotion. I. Emergence of basic gait. Biol ogical Cybern etics 73, 97 111. Thelen D. G. 2003 Adjustment of muscle mechanics model parameters to simulate dynamic contractions in older adults. J ournal of Biomech anical Eng ineering 125, 70 77.
121 Thorp, L. E. Sumner, D.R. Block, J.A. Moisio, K.C. Shott, S. Wimmer, M.A. 2006 Knee joint loading differs in individuals with mild compared with moderate medial knee osteoarthritis. Arthritis and Rheumatism 54, 3842 3849. Thorp, L.E., Sumner, D.R., Wimmer, M.A., Block, J.A. 2007. Relationship Between Pain and Medial Knee Joint Loading in Mild Radiogr aphic Knee Osteoarthritis. Arthritis and Rheumatism 57, 1254 1260. Ting, L.H., Macpherson, J.M., 2005. A limited set of muscle synergies for force control during a postural task. Journal of Neurophysiology 93, 609 613. van den Bogert, A.J., Blana, D., Hein rich, D., 2011. Implicit methods for efficient musculoskeletal simulation and optimal control. Procedia IUTAM 2, 297 316. Walter J. D. D, Colwell C. W Jr Fregly B. J. 2010 Decreased knee adduction moment does not guarantee decreased medial c ontact force during gait. J ournal of Orthop aedic Res earch 28 1348 1354. Ward, S.R., Eng, C.M., Smallwood, L.H., Lieber, R.L., 2009. Are current measurements of lower extremity muscle architecture accurate? Clin ical Orthop aedic s and Related Research 467, 1074 1082. White, S.C., Winter, D.A., 1992. Predicting muscle forces in gait from EMG signals and musculotendon kinematics. Journal of Electromyogr aphy and Kinesiol ogy 2, 217 231. Winters, J.M., Stark, L., 1988. Estimated mechanical properties of synergist ic muscles involved in movements of a variety of human joints. Journal of Biomechanics 21, 1027 1041. Woltring, H.J., 1985. On optimal smoothing and derivative estimation from noisy displacement data in Biomechanics Human Movement Science 4, 229 245. Zaja c, F.E., 1989. Muscle and tendon: properties, models, scaling, and application to Biomechanics and motor control. Crit ical Rev iews in Biomedical Engineering 17, 359 411. dial and lateral tibial loads during dynamic and high flexion activities. Journal Orthopaedic Research 25, 593 602. 2007. Correlation between the Knee Adduction Torque and Medial Contact Force for a Vari ety of Gait Patterns. Journal of Orthopaedic R esearch 25, 789 797. Zhao, D., Sakoda, H., Sawyer, W.G., Banks, S.A., Fregly, B.J., 2008. Predicting knee replacement damage in a simulator machine using a computational model with a consistent wear factor. Jou rnal of Biomechanical Engineering 130, 011004.
122 BIOGRAPHICAL SKETCH Jonathan Walter received his May, 2007. He majored in a erospace e ngineering with an emphasis in numerical methods He received his in m echanical e ngineering from the University of Florida in 2011 with an emphasis in rigid body dynamics applied to biomechanical analyses He received h is doctoral degree in mechanical engineering from the Department of Mechanical and Aer ospace Engineering at the University of Florida in December of 2012