Muscle activation underestimated

hessers · October 15, 2015, 3:15pm

Hello,

Attached model is a modification from the MoCapModel; FullBody example in version 6.0.5.(4379). Currently I’m using the model to estimate joint loading and muscle activity during tasks involving movement of the right arm.

I’ve compared the maximum estimated muscle activity to the maximum measured EMG and the model seems to underestimate the needed activity heavily (results are attached as PDF). The reason why seems unclear to me at the moment, thus I would like to have some professional help. As a remark, the motion the simulation model performs is not accurate, as the pelvis moves from the seat which did not occur in the experiments. Therefor it needs more constraints in the vertical direction and would like to know how this should be implemented.

Additional info:
Included in the Input folder you will find a Reaching Far Upper Left (RFUL, known as Contralateral), Right (RFUR, known as Ipsilateral) and Shoulder Ab- Adduction (SAA) movement recorded to .c3d file. Also, I’ve attached data analysis of the estimated joint loading and muscle activity (MA) and measured EMG (ME) in these three movements (C1_RFUR.mat, C1_RFUL.mat, C1_SAA.mat) as ‘C1_EMG.zip’ in the first folder. Six muscles are included; Biceps, Triceps, Deltoid, Trapezius, Pectoralis major and Latissimus Dorsi. Also, the GlenoHumeral Ab-Adduction, Flex-Extension and Rotation, and Elbow Flex-Extension can be found in the MA structure. Subject C1 equals to row 12 in this structure. When opening for example the Biceps, the structure shows ‘mp’, ‘pa’ and ‘mt’. The first parameter ‘mp’ shows the averaged activity of all Muscle Parts over time which correspond in the model to the location of the measured EMG. ‘pa’ means separate muscle Parts Activity and ‘mt’ is the Maximum activity of ‘mp’ and its corresponding Time. This last parameter ‘mt’ is used to compare to the measured EMG (ME), which has the same structure. Thus, I compare for example the maximum muscle activity MA.Bic.mt (first value) to ME.Bic.mt (first value) in all conditions. If that seems to be correct I also check the averaged activity pattern per muscle using ‘mp’. Individual muscle parts ‘pa’ has not been used yet. Multiple subjects have been investigated leading to the results in the reference (PDF).

Thanks in advance for the help.

Greets, Hans

ananthgk1 · October 20, 2015, 9:34am

Dear Hans,

You can constrain the pelvis by adding the following lines of code to the file “ExtraDrivers.any”. The code completely constrains all degrees of freedom of the pelvis to constant values throughout the movement, so you must remove any DoFs where you still want to retain mobility.

You would have to replace the values written as <posx> or <rotx> with the appropriate values corresponding to your experiment. This would depend on your laboratory coordinate system, the height of the seat etc,. Hope this helps.

Regards
Ananth
AnyBody Support

AnyKinEqSimpleDriver ConstrainPelvis =
{
AnyFolder &TrunkDOFs = Main.Studies.HumanModel.BodyModel.Interface.Trunk;

AnyKinMeasure &ref0 = TrunkDOFs.PelvisPosX;
AnyKinMeasure &ref1 = TrunkDOFs.PelvisPosY;
AnyKinMeasure &ref2 = TrunkDOFs.PelvisPosZ;

AnyKinMeasure &ref3 = TrunkDOFs.PelvisRotX;
AnyKinMeasure &ref4 = TrunkDOFs.PelvisRotY;
AnyKinMeasure &ref5 = TrunkDOFs.PelvisRotZ;

DriverPos = {<posx>,<posy>,<posz>,<rotx>,<roty>,<rotz>};
DriverVel = {0,0,0,0,0,0};

};

hessers · October 21, 2015, 2:58pm

Dear Ananth,

Thank you, this seems to work better than my previous solution, the Pelvis is now constrained to remain seated in medio-lateral and vertical position. I’m running a few tests now with the constrained Pelvis. Did you or someone else from Tech support had a chance to look at the disagreement in muscle activity?

Greets, Hans

ananthgk1 · October 26, 2015, 11:25am

Dear Hans,

I had a look at the PDF you had attached and I understand that ‘Graph 2’ is where you compare the estimated and measured muscle activations.

It seems like the Trapezius, Pectoralis and Latissimus muscles are where the significant differences are seen in both contra and ipsilateral cases. This could be due to two reasons:

How accurate are the MVC EMG measurements for these muscles? How confident can you be that the MVC corresponds to the maximum activity the muscle can generate?
Have you scaled your muscle strengths realistically? The maximum muscle strength in your model might not correspond to the actual patient you were simulating.

Please let me know what you think about these points, and we can further diagnose the issue from there.

Regards
Ananth
AnyBody Support

hessers · October 28, 2015, 9:57am

Dear Ananth,

In Graph 2 we compare the effect of the support device on the measured and estimated muscle activity. In Graph 1, we show the muscle activity while Unsupported and with the use of the Support device. The interesting part now is comparing the Unsupported cases for the measured and estimated muscle activity. For example, we expected to see activity in the Trapezius during both tasks, but can only see it in the measured one.

In this study only healthy participants have participated. We checked the measured muscle activity visually during the MVC (before the experiments) and activity from the experiments during data analysis, and we are confident that the measured activity is valid (at least in most cases). We did scale according to height and weight of the participant.
We have not modified the muscle strength in the model because only healthy participants were included. We also did a test where the model had 25% muscle strength, however the muscle activities did not show a different activation pattern, but just increased 4 times. We’ve used a 3rd order polynomial recruitment algorithm with the 3 element Hill muscle model. The minmax function also showed very little change in amplitude and activation patterns. Seeing there were no similarities between measured and estimated muscle activity across all participants (exception for Deltoid during Ipsilateral), I figured it might be a recruitment mismatch. The 3rd order polynomial function was chosen based on a recommendation from an Anybody article (2006).

Greets, Hans

hessers · October 29, 2015, 5:09pm

Dear Ananth,

I’ve looked at the pelvis constraints and I’ve noticed that these are very strict/hard. I’ve implemented constraints so the pelvis is unable to translate in the medio-lateral and vertical direction. Also, it is not allowed to rotate around the roll axis. The model now has trouble simulating the recorded movements during the MotionOptimization study.

AnyKinEqSimpleDriver ConstrainPelvis =
{
AnyFolder &TrunkDOFs = Main.Studies.HumanModel.BodyModel.Interface.Trunk;

AnyKinMeasure &ref0 = TrunkDOFs.PelvisPosX;
//AnyKinMeasure &ref1 = TrunkDOFs.PelvisPosY;
AnyKinMeasure &ref2 = TrunkDOFs.PelvisPosZ;

//AnyKinMeasure &ref3 = TrunkDOFs.PelvisRotX;
AnyKinMeasure &ref4 = TrunkDOFs.PelvisRotY;
//AnyKinMeasure &ref5 = TrunkDOFs.PelvisRotZ;

DriverPos = {(Main.ModelSetup.C3DFileData.Points.Markers.C7.Pos[0][0]+
Main.ModelSetup.C3DFileData.Points.Markers.T10.Pos[0][0]+
Main.ModelSetup.C3DFileData.Points.Markers.CLAV.Pos[0][0]+
Main.ModelSetup.C3DFileData.Points.Markers.STRN.Pos[0][0])/4,
Main.ModelSetup.C3DFileData.Points.Markers.C7.Pos[0][2]-
0.288*Main.BodyH/100 //David Winter Anthropometrics page 52, Shoulder to Pelvis distance = 0.288 * Body Height
,0};
DriverVel = {0,0,0};
};

The performed experiments were done while seated and participants may have had some movement in these directions. How can I implement a more soft version of these constraints to allow the Motion Optimization to change the model’s position and allow for small rotations?

I appreciate the help and therefore I thank you.

Greets, Hans

ananthgk1 · October 30, 2015, 11:22am

Dear Dan,

Regarding the EMG activities, I am not able to completely understand what may be going wrong with the current information. I don’t have a MATLAB license for opening your EMG files. I would hence request the following information:

Magnitude vs time plots for the EMG (normalized such that peak EMG = peak muscle activation) with estimated muscle activations overlaid. Since I’m interested in the temporal variation of both, I am okay with the unrealistic EMG magnitude.
Details of the functional tasks during which your MVC EMG voltage for each muscle was recorded? The third order polynomial recruitment criterion seems like an okay choice as well. Out of curiosity, at what percentage of the original muscle strengths do the peak estimated and measured muscle activation magnitude show some sort of correspondence?

With regards to the pelvis drivers, it is indeed possible to make your drivers soft by adding the following statements in your “ConstrainPelvis” folder:

CType = {Soft,Soft,Soft};

AnyFunConst ConstFun = {Value = {0.01,0.01,0.01}; };
WeightFun = {&ConstFun};

This code sets constant weights for your soft constraint. Remember that during the kinematics operation, the pelvic soft constraints will be fighting against your soft marker drivers, so remember to keep the above mentioned weights in balance with the values in “Markers.any”.

Also remember that your pelvis translations and rotations are measured with respect to the pelvis local coordinates. So it might be a good idea to visualize this before confirming the DriverPos values.

Regards,
Ananth
AnyBody Support

hessers · October 30, 2015, 1:04pm

Dear Ananth,

Please find a zipfile with the requested pictures attached. The Blue line represents the measured EMG. The Red line is the mean estimated muscle activity from the selected muscle parts (Green lines) based on the EMG sensor location. The Black lines are the excluded muscle parts of the same muscle.

I’ve attached the most relevant pictures, thus only Unsupported condition for the Reaching Far Upper Left (RFUL) and Right (RFUR). Also, in subject 13 (S13) the model was both at 100% and 25% (S13C_025) muscle strength.

The MVC was recorded while a researcher restricted the movement of the subject. This seemed to be sufficient with the exception for Trapezius where the researcher was lifted from the ground.

Biceps Brachii; The upper arm is placed relaxed next to the torso. The elbow is flexed 90° and placed in mid position (palm facing up). The elbow is supported and the wrist is fixated.
Instructions: Pull the forearm to the shoulder, while the examiner holds the arm in the same place.

Deltoid: The straight arm is placed in 60° shoulder abduction, with the thumb pointing forward. The arm is fixated just above the elbow joint and the wrist.
Instructions: Push the whole arm sideward up as hard as possible, while the examiner holds the arm in place.

Triceps brachii: Elbow flexed 90°, upper arm abducted. Instructions: Extend the forearm against the resistance of the examiner.

Trapezium Descendens: Elbow flexed 90°, upper arm abducted. Instructions: Pull the shoulder upwards

Pectoralis major: Arm abducted at 90° (horizontal) in the frontal plane and elbow flexed 90° with the hand palm facing down, adduct the arm towards the sagittal plane

Latissimus dorsi: Arm abducted at 90° (horizontal) in the frontal plane and elbow flexed 90° with the hand palm facing down. Instruction: adduct the upper arm towards the trunk

I have not checked when the altered model muscle strength would match the measured EMG. However, when looking at the 25% muscle strength for Latissimus and Pectoralis there were no noted changes in estimated muscle activity. On the other hand for the Trapezius the amplitude were comparable, but the individual muscle parts provided over 100% amplitude even with an Upperbound of 1. Does the Upperbound prevent the muscles from overloading or does it only give a warning? I’ve only seen the latter…

InverseDynamics.Criterion.UpperBoundOnOff=On;
InverseDynamics.Criterion.UpperBound=1.0;
InverseDynamics.Criterion.Type = MR_Polynomial;
InverseDynamics.Criterion.Power = 3;

Greets, Hans

ananthgk1 · November 3, 2015, 8:57am

Hi Hans,

I’m trying to get a good understanding of your problem and will get back to you as soon as possible!

Regards
Ananth

hessers · November 9, 2015, 10:38am

Hello Ananth,

I’ve checked the joint loading during Abduction movement with reference values from OrthoLoad. Attached you will find graphs where the red line represents the AnyBody estimated joint loading and the green and black line are reference values from Orthoload. Data has been normalized to % of bodyweight. Movement has been performed from arm hanging next to body to 150 degrees of abduction and returning to starting position in OrthoLoad experiments and above 150 degrees abduction for our study, but other than that the movement is the same. However, the patterns and magnitude between AnyBody and Orthoload are different. Is it possible that the joint loading is not estimated correctly in my model, which could explain the difference in muscle activity? Perhaps some settings that I’ve missed? I’ve used the MoCap Model from AnyBody demo’s and modified it to be on the upper extremity only. Other than the Pelvis constraint and gravity acting on the Z-axis (I’ve checked it for the correct direction) no other changes were made.

Greets, Hans

ananthgk1 · November 9, 2015, 1:41pm

Dear Hans,

Thank you for detailed plots and I’m sorry I have been slow in going through all of them and getting back. I was looking through your first set of EMG vs. activation plots and was trying to run the simulation for the RFUL movement in AnyBody so that I could get a more physical understanding for why the model was producing the observed muscle activation patterns.

I tried doing this by uncommenting the lines you had for the RFUL movement in your main file, but the parameter identification and motion analysis simulation seemed to fail. Could this be because I do not have the correct version of the code?

It is of course possible that the model is being inaccurate with the predictions, but before getting there I’d like to first understand the simulation properly and see if anything weird is going on with muscle wrapping. etc. Could you please check the code version or let me know if there is other information I need to uncomment/comment before running the RFUL trial?

Regards
Ananth

ananthgk1 · November 9, 2015, 1:46pm

Could you also let me know the name of the subject whose c3d files you included in the shared model?

hessers · November 9, 2015, 1:47pm

#include “…/libdef.any”

// Switch to save all results in output folder
#ifndef AutoSaveOption
#define AutoSaveOption 1
#endif

// Detailed description can be found at:
#include “Model/Description.any”

//Set this to 1 if you want to run the motion and Parameter Optimization identification
//**************************************************
#ifndef MotionAndParameterOptimizationModel
#define MotionAndParameterOptimizationModel 1
#endif

//Set this to 1 if you want to run the inverse dynamic analysis
#ifndef InverseDynamicModel
#define InverseDynamicModel 0
#endif
//Usually only have one of the two switches active so set the inactive analysis to 0
//**************************************************

// Following will include several predefined classes to generate markers and environments
// Usually, there is no need to make changes in this file:
#include “Model/Classes.any”

Main = {

// This file controls most of the Settings of your Model, the following items are defined here:
// 1. Marker Names
// 2. Marker locations

// #include “Model/ModelInput.any”

//RFUL
AnyFloatVar BodyH = 178;
AnyFloatVar BodyM = 78;
AnyFloatVar StFr = 784;
AnyFloatVar EnFr = 1411;
AnyFloatVar SL_CW = 0;
AnyFloatVar G = -9.81;
AnyString FN ="…/Input/RFUL 1";
AnyString FNM ="…/Output/C1_RFUL_Control_InvD_V6_Pelvis";
//
////RFUR
//AnyFloatVar BodyH = 178;
//AnyFloatVar BodyM = 78;
//AnyFloatVar StFr = 716;
//AnyFloatVar EnFr = 1472;
//AnyFloatVar SL_CW = 0;
//AnyFloatVar G = -9.81;
//AnyString FN ="…/Input/RFUR 1";
//AnyString FNM ="…/Output/C1_RFUR_Control_InvD_V6_Pelvis";
//
//SAA
//AnyFloatVar BodyH = 178;
//AnyFloatVar BodyM = 78;
//AnyFloatVar StFr = 1218;
//AnyFloatVar EnFr = 2192;
//AnyFloatVar SL_CW = 0;
//AnyFloatVar G = -9.81;
//AnyString FN ="…/Input/SAA 1";
//AnyString FNM ="…/Output/C1_SAA_Control_InvD_V6_Pelvis";

#include “Model/ModelSetup.any”

// Define the Environment around the human
#include “Model/Environment.any”

// This file controls most of the Trial Specific setup, the following items are defined here:
// 1. C3D file Name
// 2. Start and End Frame
// 3. Subject Anthropometric data
// 4. Initial guess on position
// If you want to run the model on your own data this is the place to start:
#include “Input/TrialSpecificData.any”

AnyFolder Studies ={
// Usually, there is no need to make changes in the following files. PLease be aware. that
// changes in the Kinematics and InverseDynamics may lead to model failure and/or unrealistic results:
#if MotionAndParameterOptimizationModel
#include “Model/Kinematics.any”
#endif

#if InverseDynamicModel
#include "Model/InverseDynamics.any"  
#endif

};
}; //Main

hessers · November 9, 2015, 1:53pm

It seems that the RFUL 1.c3d file from C1 can’t be solved with this model. I’ve attached a new file, please replace the file and text:

//RFUL
AnyFloatVar BodyH = 178;
AnyFloatVar BodyM = 70;
AnyFloatVar StFr = 802;
AnyFloatVar EnFr = 1277;
AnyFloatVar SL_CW = 0;
AnyFloatVar G = -9.81;
AnyString FN ="…/Input/RFUL 1_C3";
AnyString FNM ="…/Output/C3_RFUL_Control_InvD_V6_Pelvis";

It works fine without the name of the subject for me, does it give marker errors for you?

ananthgk1 · November 9, 2015, 1:56pm

Nope, it does not give me errors, but I wanted to know which of the graphs you sent me correspond to this subject alone. Otherwise I was getting a little lost in the data :). I could then use that as a baseline to compare any changes that occur when I alter the model.

hessers · November 9, 2015, 2:00pm

The files I’ve sent are from Subject 12, and the file I’ve sent just now is from Subject 14. I needed to alter the model a bit to get it running again, chest markers should have a weight of 1.0 and then it should work. Sorry for the confusion.

ananthgk1 · November 9, 2015, 2:00pm

This file does not work either and in fact crashes on the first time-step. It is possible that I’m missing something. Could you please re-share your entire AnyBody working folder where RFUL works, and I can use that as a new start point?

ananthgk1 · November 9, 2015, 2:01pm

No worries Hans Just share the new files and confirm which subject the files are for, and I’ll take it from there.

hessers · November 9, 2015, 2:03pm

This works for me, I’ve commented out the Pelvis constraints in ExtraDrivers.any as they were not present in the original debug version.

ananthgk1 · November 12, 2015, 10:12am

Hi Hans,

The parameter identification and motion optimization still fails in the new files that you sent me. In the final iteration which crashed, the failure occured at 4.75 seconds.

Since our debugging focus is more on the inverse dynamics, could you please run the motion optimization and just send me the text files with the joint angle data for the RFUL 1 motion for Subject 12 (I assume this is the subject whom the latest data corresponds to)?

Regards
Ananth
AnyBody Support