Estimation of passive joint torques

[Axel Koussou]

Hi Thomas,

Thanks for your answer.

As you can see in my previous post (the 7th of this topic), I did use setMinimumActivation().
Nevertheless, by comparing with your code, I did not use realizeVelocity (Could you explain to me why you add this line?). So, I added this line too.
My code is now :

Code: Select all

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
    for i=0:num_muscles-1
        
     MuscMill=Millard2012EquilibriumMuscle.safeDownCast(muscle_set.get(i));
     MuscMill.setMinimumActivation(0);
     MinimumActivation(i+1)=MuscMill.getMinimumActivation(); %Print 0 for all the muscles
     
     MuscMill.setActivation(state,0);
     model.realizeVelocity(state);
     Activation(i+1)=MuscMill.getActivation(state); %Print 0.01 for all the muscles
    end
end

But, I am still not able to set my activation to zero.
What can you suggest me? Do you see any mistake in my code?

Thanks in advance,

Regards

[Thomas Uchida]

Can you confirm that the code I posted works? If so, please check the properties of your muscle- there must be some property that is different between the example that works and the muscle in your model. For example, perhaps you have set ignore_tendon_compliance to true? I could investigate further if you provide a MWE so I can reproduce the issue.

[Axel Koussou]

Hi Thomas,

Thanks for your answer.

I adapt your code for MatLab :

Code: Select all

import org.opensim.modeling.*

model=Model();

muscle=Millard2012EquilibriumMuscle('muscle', 1.0, 0.1, 0.1, 0.0);
vec1=osimVec3FromArray([0 0 0]);
vec2=osimVec3FromArray([0.2 0 0]);
%muscle.addNewPathPoint('p1', model.getGround(), osimVec3FromArray(0,0,0))
muscle.addNewPathPoint('p1', model.getGround(), vec1)
muscle.addNewPathPoint('p2', model.getGround(), vec2)
model.addModelComponent(muscle)

muscle.set_minimum_activation(0.0)
MinActivation=muscle.get_minimum_activation(); %print 0

state = model.initSystem();

muscle.setActivation(state, 0.5)
model.realizeVelocity(state)
Activation=muscle.getActivation(state); %print 0.5

muscle.setActivation(state, 0.0)
model.realizeVelocity(state)
Activation=muscle.getActivation(state); %print 0

And, indeed, it works.

Concerning, the ignore_tendon_compliance properties. In the Rajagopal's model, this properties is sometimes set to true for some muscle and set to false for others. But, it does not seem to influence my result, because, as explained, I got an activation of 0.001 for all my muscles (after setting it to 0). I will check for the other properties.

I do not understand what you mean about providing you a MWE.
Couldn't you try my code with the Rajagopal's model?

Thanks in advance,

Regards

[Thomas Uchida]

In the Rajagopal model, the "minimum_activation" property is set to 0.01 for all muscles. You need to modify this property using the setMinimumActivation() method (https://simtk.org/api_docs/opensim/api_ ... ae3031d499) and this must be done before you call initSystem(). In the script you have uploaded, you are setting the "minimum_activation" property after calling initSystem() which modifies the model but not the underlying computational system. The script below sets all muscle activations to zero:

Code: Select all

import opensim as osim

model = osim.Model('Rajagopal2015.osim')
muscles = model.getMuscles()

# Set the minimum_activation property to 0.0
for i in range(muscles.getSize()):
    mcl = osim.Millard2012EquilibriumMuscle_safeDownCast(muscles.get(i))
    mcl.setMinimumActivation(0.0)

state = model.initSystem()

# Set the activation to 0.0
for i in range(muscles.getSize()):
    mcl = osim.Millard2012EquilibriumMuscle_safeDownCast(muscles.get(i))
    act_original = mcl.getActivation(state)
    mcl.setActivation(state, 0.0)
    act_new = mcl.getActivation(state)
    print '{}: {} --> {}'.format(mcl.getName(), act_original, act_new)

[Axel Koussou]

Hi Thomas,

Thank you for your answer.

It works !
Last question : Just for curiosity, how could I have known that I should use setMinimumActivation() before calling initSystem() ?
I have found the information on the OpenSim API page.

Thanks in advance,

Regards

[Thomas Uchida]

how could I have known that I should use setMinimumActivation() before calling initSystem() ?

The examples and test cases set properties before calling initSystem(). There's also some information on the "API Guide" page in the doxygen documentation: https://simtk.org/api_docs/opensim/api_ ... Guide.html.

[Axel Koussou]

Hi Thomas,

After all, it seem that this is not the last question.

As a reminder, I wanted 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 compute the musculo-tendinous forces after setting 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 the products between passive forces and moments arms.

Below, you can find my MatLab Code :

Code: Select all

% Initialisation
% Pull in the modeling classes straight from the OpenSim distribution
import org.opensim.modeling.*

model=Model;

% Main Code
%Initialize
muscle_set = model.getMuscles();
coordinate_set = model.getCoordinateSet();
PathActuator_set=model.getActuators();

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

%Set the minimum_activation property of all muscles to 0.0
for i=0:num_muscles-1
    MuscMill=Millard2012EquilibriumMuscle.safeDownCast(muscle_set.get(i));
    MuscMill.setMinimumActivation(0.0);
end

state=model.initSystem();

%Compute muscle force at each time point in the motion by updating the
%state and muscle activation (to 0)
muscleforces=zeros(num_timepoints,size(nameIndex,2)); 
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 all the muscles to zero
    for i=0:num_muscles-1
     MuscMill=Millard2012EquilibriumMuscle.safeDownCast(muscle_set.get(i));
     MuscMill.setActivation(state,0);
     
       %model.realizeVelocity(state);
     Activation(t,i+1)=MuscMill.getActivation(state); %Print 0 for all the muscles
    end
    
    model.equilibrateMuscles(state)
    
    %Compute Muscle Forces
    model.realizeDynamics(state);
    
    %DownCast + Passive Force
     for i=1:size(nameIndex,2) %nameIndex : Index of the muscles acting on the ankle DoF only
        MuscMill=Millard2012EquilibriumMuscle.safeDownCast(muscle_set.get(nameIndex(i)-1)); 
	muscleforces(t,i)=MuscMill.getPassiveFiberElasticForceAlongTendon(state);
    end

%DownCast + Moment Arms
    for i=1:size(nameIndex,2) %nameIndex : Index of the muscles acting on the ankle DoF only
        MuscMill=Millard2012EquilibriumMuscle.safeDownCast(PathActuator_set.get(nameIndex(i)-1));  
        momentarm(t,i)=MuscMill.computeMomentArm(state,coordinate_set.get('ankle_angle_r'));
    end
    
end

%Compute Passive Moment
Moment=muscleforces.*momentarm;

for i=1:size(Moment,1)
    TotalMoment(i,1)=sum(Moment(i,:));
end

It works and I am able to obtain the passive moment around the ankle during the mobilization (cf. Attachments).

In a second step, I wanted to validate my results with the GUI. Thus, I decided to plot the passive moments directly through the GUI. To do so, I open my model (the same used in the MatLab Code), I load the motion corresponding to the mobilization.
Then, using the Plot Tool, I sum the moment of all the muscles around the ankle. Moreover, I set the muscle activation to 0 via the Advanced Properties. In the attachments, you can find the details of my plot window.

Nevertheless, when I compare the passive moments obtained through MatLab and through the GUI, the results don't match, as you can in the attachment files.

Do you have any idea of where this difference can come from?

Thanks in advance,

Regards

[Axel Koussou]

Hi all,

Since I have not had any answer, I take the libery of reopening this topic.

Thanks in advance,

Regards