Custom Muscle Model Part One

The steps covered in part one are:

The header file (LiuThelen2003Muscle.h)

We start defining our new muscle model by creating a header file. The model will be based on Thelen2003Muscle, and will include fatigue effects as defined in the following paper:
Liu, Jing Z., Brown, Robert, Yue, Guang H., "A Dynamical Model of Muscle Activation, Fatigue, and Recovery," Biophysical Journal, Vol. 82, Issue 5, pp. 2344-2359 (2002).

So we will call our model LiuThelen2003Muscle. At the top of LiuThelen2003Muscle.h, we include the header file for the base class and derive our new class.

#include <OpenSim/Actuators/Thelen2003Muscle.h> 
namespace OpenSim { 
 
class LiuThelen2003Muscle : public Thelen2003Muscle 
{ 

Defining properties

Our new muscle model will have all of the properties of the Thelen2003Muscle model, plus two additional ones to define the rates at which active muscle fibers fatigue and the rate at which fatigued fibers recover.

class LiuThelen2003Muscle : public Thelen2003Muscle 
{
  protected:
  // the rate at which active muscle fibers become fatigued
  PropertyDbl _fatigueFactorProp;
  double &_fatigueFactor; 
 
  // the rate at which fatigued fibers recover (become active)
  PropertyDbl _recoveryFactorProp;
  double &_recoveryFactor; 

Defining states

Thelen2003Muscle has two states: activation and fiber length. We will add two more states: one for the number of active motor units (range 0.0 to 1.0), and one for the number of fatigued motor units (range 0.0 to 1.0).

In LiuThelen2003Muscle.cpp, we will set the indices for the new states. In the header file we just define the names.

  protected:
  static const int STATE_ACTIVE_MOTOR_UNITS;
  static const int STATE_FATIGUED_MOTOR_UNITS;

Required functions

The muscle base class (Muscle) has a number of pure virtual functions. These functions must be defined in all classes that derive from Muscle. To see the complete set of these functions, look for function declarations in Muscle.h that end in " = 0;". Thelen2003Muscle defines all of these functions, such as getFiberLength() and computeIsometricForce(). Because our new class derives from Thelen2003Muscle and not directly from Muscle, we only have to define the base class functions whose behavior we want to change. We will also be adding additional functions to compute the new states for fatigue. The base class functions that we will focus on in this example are:

  public:
  virtual void computeEquilibrium(SimTK::State& s ) const;
  virtual double computeActuation(const SimTK::State& s) const;
  virtual double computeIsometricForce(SimTK::State& s, double act) const;
  virtual void equilibrate(SimTK::State& state) const;

  private:
  void setNull();
  void setupProperties(); 

The source file (LiuThelen2003Muscle.cpp)

We will put the function definitions for our class into LiuThelen2003Muscle.cpp. This source file will include all of the required functions that we need to override in Thelen2003Muscle, as well as new functions that describe the behavior of the fatigue states.

Enumerating the states

At the top of the source file, we define the indices of the muscle states that we are adding to the base class. Thelen2003Muscle defines 0 as the activation state and 1 as the fiber length state. It is important to use 2 and 3 for the states we are adding because the indices for all of the muscle states must start at 0 and be contiguous. Once we define the integer indices STATE_ACTIVE_MOTOR_UNITS and STATE_FATIGUED_MOTOR_UNITS, those names will be used throughout the code to access the state variables, rather than the numbers 2 and 3.

// States 0 and 1 are defined in the base class, Thelen2003Muscle.
const int LiuThelen2003Muscle::STATE_ACTIVE_MOTOR_UNITS = 2;
const int LiuThelen2003Muscle::STATE_FATIGUED_MOTOR_UNITS = 3; 

setNull()

The function setNull() is used to set some of the basic elements of the muscle class, such as its type (name) and its number of states. It is called by OpenSim when an object of this class is constructed.

void LiuThelen2003Muscle::setNull()
{
  setType("LiuThelen2003Muscle"); 
  setNumStateVariables(4); 
  _stateVariableSuffixes[STATE_ACTIVE_MOTOR_UNITS]="active_motor_units";
  _stateVariableSuffixes[STATE_FATIGUED_MOTOR_UNITS]="fatigued_motor_units";
} 

setupProperties()

The function setupProperties() is used to define the properties of the muscle class and add them to the set of all OpenSim properties. This enables you to define them in an OpenSim model file. This function is called by OpenSim when an object of this class is constructed. As we did in the header file, we need to define a property for the rate at which active muscle fibers fatigue (which we will call fatigue_factor) and one for the rate at which fatigued fibers recover (recovery_factor). These factors have a default value of 0.0, are assumed to be in the range 0.0 to 1.0, and are normalized (so they are usually the same for all muscles).

void LiuThelen2003Muscle::setupProperties()
{
  _fatigueFactorProp.setName("fatigue_factor");
  _fatigueFactorProp.setValue(0.0);
  _fatigueFactorProp.setComment("percentage of active motor units that fatigue in unit time");
  _propertySet.append(&_fatigueFactorProp, "Parameters"); 
  _recoveryFactorProp.setName("recovery_factor");
  _recoveryFactorProp.setValue(0.0);
  _recoveryFactorProp.setComment("percentage of fatigued motor units that recover in unit time");
  _propertySet.append(&_recoveryFactorProp, "Parameters");
} 

equilibrate() and computeEquilibrium()

The function equilibrate() computes values of the muscle states assuming the muscle is in an equilibrium state. It should be called for each muscle before you begin a dynamic simulation. In our muscle class, this function initializes the states to reasonable values, realizes the Simbody model to the velocity stage, and then calls computeEquilibrium(). computeEquilibrium() gets the current activation of the muscle and calls computeIsometricForce(), which is described in Section 4.4.3.6.

void LiuThelen2003Muscle::equilibrate(SimTK::State& state) const
{
  // Reasonable initial activation value
  setActivation(state, 0.01);
  setFiberLength(state, getOptimalFiberLength());
  setActiveMotorUnits(state, 0.0);
  setFatiguedMotorUnits(state, 0.0);
  _model->getMultibodySystem().realize(state, SimTK::Stage::Velocity); 

  // Compute isometric force to get starting value of _fiberLength.
  computeEquilibrium(state);
} 

 void LiuThelen2003Muscle::computeEquilibrium(SimTK::State& s) const
{
  double force = computeIsometricForce(s, getActivation(s));
}  

computeActuation()

computeActuation() computes the values of the muscle states and their derivatives. It is called by the integrator or other code whenever the musculoskeletal model's state has changed and the muscle must be updated accordingly. The code for this function in the LiuThelen2003Muscle class is very similar to the code in Thelen2003Muscle::computeActuation(), but includes calculations of the fatigue states. This function performs four basic steps:

• Calculate normalized values of the muscle states from the State object passed in. The State object contains the current values and derivatives of the model's coordinates, muscle states, and any other states in the musculoskeletal model. During a forward dynamics simulation, the integrator operates on the state; it uses the current state, their derivatives, and the equations of motion to determine the state at the next time step. So that the same differential equations can be used for every muscle of the same type, the equations are normalized to certain muscle parameters. For example, tendon length is usually normalized by dividing it by the resting length of the tendon for that muscle, thus producing tendon strain. Tendon force is usually normalized by dividing it by the maximum isometric force of the muscle fibers. Once this is done, the same force vs. tendon strain equation can be used for all muscles.

• Compute normalized derivatives of the muscle states, using the current values of the muscle states and the muscle state equations. For the fatigue and recovery states in our LiuThelen2003Muscle, the equations look like this:

normStateDeriv[STATE_ACTIVE_MOTOR_UNITS] =
  normStateDeriv[STATE_ACTIVATION] –
  getFatigueFactor() * getActiveMotorUnits(s) +
  getRecoveryFactor() * getFatiguedMotorUnits(s);
  normStateDeriv[STATE_FATIGUED_MOTOR_UNITS] =
  getFatigueFactor() * getActiveMotorUnits(s) –
  getRecoveryFactor() * getFatiguedMotorUnits(s); 

The first equation means that the rate of change in the percentage of active motor units is equal to the rate of change in the activation level minus the fatigue factor times the number of active motor units, plus the recovery factor times the number of fatigued motor units. The second equation means that the rate of change in the percentage of fatigued motor units is equal to the fatigue factor times the number of active motor units minus the recovery factor times the number of fatigued motor units.

• Un-normalize the state derivatives and store them in the State object, so they can be accessed by the integrator.

• Store the muscle force and other state-dependent variables so they can be accessed by the integrator, muscle analysis, or other components of the program.

computeIsometricForce()

To compute the isometric force in the muscle, we assume that the muscle has reached an equilibrium position in which the fiber velocity is zero, and the fatigue and recovery rates are constant. According to the Liu fatigue model, in this equilibrium position the number of active motor units is independent of the activation level (as long as the activation is greater than some threshold that depends on the fatigue and recovery rates).So we can use the fatigue and recovery rates to calculate the steady-state activation level, and then pass that activation to the computeIsometricForce function in the base class to compute the force in the muscle. 

double LiuThelen2003Muscle::computeIsometricForce(SimTK::State& s, double aActivation) const
{
  if (_optimalFiberLength < ROUNDOFF_ERROR) {
  return 0.0;
} 
 
// This muscle model includes two fatigue states, so this function
// assumes that t=infinity in order to compute the [steady-state]
// isometric force. When t=infinity, the number of active motor
// units is independent of the activation level (as long as activation
// > _recoveryFactor / (_fatigueFactor + _recoveryFactor)). So the
// passed-in activation is not used in this function (unless the fatigue
// and recovery factors are both zero which means there is no fatigue).
if ((_fatigueFactor + _recoveryFactor > 0.0) && (aActivation >=
       _recoveryFactor / (_fatigueFactor + _recoveryFactor))) {
  setActiveMotorUnits(s, _recoveryFactor / (_fatigueFactor + _recoveryFactor));
  setFatiguedMotorUnits(s, _fatigueFactor / (_fatigueFactor + _recoveryFactor));
} 

else {
  setActiveMotorUnits(s, aActivation);
  setFatiguedMotorUnits(s, 0.0);
} 
 
aActivation = getActiveMotorUnits(s); 
 
// Now you can call the base class's function with the steady-state
// activation.
return Thelen2003Muscle::computeIsometricForce(s, aActivation);
}