The steps covered in part one are:
The header file (FatigableMuscle.h)
We start defining our new muscle model by creating a header file. The model extends the Millard2012EquilibriumMuscle, and will include fatigue effects that are loosely based on 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).
The model defines active and fatigued motor unit states which are the percentages of motor units in either pool. The motor units not in the active or fatigued pools are considered "recovered" and are able to become active. The actual activation used to compute the contraction dynamics is thus a active motor unit pool percentage times the desired activation resulting from the neural excitation of the muscle. Refer to the Doxygen documentation for a complete description of the Muscle interface and details of theMillard2012EquilibriumMuscle model.
We will call our model FatigableMuscle. At the top of Fatigable Muscle.h, we include the header file for the base class and derive our new class.
FatigableMuscle will have all of the properties of the Millard2012EquilibriumMuscle model, plus two four additional ones to define the rates at which active muscle fibers fatigue and the rate at which fatigued fibers recover and defaults for the active and fatigued motor unit states.
The fatigue and recovery factors (new properties) define the behavior of the FatigableMuscle and we therefore define a constructor to initialize these properties upon creation.NOTE that the constructors also construct the properties and initialize the methods and this is typically done by calling a private constructProperties() methods.
Provide accessor methods to get and set the new properties. Note, specified properties automatically provide methods get/set_<property_name>() but we include these methods for completeness.
Millard2012EquilibriumMuscle has two states: activation and fiber length. We will add two three more states: one for the target activation, the fraction of active motor units (range 0.0 to 1.0), and one for the fraction of fatigued motor units (range 0.0 to 1.0). These states will be added to the system of equations representing the dynamics of the muscle in FatigableMuscle.cpp, but in the header we specify the accessor methods:
As a Muscle and a ModelComponent, the FatigableMuscle must satisfy the required interfaces put in place by the Muscle and ModelComponent classes from which it derives. Because our new class derives from a concrete (meaning it implements all the methods to satisfy the interfaces) Muscle class, the Millard2012EquilibriumMuscle, and not directly from Muscle, we only have to specify (override) the methods for which we want to change the behavior. We are adding new states to model fatigue and this is done by overriding the ModelComponent method addToSystem(). Methods to initialize the new states from properties and to store state values in properties are also part of the ModelComponent interface.
Finally we want to change how the Millard2012EquilibriumMuscle computes the activation (and thus force) of the muscle based on available active motor units. In this case we need to: 1) specify how the new states of the model are initialized and 2) specify the dynamics (derivatives) of all the muscle states by overriding the methods that perform these computations.
The source file (FatigableMuscle.cpp)
We implement the methods that we defined for our class in FatigableMuscle.cpp. This source file will include the methods that we need to override in Maillard2012EquilibriumMuscle that implement the fatigue behavior of the muscle.
The function constructProperties() is used to create the properties of the muscle class, specify their default values, and add them to the set of all OpenSim properties. This enables you to define them in an OpenSim model file. FatigableMuscle constructors must call this method. As we did in the header file, we need to construct 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). Two additional properties define the default values for the state variables that describe the fraction of motor units that are active and fatigued.
Adding state variables
Adding new variables to the underlying computational system of equations, whether they be continuous or discrete state variables, cache variables, or modeling options (flags) must be specified in the addToSystem method implemented by all ModelComponents. In the FatigableMuscle example, we will have the base, Millard2012EquilibriumMuscle, add its variables and then add the state variables necessary to model fatigue. Specifically, we add a targetActivation state, which represents the activation state if fatigue was not present and is equivalent to the activation state of the base class. The remaining two states are the fraction of active and fatigued motor units of the muscle.
The function computeInitialFiberEquilibrium() computes values of the muscle states assuming the muscle is in an equilibrium state. It is typically called for each muscle before beginning a dynamic simulation. In our muscle class, this function initializes the new states to reasonable values prior to simulation.
Specify the dynamics of the FatigableMuscle
With the states defined and initialized appropriately, all that remains is to specify how the states should evolve in time. The base class (Millard2012EquilibriumMuscle) already defines activation and fiber_length state variables and we added target_activation, active_motor_units and fatigued_motor_units, for a total of five muscle states and their derivatives must be provided by computeStateVariableDerivatives() in the same order. The dynamics for the target_activation state are already defined by the activation dynamics model employed by the Millard2012EquilibriumMuscle and we can use it to compute the targetActivationRate. Given that the actual activation the muscle experiences is the product of the fraction of active motor units (Ma) and the target activation (a) we can compute the activation rate of the muscle knowing the time derivative of the two quantities such that activationRate = dMa/dt*a + Ma*da/dt. The rate of conversion to/from active and fatigued motor units is based on Liu et al. 2008.