Custom Muscle Model Part Two

The steps covered in part two are:

Example program

To illustrate the new muscle model, we will modify the tug-of-war example described in Performing a Simulation so that a LiuThelen2003Muscle pulls against a Thelen2003Muscle. We will run two five-second forward dynamics simulations, one with the fatigue and recovery parameters set to zero, and one with them set to abnormally high values, to exaggerate the fatigue effect. We will quantify the fatigue effect by plotting the lengths of the two muscles during the simulations.

Modified tug-of-war example

The source file mainFatigue.cpp contains the main() function from the tug-of-war simulation, with a few modifications so we can use it with our new muscle model. The source file in the MuscleExample folder already contains all of the necessary modifications, but here is a description of how it was changed. First, we want to lock the rotational degrees of freedom in the ground->block joint so that the block does not twist as the muscles pull on it (twisting interferes with the straight back-and-forth length changes).

// Free joint states
CoordinateSet &coordinates = osimModel.updCoordinateSet();
coordinates[0].setValue(si, 0, true);
coordinates[1].setValue(si, 0, true);
coordinates[2].setValue(si, 0, true);
coordinates[3].setValue(si, 0, true);
coordinates[4].setValue(si, 0, true);
coordinates[5].setValue(si, 0, true);
coordinates[0].setLocked(si, true);
coordinates[1].setLocked(si, true);
coordinates[2].setLocked(si, true);

Next, we want to increase the simulation time from 4.0 seconds to 5.0 to better see the fatigue effect.

// Define the initial and final simulation times
double initialTime = 0.0;
double finalTime = 5.0;

Lastly, we want to change the names of the output files. 

// Save the simulation results
Storage statesDegrees(manager.getStateStorage());
osimModel.updSimbodyEngine().convertRadiansToDegrees(statesDegrees);
statesDegrees.print("tugOfWar_fatigue_states_degrees.mot"); 
 
// Save the OpenSim model to a file
osimModel.print("tugOfWar_fatigue_model.osim"); 

Using LiuThelen2003Muscle

The final change we need to make to the original example is to use a LiuThelen2003Muscle instead of a Schutte1993Muscle. We are also going to set the maximum isometric force of the Liu muscle to a value twice that of the Thelen muscle, so that it initially pulls the block towards its side before fatiguing and letting the Thelen muscle pull it back. After including the header file, LiuThelen2003Muscle.h, at the top of mainFatigue.cpp, we can modify the code that creates the muscles to look as shown below. Note that the fatigue and recovery factors are both 0.0. We will change them to non-zero values before running the second simulation.

// Create two muscles
double maxIsometricForce1 = 4000.0, maxIsometricForce2 = 2000.0,
optimalFiberLength = 0.2, tendonSlackLength = 0.2;
double pennationAngle = 0.0, activation = 0.0001,
deactivation = 1.0, fatigueFactor = 0.0, recoveryFactor = 0.0; 
 
// muscle 1 (model with fatigue)
LiuThelen2003Muscle muscle1("Liu", maxIsometricForce1, optimalFiberLength, tendonSlackLength, pennationAngle,
    fatigueFactor, recoveryFactor);
muscle1.setActivationTimeConstant(activation);
muscle1.setDeactivationTimeConstant(deactivation);
 
// muscle 2 (model without fatigue)
Thelen2003Muscle muscle2("Thelen", maxIsometricForce2, optimalFiberLength, tendonSlackLength, pennationAngle);
muscle2.setActivationTimeConstant(activation);
muscle2.setDeactivationTimeConstant(deactivation); 
 
// Define the path of the muscles
muscle1.addNewPathPoint("Liu-point1", ground, SimTK::Vec3(0.0,0.05,-0.35));
muscle1.addNewPathPoint("Liu-point2", block, SimTK::Vec3(0.0,0.0,-0.05)); 
muscle2.addNewPathPoint("Thelen-point1", ground, SimTK::Vec3(0.0,0.05,0.35));
muscle2.addNewPathPoint("Thelen-point2", block, SimTK::Vec3(0.0,0.0,0.05));

We also want to change the activation controls for the muscles so that they both start inactivated and ramp up to full activation at the same rate.

// Define the initial and final controls
double initialControl[2] = {0.0, 0.0};
double finalControl[2] = {1.0, 1.0}; 

Analyzing the fatigue effect

After compiling and running the simulation, we can load the model into the OpenSim GUI and plot muscle lengths to analyze the fatigue effect. The model that is used by the example program is saved to the file tugOfWar_fatigue_model.osim. This model contains a Thelen2003 muscle and a LiuThelen2003 muscle. However, because LiuThelen2003 is a new muscle model that OpenSim does not know about, we cannot load that model into OpenSim without first creating a plug-in containing the LiuThelen2003Muscle class. Rather than create a plug-in for this example, we will instead make a small modification to the model file so that it can be loaded into OpenSim. All we need to do is change the XML tag for the muscle named "Liu" from LiuThelen2003Muscle to Thelen2003Muscle. This changes its type to a Thelen2003 muscle, meaning that it will not include any fatigue effects. But this is OK for our example. The example program still uses the correct muscle types, but the OpenSim GUI will load the modified model. The muscles in this modified model have the same paths as in the original model, so we can still plot muscle-tendon lengths to analyze the fatigue effects. This modified model is provided for you, and is in the file tugOfWar_fatigue_model_GUI.osim. So we now want to load the model tugOfWar_fatigue_model_GUI.osim into OpenSim along with the motion file with the simulation results, tugOfWar_fatigue_states_degrees.mot. Plotting muscle-tendon length for both muscles as a function of the motion should generate a plot like the one below.

Both muscles start with zero activation at a length of 0.3 meters. Because the Liu muscle is twice as strong as the Thelen muscle, as they both develop force the Liu muscle "wins" the tug-of-war and stretches the Thelen muscle.

We can now change the fatigue and recovery parameters and see what difference this makes in the tug-of-war. Normal values for the parameters are fatigue = 0.02 and recovery = 0.008. But we will exaggerate the effect for this example and use 0.1 and 0.04, respectively.

// Create two muscles
double maxIsometricForce1 = 4000.0, maxIsometricForce2 = 2000.0, optimalFiberLength = 0.2, tendonSlackLength = 0.2;
double pennationAngle = 0.0, activation = 0.0001, deactivation = 1.0, fatigueFactor = 0.1, recoveryFactor = 0.04;

After compiling and running a second simulation, we can load the same (overwritten) motion file into the OpenSim GUI and plot the muscle lengths as a function of the new motion, adding the curves to the first plot. The result looks like this:

The green and magenta curves show the Liu and Thelen muscle lengths, respectively, during the second simulation. In this case, the Liu muscle is still stronger initially, and stretches the Thelen muscle. But it quickly starts to fatigue, and at about one second, has weakened enough so that the Thelen muscle begins to lengthen it.