Joint Reaction Analysis

[Nicholas Vandenberg]

Hi all,
I've been using Moco, working with simulating rehab in TFA models during stance phases. Previously I'd been looking into how different interventions affected changes in movement patterns, I've recently started looking at how the JRFs change (using code shared by Pasha). I'm wondering if anyone can explain what's goin on "under-the-hood" with this though, particularly 'analyzeSpatialVec'? I see it using the states and controls of the Moco trajectory to find the JRFs, but I'm curious how it does that. In previous OpenSim work I've used the analyze tool to find JRFs with results from a muscle optimization, some kinematics, and an external loads file. I wouldn't expect a 1:1 outcome between Moco and the other workflow but the differences over the same analysis timeframe are a bit alarming.

HJRF.jpg (30.3 KiB) Viewed 2571 times

The JRF output from Moco seems to missing the loading response from early stance, even though the kinematic solution is pretty close to the experimental data that I use as the input kinematics,

Kin.jpg (119.2 KiB) Viewed 2571 times

and the predicted GRFs from the contact spheres also look pretty good.

GRF.jpg (40.21 KiB) Viewed 2571 times

I'm wondering if analyzeSpatialVec is doing some different type of analysis compared with the standard analyze tool and what could be going on to miss the loading during that phase of gait?

[Ross Miller]

Hi Nicholas,

What do the hip muscle forces in early stance for the Moco solution look like?

I have not done hip but I get reasonable-looking knee contact forces in Moco solutions with the method here. Some small changes in syntax may be necessary, this code is from an early version of Moco:

Hi Nick,

Thank you so much, this works! 'flatten()' provides you with a TimeSeriesTable with 6 + 1 columns for each column of TimeSeriesTableSpatialVec. Essentially, this means that each element of a SpatialVec (3 Torque and 3 Force elements), gets unpacked to its "own column', and the extra column is the time column.

As a final question: the default behaviour of "analyzeSpatialVec" is to express the SpatialVecs in the "Ground" (i.e. an inertial) frame of reference, is this correct?

If so, to transform each spatialVec to the desired frame, I must implement something like the following (pseudo) code:

Code: Select all

 state = model.realizeDynamics(statesTraj.get(0)) ; %statesTraj is the StatesTrajectory from a prediction solution

    Torque_ground = Vec3(JR(1,1), JR(1,2), JR(1,3)); 
    Force_ground  = Vec3(JR(1,4), JR(1,5), JR(1,6));
    
    %JR is an array obtained after flattening a TimeSeriesTableSpatialVec and converting 
    % it to a struct using osimTableToStruct
    
    ground = model.getGround();
    child_frame = model.getJointSet().get(0).getChildFrame; %change this to the desired frame for each joint
        
    Torque_local =  ground.expressVectorInAnotherFrame(state, Torque_ground, child_frame);
    Force_local  =  ground.expressVectorInAnotherFrame(state,  Force_ground, child_frame);
        
    joint_reactions_local(1,1) = Torque_local.get(0);
    joint_reactions_local(1,2) = Torque_local.get(1);
    joint_reactions_local(1,3) = Torque_local.get(2);
    
    joint_reactions_local(1,4) = Force_local.get(0);
    joint_reactions_local(1,5) = Force_local.get(1);
    joint_reactions_local(1,6) = Force_local.get(2);

Thank you so much, I think I'm all set to implement my project now!
Pasha

EDIT: different implementation to transform the vectors

Ross

[Nicholas Vandenberg]

Hi Ross,

Some of the stance limb hip abductors and flexors are spiking around heel-strike, a couple are starting at a pretty high activation (~0.8) during the first frame that drops off right after.

The code you shared is what I've been using (or adapted from that same discussion chain). My question is more related to how that is doing the JRF analysis, are the GRF produced by the contact spheres affecting that at all? Or is the "weirdness" in the output JRF stemming from a lack of some initial activation bound in the muscles?

Thanks,
Nick

[Ross Miller]

There is not a lot of documentation on the technical side of JointReaction. You would probably have to find the CPP code for it to really figure out what it is doing. It reads and looks to me like it sums the net joint force with the actuator forces at the joint e.g. muscle forces.

Since your GRF look pretty good I think they or the contact spheres are probably not the source.

Does your simulation have a periodicity constraint? Without that, with at least some of the collocation schemes, the first timestep is unconstrained and you can get strange results there e.g. very-high activations. Can probably also use the various initial-state MocoGoals if periodicity is undesirable.

Ross

[Nicholas Vandenberg]

Thanks, I haven't been able to find much yet, but I'll keep looking into that.

I know I'd had some weirdness with my periodicity constraint previously, because I'm only running sims over the amputated limb stance phase. I'd tried to adapt the periodicity code from the 'example2Dwalker' since that was also over just a stance but had some issues with adapting to 3D as far as state names. I've been hard-coding my periodicity constraint (manually adding each goal pair) but could very well have done it wrong. I haven't really done anything with initial-state MocoGoals yet, could that be a good way to ensure the initial activations in addition to periodicity?

Nick

[Ross Miller]

If you are simulating the stance only (and not a full stride) then the periodicity goal will not be appropriate. I would look into the various goals with Initial in the name.

For limb loss since it is presumably asymmetric, the periodicity over even a full step is not applicable, would need to be a full stride.

Nick B can probably give a more specific recommendation. I have only ever done strides!

Ross

[Nicholas Bianco]

Hi Nick,

If the initial activations are responsible for the large joint loads, you could try adding a MocoInitialActivationGoal to the problem, which constrains the initial activations to be equal to the initial excitations. This way, the initial activations will be appropriately penalized since they are effectively tied to the MocoControlGoal. (MocoInverse problems use this goal by default).

-Nick

[Nicholas Vandenberg]

That's good to know, thank you both! I appreciate the help.

Nick

[Ross Miller]

The "initial states" problem here can be a big problem. The way that the system dynamics are enforced as constraints, at least in most collocation schemes that I know of, are with finite difference calculations that involve the current state and the previous state, or something similar. But for the first timestep there is no previous timestep to use, so the initial state is "unconstrained" in this regard. You can specify a specific value or bounds on that value for those initial states, but even then it gets a little thorny. For example if doing a "walk with minimum effort" problem, the optimizer will generally pick the highest allowed value for the initial horizontal velocity, this is free movement for zero effort.

For motions that are, on average, periodic, like gait, this is where the periodicity constraint is helpful. There is still no previous timestep to constrain the initial states, but they are constrained to be equal to the final states, which do have a previous timestep to constrain their dynamics. Basically it is using the second-to-last timestep to inform the dynamics constraint of the initial timestep, in a roundabout way.

Nick, not sure if you have data (tracking targets) on a full stride, but if you do, you might get better results and simpler problems by simulating a full stride, and can avoid making questionable assumptions about the initial states e.g. is assuming the initial activations are equal to the initial excitations a good assumption? Not necessarily.

Ross

[Nicholas Vandenberg]

Hi Ross,

Sorry for the delayed response, went down a bit of a rabbit hole. I do have data from a full IK trial that I select windows of for the tracking problems, where I've been isolating to a stance phase previously. I'm wondering if I'm setting up my periodicity wrong or if I've got the wrong idea with expanding my time range? I've been hard-coding the periodicity goal (syntax for the lumbopelvic DOFs seems weird but I tried to follow this: https://opensim-org.github.io/opensim-m ... ml#details). If I'm understanding the goal correctly this should be the way to define periodicity?

Code: Select all

periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/groundPelvis/pelvis_tilt/value'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/groundPelvis/pelvis_tilt/speed'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/groundPelvis/pelvis_list/value'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/groundPelvis/pelvis_list/speed'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/groundPelvis/pelvis_rotation/value'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/groundPelvis/pelvis_rotation/speed'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/hip_l/hip_flexion_l/value','/jointset/hip_l/hip_flexion_l/value'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/hip_l/hip_flexion_l/speed','/jointset/hip_l/hip_flexion_l/speed'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/hip_l/hip_adduction_l/value','/jointset/hip_l/hip_adduction_l/value'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/hip_l/hip_adduction_l/speed','/jointset/hip_l/hip_adduction_l/speed'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/hip_l/hip_rotation_l/value','/jointset/hip_l/hip_rotation_l/value'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/hip_l/hip_rotation_l/speed','/jointset/hip_l/hip_rotation_l/speed'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/knee_l/knee_angle_l/value','/jointset/knee_l/knee_angle_l/value'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/knee_l/knee_angle_l/speed','/jointset/knee_l/knee_angle_l/speed'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/ankle_l/ankle_angle_l/value','/jointset/ankle_l/ankle_angle_l/value'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/ankle_l/ankle_angle_l/speed','/jointset/ankle_l/ankle_angle_l/speed'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/hip_r/hip_flexion_r/value','/jointset/hip_r/hip_flexion_r/value'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/hip_r/hip_flexion_r/speed','/jointset/hip_r/hip_flexion_r/speed'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/hip_r/hip_adduction_r/value','/jointset/hip_r/hip_adduction_r/value'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/hip_r/hip_adduction_r/speed','/jointset/hip_r/hip_adduction_r/speed'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/hip_r/hip_rotation_r/value','/jointset/hip_r/hip_rotation_r/value'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/hip_r/hip_rotation_r/speed','/jointset/hip_r/hip_rotation_r/speed'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/knee_r/knee_angle_r/value','/jointset/knee_r/knee_angle_r/value'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/knee_r/knee_angle_r/speed','/jointset/knee_r/knee_angle_r/speed'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/ankle_r/ankle_angle_r/value','/jointset/ankle_r/ankle_angle_r/value'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/ankle_r/ankle_angle_r/speed','/jointset/ankle_r/ankle_angle_r/speed'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/lumbar/lumbar_extension/value'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/lumbar/lumbar_extension/speed'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/lumbar/lumbar_bending/value'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/lumbar/lumbar_bending/speed'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/lumbar/lumbar_rotation/value'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/lumbar/lumbar_rotation/speed'));

However, my output model is jumping backwards. When I comment out the lumbopelvic DOFs, I can move 'normally' but then am I going to have similar issues with not being constrained?

I've also changed the time window to look at a full gait cycle (HS --> HS) but I'm realizing that may not be analogous to a stride (in my head, HS --> contralateral TO)?

Thanks,
Nick