DeGrooteFregly2016 muscle activation and deactivation speeds

[tvdbogert]

The activation dynamics is implemented exactly as in the 2016 paper, but this results in activation and deactivation happening at almost the same speed. The function f is intended as a "soft switch" between activation (a-e < 0) and deactivation (a-e > 0). But it's too soft, the value of f stays very close to 0.5 over the whole range of inputs (a-e ranging from -1 to 1).

I suspect that the tanhSteepness parameter should have been 10 instead of 0.1.

See:

opensim-core/OpenSim/Actuators/DeGrooteFregly2016Muscle.cpp

Line 208 in b4ba7ca

static const double tanhSteepness = 0.1;

I ran some simulations to see the response to a 1 Hz square wave excitation, and this confirms that activation and deactivation speeds are very similar.

I contacted Friedl for input, a few weeks ago, but received no reply yet.

[nickbianco]

Hi @tvdbogert, thanks for bringing this to our attention. Perhaps we could add a property so users could control the amount of smoothing applied by tanhSteepness.

[tvdbogert]

That would be a good solution.

For backward compatibility, you probably want to leave the default value as it is now, so previous results are reproducible, but allow users to change the value.

To make users aware of this issue, maybe put something in the documentation, saying it's recommended to change the value? Or generate a warning when the model is used with the default parameter value?

There is more to this issue though! When you have a substantial difference between activation and deactivation speeds (as intended), a cost function with squared excitation is no longer good, because it favors a solution with infinitely fast on-off switching of the excitation. Indeed that is what I see emerge during mesh refinement. That's not great for accuracy and solution speed of collocation methods. Instead, activation state should be used in the cost function. I'm working on a report about this and will share soon.

Does MOCO use excitation or activation in the effort cost?

[nickbianco]

For backward compatibility, you probably want to leave the default value as it is now, so previous results are reproducible, but allow users to change the value.

To make users aware of this issue, maybe put something in the documentation, saying it's recommended to change the value? Or generate a warning when the model is used with the default parameter value?

Agreed on all this.

Does MOCO use excitation or activation in the effort cost?

You can use MocoControlGoal for an excitation effort cost and MocoSumSquaredStatesGoal if you want an activation effort cost (or both if you want both). I don't think its clear to some users that MocoSumSquaredStatesGoal allows you do have an activation effort cost, I might consider a MocoActivationGoal or something similar.

Interested to read your report on this!

[mrrezaie]

I ran some simulations to see the response to a 1 Hz square wave excitation, and this confirms that activation and deactivation speeds are very similar.

Hi @tvdbogert, I'm very interested in this topic, and wondering if you could share the code for this simulation, so I will learn more. Thanks!!!

I had tried to understand where that equation has come from. Based on the cited references in the DeGroote2016 paper, that is the smoothed approximation of the first-order activation dynamics:

And based on this guide in OpenSim documentations, the approximated equation should be:

But they lead to different outputs:

import numpy as np

# parameters
t_a = 0.015 # (s)
t_d = 0.060 # (s) 
b = 0.1
e = 1
a = 0.5

# De Groote et al. 2016 Eq(1) and Eq(2)
f = 0.5*np.tanh(b*(e-a))
const = (f+0.5)/(t_a*(0.5+1.5*a)) + ((0.5+1.5*a)/t_d)*(-f+0.5)
print((e-a)*const) # = 18.94757846319944

# OpenSim guide
f = 0.5 + 0.5*np.tanh(b*(e-a))
const = (t_d/(0.5+1.5*a))*(1.0-f) + (t_a*(0.5+1.5*a))*f
print((e-a)/const) # = 15.316582064925484

Thank you in advance.

[tvdbogert]

My code and report are here: https://github.com/csu-hmc/activation-dynamics

Comments or questions are welcome!

I coded De Groote's activation dynamics directly from equation 1 and 2 in the paper.

Indeed the OpenSim code does not exactly the same as in the paper. Good observation! OpenSim computes the weighted average of the activation and deactivation time constants. The paper computes the weighted average of the inverse of the time constants (the activation and deactivation rates). The result is similar but not the same.

@nickbianco please take a look, this may also need to be corrected if you want to be exactly consistent with De Groote's paper.

The factor 0.5+1.5a was somewhat mysterious to me, but is explained by Winters (1995). This factor will slow down the activation, by as much as a factor 2, when the activation gets close to 1. And the deactivation will also slow down by a factor 2 when the activation gets close to 0. It seems more complex than we need, or can justify.

At the end of my report, I propose an improvement to the McLean2003 activation model. It simply uses the sigmoid function f to interpolate between the activation and deactivation rates (like De Groote, but without the 0.5+1.5a factor). It's also robust for cases where activation and excitation go above 1, as you sometimes might do to compensate for insufficient muscle strength in the model.

[nickbianco]

@mrrezaie that guide is for the Bhargava metabolics model, not the smoothed activation dynamics model. The activation dynamics model is implemented according to DeGroote et al. 2016 here:

opensim-core/OpenSim/Actuators/DeGrooteFregly2016Muscle.cpp

Line 209 in c46cc98

// f = 0.5 tanh(b(e - a))

[mrrezaie]

@tvdbogert So Nice!!! Thanks a lot for sharing your codes and the report.

@nickbianco Thanks for your comment. Do you mean the approximation method using hyperbolic tangent function (tanh) differs between the models? I thought they should follow a same rule. Sorry!!!

(I compared the original DeGroote2016 activation dynamics (with b=0.1 parameter) to the McLean2003Improved for solving muscle redundancy on real walking data (75 frames), and the later model converged significantly faster (7.3 vs. 4.9). Also, the DeGroote2016 with b=10 took much longer to converge (68 s) and there were so many bumps/noises. This test might be data-specific).

Thanks again.

[tvdbogert]

@mrrezaie That is good to know. With b=10 the model has a stronger nonlinearity in the switching function f, which can make it harder for the optimal control solver to find the solution.

[nickbianco]

Do you mean the approximation method using hyperbolic tangent function (tanh) differs between the models? I thought they should follow a same rule. Sorry!!!

No worries. The both use hyperbolic tangent functions, just in slightly different ways to achieve the desired smoothing.