Computing muscle fiber lengths using the Python API.

[Pranav Mamidanna]

Hi all,

I am using OpenSim 4.0 API (Python) in order to calculate muscle length changes given a state trajectory (sequence of coordinate values through a motion file). As a first test, I used the Arm26 model and the motion file within to verify if my results match. The model can be found at '<source>/Models/Arm26/arm26.osim' and the motion file I use can be found at '<source>/Models/Arm26/OutputReference/InverseKinematics/arm26_InverseKinematics.mot'.

1. In order to compute the muscle lengths at each time point in the motion file, I use the Analyze->MuscleAnalysis tool through the OpenSim 3.3 GUI, and solve for equilibrium states which generates a file 'arm26_MuscleAnalysis_FiberLength.sto'. Now, using the following code, I compute the same through the Python API:

Code: Select all

# Load the corresponding opensim model and motion file; initialize.
arm26 = osim.Model("C:/OpenSim 3.3/Models/Arm26/arm26.osim")
motion_filepath = "C:/OpenSim 3.3/Models/Arm26/OutputReference/InverseKinematics/arm26_InverseKinematics.mot"
init_state = arm26.initSystem()

# Declare generalized coordinate values for the trajectory
gen_coordinates = osim.STOFileAdapter.read(motion_filepath) # returns a TimeSeriesTable object
num_timepoints = gen_coordinates.getNumRows()
num_coordinates = gen_coordinates.getNumColumns()
muscle_set = arm26.getMuscles()
num_muscles = muscle_set.getSize()

# Compute muscle lengths at each time point in the motion by updating the state
musclelengths = np.zeros((num_muscles,num_timepoints))
for timepoint in range(num_timepoints):
    for i in range(num_coordinates): # I make sure that the order of specifying generalized coordinates is preserved
        arm26.getCoordinateSet().get(i).setValue(init_state, 
                np.radians(gen_coordinates.getDependentColumnAtIndex(i)[timepoint]), False) 
    arm26.assemble(init_state)
    arm26.equilibrateMuscles(init_state)
    
    for muscle_num in range(num_muscles):
        muscle = muscle_set.get(muscle_num)
        musclelengths[muscle_num, timepoint] = muscle.getFiberLength(init_state)

Here, I step through each time point, update the coordinate values, assemble the kinematics of the state and use the equilibriateMuscles() routine to find the muscle states. Is this a sensible way to proceed? This gives me almost identical results to the fiber length values obtained through the GUI. However, when I modify the motion file to discard 'r_elbow_flexion' and only keep 'r_shoulder_elev', I notice a clear difference in the computed muscle lengths. But I do not know why this is so!

MuscleAnalysis_API.jpg (44.23 KiB) Viewed 1340 times

Figure 1: Muscle analysis through the GUI vs through the API. (Left) When the specified motion file has values for both coordinates: 'r_shoulder_elv' and 'r_elbow_flexion'. Notice that the results obtained through MuscleAnalysis from GUI and through my code overlap almost perfectly (there are two lines, on top of each other)! (Right) The specified motion file only has values for 'r_shoulder_elv'. Notice now that the fiber-length curves are not identical, but only so for one of the muscles. <Sorry for not creating dashed lines for easier comparison (like in Figure 2)> the top blue line with two bumps is obtained through my code above, while the bottom line is obtained through the Analyze tool.

When I repeat this with the Stanford VA UpperLimb model, where I specify values for only 4 coordinates, things are much worse and very different from the results I obtain using the GUI.

2. Apart from this, I also noticed that using the plot tool's motion curves (where one can specify a motion file associated with the model to plot muscle length variations) gives quantitatively different values for the muscle lengths, although the fiber-length vs time curves look qualitatively similar. Why could this be so?

MuscleAnalysis_Plot.jpg (66.14 KiB) Viewed 1340 times

Any help is much appreciated.

Thanks a lot.

[Dimitar Stanev]

Hi,

This is from the API (maybe you can try to assign high values to the weight) :

Code: Select all

//Find the kinematic state of the model that satisfies constraints and coordinate goals If assemble is being called due to a coordinate set value, provide the option to weight that coordinate value more heavily if specified. 
void OpenSim::Model::assemble(SimTK::State& state, const Coordinate *  coord = NULL, double  weight = 10 ); 

You can also try something else. Update the coordinates of the model arm26.updCoordinateSet().get(i).setValue(init_state, value). Then call model.realizePosition(init_state) and recompute the muscle lengths without calling equilibrateMuscles.

Best

[Pranav Mamidanna]

Hi Dimitar,

Thanks for the quick response. I have tried all the ways you mentioned above, but it seems like I could not make them work! There is no change in the results by altering the weights in the assemble method, while using realizePosition does not update the muscle lengths at all - in fact by doing so I get flat fiber-length vs time curves. Any thoughts on this?

[Dimitar Stanev]

Could you post the code of my suggestion?

[Pranav Mamidanna]

1. As per your first suggestion, I have changed the weights for the coordinates in the assemble method between 20 and 100 and found no difference in the fiber-length vs time curves.

2. Using realizePosition(init_state) instead of equilibrateMuscles(init_state) somehow leads to no variation in the fiber-lengths i.e., flat profiles (the state doesn't get updated? or muscle lengths are not updated by using realizePosition()?). I surely must be missing something here. See the code below:

Code: Select all

# Compute muscle lengths at each time point in the motion by updating the state
musclelengths = np.zeros((num_muscles,num_timepoints))
for timepoint in range(num_timepoints):
    for i in range(num_coordinates):
        arm26_2.updCoordinateSet().get(i).setValue(init_state, 
                np.radians(gen_coordinates.getDependentColumnAtIndex(i)[timepoint]), False)
    arm26_2.realizePosition(init_state)
    
    for muscle_num in range(num_muscles):
        muscle = muscle_set.get(muscle_num)
        musclelengths[muscle_num, timepoint] = muscle.getFiberLength(init_state)

The motion file I use for the arm26 model other than the one I mention in my first post is:

arm26_GC_ShoulderOnly.mot

(2.98 KiB) Downloaded 14 times

3. Just to show you the magnitude of difference I obtain using my code as against the MuscleAnalysis Tool from the GUI, I attach the following figure: here, I use the Stanford VA UpperLimb model with 50 muscles to generate fiber-lengths vs time.

UpperLimbModel_Difference.jpg (87.82 KiB) Viewed 1263 times

The code used to generate this figure is:

Code: Select all

# Load the corresponding opensim model and motion file; initialize.
home_dir = '<home>'
model_dir = '/UpperExtremityModel-latest/Original Files - UpperExtremityModel/'
ulmodel = osim.Model(home_dir + model_dir + 'Stanford VA upper limb model_0.osim')
init_state = ulmodel.initSystem()

# Declare generalized coordinate values for the trajectory
gen_coordinates = osim.STOFileAdapter.read(home_dir + model_dir + 'test.mot') # returns a TimeSeriesTable object
num_timepoints = gen_coordinates.getNumRows()
num_coordinates = gen_coordinates.getNumColumns()
muscle_set = ulmodel.getMuscles()
num_muscles = muscle_set.getSize()

# Compute muscle lengths at each time point in the motion by updating the state
coordinates = ['elv_angle', 'shoulder_elv', 'shoulder_rot', 'elbow_flexion']
musclelengths = np.zeros((num_muscles,num_timepoints))
for timepoint in range(num_timepoints):
    for i in range(num_coordinates): # I make sure that the order of specifying generalized coordinates is preserved
        ulmodel.getCoordinateSet().get(coordinates[i]).setValue(init_state, 
                np.radians(gen_coordinates.getDependentColumnAtIndex(i)[timepoint]), False) 
    ulmodel.assemble(init_state)
    ulmodel.equilibrateMuscles(init_state)
    
    for muscle_num in range(num_muscles):
        muscle = muscle_set.get(muscle_num)
        musclelengths[muscle_num, timepoint] = muscle.getFiberLength(init_state)

The motion file I use in this case is:

test.mot

(9.23 KiB) Downloaded 13 times

Any idea where things go awry?

[Dimitar Stanev]

Hi,

You need to use the getLength() function. This method calculates the actual length of the musculotendon. The getFiberLength methods return the state variable of the muscle length that is used for the muscle dynamics. That is why you had to equilibrate the muscles in order to get valid results. This agrees with the GUI plot tool (muscle-tendon-length). Arm26:

Code: Select all

import opensim
import numpy as np

model = opensim.Model('arm26.osim')
state = model.initSystem()

coordinate = 'r_elbow_flex'
elbow = np.arange(0, 1.57, 0.1)
muscle = 'BIClong'
BIClong = []

for q in elbow:
    model.updCoordinateSet().get(coordinate).setValue(state, q)
    model.realizePosition(state)
    BIClong.append(model.getMuscles().get(muscle).getLength(state))

Best

[In Bae Chung]

Hi,

Thanks for the codes you guys provided

You can also try something else. Update the coordinates of the model arm26.updCoordinateSet().get(i).setValue(init_state, value). Then call model.realizePosition(init_state) and recompute the muscle lengths without calling equilibrateMuscles.

As you said, I tried using model.realizePosition(state) instead of equilibrateMuscles.

You need to use the getLength() function. This method calculates the actual length of the musculotendon. The getFiberLength methods return the state variable of the muscle length that is used for the muscle dynamics.

Also, I used the getLength() function instead of getFiberLength and was able to get the value for fiber length.
However, when I activated the muscles using the fiber lengths I got and tried to get the fiber force, the result is 0.

Code: Select all

import org.opensim.modeling.*
model = Model('C:\Users\In Bae Chung\Desktop\Arm swing simulation data (1)\arm_swing\scaledmodel.osim');
state = model.initSystem(); 
elbowangle=model.getCoordinateSet().get(1);
elbowangle.setValue(state,45);
model.realizePosition(state)

muscle=model.getMuscles().get(0);
musclegeopath=muscle.getGeometryPath();

FiberLength=muscle.getLength(state);

muscle.setActivation(state, 0.1);
afl = ActivationFiberLengthMuscle.safeDownCast(muscle);
afl.setFiberLength(state, FiberLength);
model.realizePosition(state)
FiberForce=muscle.getFiberForce(state);

Is it because getFiberForce method returns the muscle force that is used for the muscle dynamics just like getFiberLength?
In that case, is there another function I could use to get the correct fiber force value?

Thanks.

[Thomas Uchida]

I used the getLength() function instead of getFiberLength and was able to get the value for fiber length.

getLength() is a method on GeometryPath (https://simtk.org/api_docs/opensim/api_ ... 55bc15dc8c); you're getting the entire length of the actuator (muscle plus tendon), not just the length of the muscle fibers.

However, when I activated the muscles using the fiber lengths I got and tried to get the fiber force, the result is 0.

Fiber force requires dynamics-level calculations because it is a force; "model.realizePosition(state)" would need to be "model.realizeDynamics(state)" (https://simtk.org/api_docs/opensim/api_ ... 96e8844fe7).

[In Bae Chung]

Thank you Thomas.
I was able to run the code without error as follows

Code: Select all

for i = 0 : nt - 1
    model.updCoordinateSet().get('elbow_flexion').setValue(s, deg2rad(motion.getDependentColumnAtIndex(38).get(i) ));
    coord = model.getCoordinateSet().get('elbow_flexion');

    model.realizeDynamics(s);
    
    % Get the computed force from each muscle
    for u = 0 : nm - 1 
         % Get a muscle
        m = ms.get( char(muscleNames.get(u)));
        % Compute muscle force
        mf(i+1,u+1) = m.getFiberForce(s);
    end
end

However, the results for the muscle force are all 'NaN.'
Meanwhile, when I used 'getTendonForce' instead of 'getFiberForce' it does give the results.

1. Why does the API give 'NaN' as results for 'getFiberForce' ? Is there a possible way to get the correct results by changing the code?

I found that fiber force is the force applied to the tendon and tendon force is the force applied to the bones.
I want to compute the joint moments by multiplying the moment arm of each muscles with the muscle forces.

2. Which one out of fiberforce and tendonforce should be multiplied with the moment arm to get the correct joint moment?
I'm not sure if this is the right way to compute joint moment. If so, is there any document that I can look into to find the right way to compute joint moment?

I would really appreciate your help.
Thanks a lot.

[Thomas Uchida]

is there any document that I can look into to find the right way to compute joint moment?

Please see Sherman, M.A., Seth, A., Delp, S.L. What is a moment arm? Calculating muscle effectiveness in Biomechanical models using generalized coordinates, Proceedings of the ASME 2013 International Design Engineering Technical Conferences. https://nmbl.stanford.edu/wp-content/up ... -13633.pdf