Estimation of passive joint torques

[Axel Koussou]

Dear OpenSim experts,

For a project, I would like to estimate the passive joint torques during a passive (without muscle activity) sollicitation.
Experimentally, we mobilize the ankle of a subject with an handheld dynamometer, and we measure the kinematics thanks to VICON.

Thus, I am able to animate my model (Rajagopal's model) - I mean I impose the kinematic on it - and I would like to know the joint torque during this passive movement.

To do so, I would like to compute the musculo-tendinous forces after setting (I don't know how) the activations of all the muscles to 0. Thus the calculated musculo-tendinous forces should be the passive forces.

Then, using the moment arms, I should be able to obtain my passive joint torque by summing, for each degree of liberty, the products between passive forces and moments arms.

Do you think it is possible? And, how can I achieve such a process?

I've also got the forces during the sollicitation through the dynamometer but I don't think it would be necessary.

Thanks in advance,

Regards

[Carmichael Ong]

The doxygen page for the Class you're working with is a helpful place to start when looking for how to perform some actions and get some quantities of interest. I think that model uses the Millard2012EquilibriumMuscle: https://simtk.org/api_docs/opensim/api_ ... uscle.html

In particular, make sure to check under the functions inherited from parent classes, such as setActivation() from Muscle (https://simtk.org/api_docs/opensim/api_ ... ea7c9c712a) and computeMomentArm() from PathActuator (https://simtk.org/api_docs/opensim/api_ ... 6c06767289)

[Axel Koussou]

Hi Carmichael,

Thank you for your reply.

I'm already aware of this page, which indeed can help me. Nevertheless, there is a grey area on how I should operate to get the passive moments.
The thing is we compute experimentally by inverse dynamics this passive moments during an ankle mobilization and we want to compare what we obtain to what the model predicts. So, I want ot be sure that what I compute on OpenSim correspond at what we obtain experimentally.

If I understand well what you said, I should be able to obtain the passive moments by following, for every frames of my motion, such a process :

- Update the state of my model
- Set the activation of my muscle at 0
- Compute the moment arms from PathActuator
- Compute the forces
- Multiply the forces and the moment arms

Am I right until here?

But, one of my problem concern the function I should use to obtain the forces. It seems that I got the choice between getPassiveFiberElasticForce or getPassiveFiberElasticForceAlongTendon. (There is also getPassiveFiberDampingForce and getPassiveFiberDampingForce AlongTendon)
Intuitively, I would use getPassiveFiberElasticForceAlongTendon.
"AlongTendon" so I could then multiply by the Moment Arms. But why the "PassiveFiberElasticForce" only ?

Maybe, I'm misunderstanding but I would have thought that I should use the forces all the MTU, not only the Passive Fiber Elastic. So, I thought that I should use : (Tendon Force + Muscle Force) * Moment Arms.

So, here is my first question : Which forces do you advice me to use in order to obtain the passive moment at the ankle during a passive sollicitation ?

Another point is that when I run an MuscleAnalysis in the GUI, the PassiveFiberForceAlongTendon that I obtain are 0 for all the plantarflexor muscles. Why is this? Will I get the same result if I run the process explained above?

Finally, I'm able to plot a curve that is near of what I want, that you can see in the Attachments. This curve is obtained, in the Plot Tool, by setting the activation of all the muscles to zero, and computed the ankle moment by summing the moments generated by all the muscles that could contribute to the joint moment. This is what Arnold et al. 2010, or Rajagopal et al. 2016, do when they compare the passive moments determined by the model to experimental ones.
So in theory, I should be able to use this curve to compare the moments that we obtain experimentally.
But I do not understand how this curve is obtained, I mean what are the equations behind it. If, I run the process explained above (with getPassiveFiberElasticForceAlongTendon and so on..) for the same angle range, and with all the other joints set to the same states, would I obtain the same curve?

Sorry for all these question, I would be really grateful if you can answer to all. I am starting a PhD so there is still a lot of things that I do not understand well.
I hope I've been clear enough. If not, do not hesitate to ask me for more informations.

Thank you in advance,

Regards

[Axel Koussou]

Hi,
I try to implement the code suggest by Carmichael, using the function from the Muscle parent class, to compute the muscle forces at each time point in the motion by updating the state and setting the muscle activation to 0.
Unfortunately my code doesn't work. Below you can find a snippet of it :

Code: Select all

import org.opensim.modeling.*

model=Model;
state=model.initSystem();

%Initialize
muscle_set = model.getMuscles();
Millard=Millard2012EquilibriumMuscle();

% Data : Values of coordinates in radian. (pelvis_tilt, pelvis_list...)
num_timepoints = size(Data,1);
num_coordinates = size(Data,2);

num_muscles = muscle_set.getSize();

muscleforces=zeros(num_timepoints,num_muscles);  

%Compute muscle force at each time point in the motion by updating the
%state and muscle activation (to 0)

for t=1:num_timepoints
    
    %Update state
    for i=0:num_coordinates-1
        model.updCoordinateSet().get(i).setValue(state,Data(t,i+1),false);
    end
    clear i
    model.assemble(state)
    
    % Set activation of the muscles to zero
    % The first line is supposed to set the activation of the muscle to
    % zero. 
    % But if I run the second line, I do not get 0 but 0.01. 
    for i=0:num_muscles-1
     muscle_set().get(i).setActivation(state,0);
     %muscle_set().get(i).getActivation(state);
    end
    
    model.equilibrateMuscles(state)
    
    %Compute Muscle Forces
    model.realizeDynamics(state);
    
    %This line doesn't run. 
    %Error : Check for missing argument or incorrect argument data type in call to function 'getPassiveFiberElasticForceAlongTendon'.
    for i=0:num_muscles-1
        muscleforces(t,i+1)=muscle_set.get(i).getPassiveFiberElasticForceAlongTendon(state);
    end
    
    %Compute Moment Arms
    
    %Work in progress..
    
end

First, I don't know why but I can't set the muscle activation to 0. As explained in the snippet, when I try to set the activation to 0, it appears that the activation is in fact equal to 0.01. Am I missing something ? How can I proceed?

Secondly, as quoted in the snippet, I got a message error when I try to compute my PassiveFiberElasticForceAlongTendon. Do you know where my mistake is?

Finally, a more general question. I know that the Rajagopal's model uses Milliard2012EquilibriumMuscle. So, I would like to know if it is ok if I use function from Muscle class and not Millard2012EquilibriumMuscle()? Have I to precise elsewhere that the model use Millard2012EquilibriumMuscle()?

Thanks in advance,

Regards

[Thomas Uchida]

getPassiveFiberElasticForceAlongTendon() is a method in the Millard2012EquilibriumMuscle class (https://simtk.org/api_docs/opensim/api_ ... 56898bdd2c), but the reference returned by "muscle_set().get(i)" is of type Muscle. You need to downcast to access getPassiveFiberElasticForceAlongTendon(). Please search for "downcast" on the "Common Scripting Commands" page in the Confluence documentation: https://simtk-confluence.stanford.edu/d ... g+Commands.

[Axel Koussou]

Hi Thomas,

Thank you for your answer. It works well now, I can obtain the passive fiber elastic force along the tendon.

Nevertheless, I am still not able to set my muscle activation to 0.

Do you have any idea of what I'm doing wrong ?

Thanks in advance,

Regards

[Axel Koussou]

Here is the code that I'm using:

Code: Select all

    for i=0:num_muscles-1
        
     MuscMill=Millard2012EquilibriumMuscle.safeDownCast(muscle_set.get(i));
     MuscMill.setMinimumActivation(0);
     MuscMill.setActivation(state,0);
     MuscMill.getActivation(state);
     %muscle_set().get(i).setActivation(state,0);
     %muscle_set().get(i).getActivation(state);
    end

Thanks

[Thomas Uchida]

Please read the documentation for the Millard2012EquilibriumMuscle class: https://simtk.org/api_docs/opensim/api_ ... ml#details. "When an elastic tendon without fiber damping is selected, the minimum active-force-length value is set to 0.1, the minimum permissible activation is set to 0.01, and the maximum permissible pennation angle is set to acos(0.1) or 84.3 degrees. This is done as a convenience for the user to prevent the model from taking an unreasonable amount of time to simulate." There is more detail on that page and in the accompanying paper: Millard et al. (2013) "Flexing computational muscle: modeling and simulation of musculotendon dynamics" http://nmbl.stanford.edu/publications/p ... rd2013.pdf.

[Axel Koussou]

Hi Thomas,

Thanks for your answer.
Of course, I read the documentation for the Millard2012EquilibriumMuscle class. But may be, I did not understand it well.
The quote you made specify that this condition of a non-zero activation is for "an elastic tendon without fiber damping".
Nevertheless in my model (the Rajagopal's one), the fiber damping is included. Indeed, the value setting this parameter is greater than 0.001. The fiber damping parameter of all the muscles is set to 0.01 and the documentation for the Millard2012EquilibriumMuscle class define :

"fiber_damping: set to a value greater than 0.001 to include fiber damping in the model. The addition of damping reduces simulation time while allowing the muscle model to be more physiological (it can have an activation of zero, its active-force-length curve can go to zero, and its force-velocity curve can be asymptotic)."

That's why I do not understand why I am not able to set my activation to zero. So, do you have any ideas to solve my problem?

Thanks in advance,

Regards

[Thomas Uchida]

I don't see a MWE (https://en.wikipedia.org/wiki/Minimal_working_example) to reproduce the issue so I have written the Python script below, which works as expected. Perhaps you have not called muscle.set_minimum_activation(0.0)?

Code: Select all

import opensim as osim

model = osim.Model()

muscle = osim.Millard2012EquilibriumMuscle('muscle', 1.0, 0.1, 0.1, 0.0)
muscle.addNewPathPoint('p1', model.getGround(), osim.Vec3(0,0,0))
muscle.addNewPathPoint('p2', model.getGround(), osim.Vec3(0.2,0,0))
model.addModelComponent(muscle)

muscle.set_minimum_activation(0.0)
print muscle.get_minimum_activation() #prints 0.0

state = model.initSystem()

muscle.setActivation(state, 0.5)
model.realizeVelocity(state)
print muscle.getActivation(state) #prints 0.5

muscle.setActivation(state, 0.0)
model.realizeVelocity(state)
print muscle.getActivation(state) #prints 0.0