Force-length-velocity curve, API

[Natalia Kosterina]

Hi all,

I am planning to make some changes to Active Fiber Force and then calculate activations for the new force. In Thelenetal2003 there is an equation (10) showing how force and activation are related. Using OpenSim terms, I simplified this equation:
Force_active_fiber=Activation*F_lv (*)

I wrote a plug-in and test it on Gait2354 model (example). We can get the force-length-velocity surface (F_lv) using function Muscle::evaluateForceLengthVelocityCurve(aActivation, aNormalizedLength, aNormalizedVelocity). However, this F_lv does not match with the one calculated from (*).

Could you please explain why the evaluated values do not agree with the equation describing Thelen2003 muscle model?

Thanks for any help,
Natalia

p.s.: I attached the curves for mentioned parameters.

[Matthew Millard]

Hello Natalia,

I just started here at OpenSim, and one of my first tasks was to benchmark the muscle models. The first one I've tested was the Thelen2003 model, and like you I found that the force velocity curves were different than the published paper. We're presently testing out an implementation that is faithful to the Thelen et al 2003 paper, and benchmarking the other muscle models.

As a note, during my bench mark simulations I have found that it is necessary to use the stiff integrator in OpenSim (CPodes), the non-stiff integrators will have convergence issues when activation approaches low values (for example RungeKuttaFeldberg starts to have problems when activation approaches 0.2). I haven't fully tested the range of settings that will work, but the following settings appear to work well with the CPodes integrator provided the activation of the muscle is not allowed to go below 0.05:

maxstep = 0.001;
minstep = 1e-6;
accuracy = 1e-5;
tolerance = 1e-5;
CPodesIntegrator integrator(osimModel.getMultibodySystem(),CPodes::BDF,CPodes::Newton);

If the activation level is allowed to go below 0.05 you may observe a little numerical instability in the plots of d lce/dt. I should note that my benchmark simulation is aggressive, so you may be able to use lower activation values without incurring any numerical instability.

[Natalia Kosterina]

Hello Matthew,

Thanks for the answer. Interesting that you work on solving this discrepancy. When do you think the corrected muscle model will be accessible for the OpenSim users?

I am now confused how to apply my changes to the model. I planned to modify force values at every time step and recalculate Activations using formula from the Thelen et al 2003, then compare the modified activations with the initial ones and with the EMG data. Do you think I can use this way of activation modification if instead of evaluateForceLengthVelocityCurve function I will calculate F_lv as F_act_fib/Activation?

Best,
Natalia

[Matthew Millard]

Hello Natalia
***********************************
Q1: When do you think the corrected muscle >model will be accessible for the OpenSim users?
***********************************

This is a top priority for us, so we're working quickly to correct the muscle models and thoroughly bench mark them.

***********************************
Q2: I am now confused how to apply my changes to the model. I planned to modify force values at every time step and recalculate Activations using formula from the Thelen et al 2003, then compare the modified activations with the initial ones and with the EMG data. Do you think I can use this way of activation modification if instead of evaluateForceLengthVelocityCurve function I will calculate F_lv as F_act_fib/Activation?
***********************************

Just so I'm clear on what you're doing, are these the steps you are trying to follow for you investigation:

1. Simulate some movement and save the muscle forces
2. Perturb the simulation-generated muscle forces
2. Back solve for activation given the perturbed muscle forces using the Thelen model.

Is this correct? I'll wait until I hear from you to answer the remainder of your question, just to avoid undue confusion if I've misunderstood.

[Natalia Kosterina]

Hello,

yes, this is correct. Let us use some symbols for these steps:

1) Simulation -> Force, Activation, F_lv=Force / Activation
2) Force_mod=Force+modification
3) Activation_mod=Force_mod / F_lv.

Then I will compare Activation and Activation_mod with the EMG data. After all this can I refer to the Thelen et al 2003 paper?