In the last tutorial we extended the single link pendulum to three links. This made for slightly more interesting dynamics but we could still only stand by helplessly while it flopped around. In this tutorial we will add muscles to the three-link pendulum. We'll add two mono-articular (crosses one joint) muscles for each rigid body, one to "flex" and one to "extend". We'll add four bi-articular (crosses two joints) muscles. Two to span the first and second joints and two more to span the second and third joints.

In a future tutorial we will learn some fancy ways to control these muscles, but for the time being we will just learn how to use Neuromechanic to "equilibrate" the model. That is, Neuromechanic will automatically choose a pattern of muscle activation such that the muscles generate forces that balance all the external forces acting on the pendulum.

## Muscle Modeling

Neuromechanic has a few different Hill-type models which are used to approximate the dynamics of muscle excitation and contraction. We're not going to get into the theory of muscle modeling here, but I will summarize the muscle model we will use throughout the tutorial. The models are all point-to-point muscle models, meaning that the muscle is modeled as a linear element that traverses through a sequence of defined points.

The muscle model we will use (Zajac) has an inelastic tendon and no series elastic element. The force generated by the muscle is a function of three time dependent quantities (state variables), six quantities which are independent of time (parameters), and three curves which we will define using cubic splines.

Variable description Symbol
The length of the muscle $L_M$
The velocity of contraction of the muscle $\dot{L}_M$
The activation level of the muscle $a$
Parameter description Symbol Body File Tag
The pennation angle $α$ <PennationAngle>
The optimal fiber length $L_F^0$ <OptimalFiberLength>
The tendon slack length $L_{TS}$ <TendonSlackLength>
The passive viscosity $η$ <PassiveDamping>
The maximum muscle contraction timescale $τ_{max}$ <Timescale>
The maximum isometric force of the muscle $F_{max}$ <MaxForce>
Curve description Symbol Body File Tag
The active force - length relationship $afl$ <ActiveForceLength>
The passive force - length relationship $pfl$ <PassiveForceLength>
The force - velocity relationship $fv$ <VelocityForceLength>

These parameters and variables are used to define

• The fiber length $L_F=\cfrac{L_M-L_{TS}}{cos \left(α\right)}$
• The fiber velocity $\dot{L}_F=\dot{L}_M$
• The maximum shortening velocity $V_F^{max}=\cfrac{L_F^0}{τ_{max}}$

And together, all of this gives us the force generated by the muscle

$F_{M} = F_{max} cos \left(α\right)\left( afl \left( \cfrac{L_F}{L_F^0} \right) fv \left( \cfrac{\dot{L}_F}{V_F^{max}} \right) a + pfl \left( \cfrac{L_F}{L_F^0} \right) + η \cfrac{\dot{L}_F}{V_F^{max}} \right)$

The parameters and curves described above will need to be defined in the body file. We will also have to define a source of muscle activation. In addition we need to define the muscle path with a sequence of points.

## The active force length, passive force length, and active force velocity curves

The $afl$, $pfl$ and $fv$ curves are defined in the <Functions> element. Below we show how to populate <Spline> elements to define active/passive force-length and force-velocity curves. As shown in the figure, each pair of numbers in the <ControlPoints> element give the interpolating values for the spline. These functions will be called later when we define muscle <Model> elements in the <FunctionReferences> section.

<NeuromechanicFile>
...
<Functions>
<Spline Name="activeforcelengthcurve" Type="NaturalCubic">
<ConstraintA>0.</ConstraintA>
<ConstraintB>0.</ConstraintB>
<ControlPoints>-5. 0.
0. 0.
4.01E-1 0.
4.02E-1 0.
4.035E-1 0.
5.2725E-1 2.2667E-1
6.2875E-1 6.3667E-1
7.1875E-1 8.5667E-1
8.6125E-1 9.5E-1
1.045 9.9333E-1
1.2175 7.7E-1
1.4387 2.4667E-1
1.6187 0.
1.62 0.
1.621 0.
2.2 0.
5. 0.
</ControlPoints>
</Spline>
<Spline Name="passiveforcelengthcurve" Type="NaturalCubic">
<ConstraintA>0.</ConstraintA>
<ConstraintB>0.</ConstraintB>
<ControlPoints>
-5. 0.
9.98E-1 0.
9.99E-1 0.
1. 0.
1.1 3.5E-2
1.2 1.2E-1
1.3 2.6E-1
1.4 5.5E-1
1.5 1.17
1.6 2.
1.601 2.
1.602 2.
5. 2.
</ControlPoints>
</Spline>
<Spline Name="velocityforcecurve" Type="NaturalCubic">
<ConstraintA>0.</ConstraintA>
<ConstraintB>0.</ConstraintB>
<ControlPoints>
-1. 0.
-9.5E-1 1.0417E-2
-9.E-1 2.1739E-2
-8.5E-1 3.4091E-2
-8.E-1 4.7619E-2
-7.5E-1 6.25E-2
-7.E-1 7.8947E-2
-6.5E-1 9.7222E-2
-6.E-1 1.1765E-1
-5.5E-1 1.4062E-1
-5.E-1 1.6667E-1
-4.5E-1 1.9643E-1
-4.E-1 2.3077E-1
-3.5E-1 2.7083E-1
-3.E-1 3.1818E-1
-2.5E-1 3.75E-1
-2.E-1 4.4444E-1
-1.5E-1 5.3125E-1
-1.E-1 6.4286E-1
-5.E-2 7.9167E-1
0. 1.
5.E-2 1.482
1.E-1 1.6016
1.5E-1 1.6558
2.E-1 1.6867
2.5E-1 1.7068
3.E-1 1.7208
3.5E-1 1.7311
4.E-1 1.7391
4.5E-1 1.7454
5.E-1 1.7505
5.5E-1 1.7547
6.E-1 1.7583
6.5E-1 1.7614
7.E-1 1.764
7.5E-1 1.7663
8.E-1 1.7683
8.5E-1 1.7701
9.E-1 1.7717
9.5E-1 1.7732
1. 1.7745
</ControlPoints>
</Spline>
</Functions>
...
<NeuromechanicFile>

## The Motor Neuron

All neurons (including motor neurons) are defined in the <Neurons> element of the body file with each individual neuron being defined in its own <Neuron> element. We'll save the discussion of <Neuron> elements for another tutorial. For now, make your first neuron look like this:

<NeuromechanicFile>
<Neurons>
<EqCCost>0.</EqCCost>
<EqBounds>1.E-2 9.5E-1</EqBounds>
<OutputBounds>0. 1.</OutputBounds>
<Synapse Name="synapse1" Type="constant">0.11 0.</Synapse>
</Neuron>
</Neurons>
</NeuromechanicFile>

The output of this particular neuron will be the first value of the <Synapse> element (0.11) and is used as the activation of the corresponding muscle. Any value outside the range set by <OutputBounds> will be saturated. For example, if the first element of the <Synapse> element was 1.1 the output of the motor neuron (and the activation of the corresponding muscle) would be 1.

The <EqBounds> and <EqCCost> elements are used during equilibration and will be discussed in the Equilibration section below.

In the next section we will be defining muscles. After you learn how to define a muscle, return to this section and make sure that there is a neuron for each muscle that we will define in the <Muscles> section. For motor neurons the Name attribute of the <Neuron> should match the contents of the corresponding <MotorNeuron> element.

## Defining muscle parameters in the <Muscles> section of the body file

All muscles are defined in the <Muscles> element of the body file with each individual muscle being defined in its own <Muscle> element. The contents of <Muscle> depend on the muscle model used. For the muscle model in this tutorial the <Muscle> element looks like this:

<NeuromechanicFile>
<Muscles>
...
<MaxForce>80</MaxForce>
<Model Type="Zajac">
<TendonSlackLength>.01</TendonSlackLength>
<OptimalFiberLength>.11</OptimalFiberLength>
<PennationAngle>0.</PennationAngle>
<PassiveDamping>1.E-4</PassiveDamping>
<Timescale>.1</Timescale>
<FunctionReferences>
<ActiveForceLength>activeforcelengthcurve</ActiveForceLength>
<PassiveForceLength>passiveforcelengthcurve</PassiveForceLength>
<ActiveForceVelocity>velocityforcecurve</ActiveForceVelocity>
</FunctionReferences>
<MusclePath>
<Point Parent="base">-.03 0 0</Point>
</MusclePath>
</Model>
</Muscle>
...
</Muscles>
</NeuromechanicFile>

Our first muscle is defined on lines starting and ending with <Muscle>. Each muscle is given a unique Name attribute. Next, each muscle must have an associated <MotorNeuron> which supplies the excitation (and in this case activation) characteristics of the muscle. Notice that the name matches the <Neuron> we specified in defining a neuron above. The next element <MaxForce> specifies the maximum isometric force of the muscle.

The <MotorNeuron>, <MaxForce>, and Name attributes are required for every muscle. Parameters specific to the type of muscle model are specified use with the <Model> element. <Model> has an attribute Type to designate the model type which is, in this case, "Zajac" (the Hill-type model with no series elasticity and pennation angle used here). Within the <Model> element we define the necessary parameters for the model, as mentioned in the table above. We also supply the model with functions that relate the time dependent variables to the output force using references to spline functions using <FunctionReferences>. The values of these elements (e.g. <ActiveForceLength>) are strings referring to the name of a function defined in the <Functions> section.

In addition, we define the line of action of the muscles by a sequence of <Point> elements defined within the <MusclePath> element. Each <Point> must have a Parent attribute which is the name of the <RigidBody> in which the point is defined in.

Ok, so lets add muscles to the model! Create a <Muscles> and <Neurons> element as children of the <NeuromechanicFile> element. Copy the muscle above so that you have six muscles in the <Muscles> section and change the Name, MotorNeuron and MusclePath as shown in the table.

Name MotorNeuron First Point Second Point

Also, make copies of the neuron you specified in the neuron above so that you have six neurons in the <Neurons> section. Each neuron should have a name that corresponds to each of the <MotorNeuron> references in the <Muscle> elements.

Finally we'll add some biarticular muscles

Name MotorNeuron OptimalFiberLength First Point Second Point

Our models are getting a bit large to embed here in the tutorial but if you get stuck you can download the file here.

## Simulation and Equilibration

Load the file into Neuromechanic and set up the view to suit. Muscles are rendered as lines unless the preferences "MusclePathVisible" is set to false in Preferences. If "MusclePathVisible" is false then the muscles are rendered with a fusiform shape. The maximum cross-sectional area of the muscle is proportional to the maximum isometric force of the muscle (<MaxForce>). The proportionality constant is the specific tension (<MuscleSpecificTension>) of muscle and is set in the <Parameters> section of the body file. By default this value is 225000. If you are using the MKS system then this is 22.5N/cm^2.

 Edit joint angles Set equilibration type

• Adjust the joint angles of the pendulum
• Run a simulation. Notice that the pendulum with muscles falls very differently than without muscles.
• After the simulation click the reset button on the simulation toolbar
• Expand parameters in the ModelTree and click "EquilibrationType"

Equilibration is the process of choosing a muscle activation pattern that will balance the external forces acting on the model. Equilibration is accomplished by solving a quadratic programming problem with respect to muscle forces. The quadratic programming problem has linear equality constraints (the part that actually balances forces) and linear inequality constraints (upper and lower bounds on each muscle's force). The equilibration editor lets you adjust parameters that affect the terms of the quadratic programming problem.

First the <EquilibrationType> element in the <Parameters> section can be set to one of three values in the body file or selected in the drop-down menu at the top of the equilibration editor.

• If it is set to "MinimumActivation" then the quadratic term of the cost function is a diagonal matrix with diagonal terms = 1/UB^2. UB is the upper bound of the corresponding muscle.
• If it is set to "MinimumForce" then the quadratic term of the cost function is an identity matrix
• If it is set to "UserDefined" then all elements of the quadratic term of the cost function are zero except those explicitly defined in <EqQCost> elements in each motor neuron element.

The elements of the linear term of the cost function are set in <EqCCost> elements in each motor neuron element and can be adjusted in the equilibration editor.

The lower and upper bounds of the inequality constraint are set in <EqBounds> elements in each motor neuron element and can also be adjusted in the equilibration editor.

The quadratic problem is solved with respect to muscle force. However, each term in an <EqCCost> or <EqBounds> element is in units of the neuron output. These are converted internally to units of force before the quadratic programming problem is solved. Each <EqQCost> term is assumed to be unitless. Therefore the cost function is in units of Force^2.

In our model the <EqCCost> has been set to zero for each muscle and <EqBounds> are set to 1% and 95% of neuron output.

• From the menu select Tools->Equilibrate (shortcut key F6)

In the log window you will see the following output (I used 20 degrees for each unlocked degree of freedom to get these results). First the imbalance at each state before the equilibration is shown. In this case the states are the velocity and acceleration of the 9 degrees of freedom. Velocities are listed first and are zero. If these are non-zero the equilibration will proceed and try to equilibrate the accelerations but the velocities will still be non-zero. So, for a true equilibration the velocities need to be zeroed before equilibrating.

After the velocities the accelerations are listed. Notice that the values are non-zero for the three unlocked degrees of freedom.

 Index State                                  Value        Locked
==== ======================================= ============ ===
1                                     base1tx  0.000E+000 *
2                                     base1ty  0.000E+000 *
3                                     base1tz  0.000E+000 *
4                                     base1rx  0.000E+000 *
5                                     base1ry  0.000E+000 *
6                                     base1rz  0.000E+000 *
10                                    base1tx  0.000E+000 *
11                                    base1ty  0.000E+000 *
12                                    base1tz  0.000E+000 *
13                                    base1rx  0.000E+000 *
14                                    base1ry  0.000E+000 *
15                                    base1rz  0.000E+000 *
=====================================================================

Next all the terms of the quadratic programming problem (except the quadratic term Q) and the solution are printed.

********************************************************************************
********************************************************************************
Minimize:           f = 0.5*([x-x0]'*Q*[x-x0]) + c'*[x-x0]
Subject to:         Aeq*x = beq
LB <= x <= UB
Muscle ------------------------------------- C ------- LB ------ UB ------ X ------
1                             MuscleLink1A     0.000     0.787    74.800     0.787
2                             MuscleLink1B     3.944     4.689    74.668     4.689
3                             MuscleLink2A     0.000     0.787    74.800    14.606
4                             MuscleLink2B     3.944     4.689    74.668     4.689
5                             MuscleLink3A     0.000     0.787    74.800     0.787
6                             MuscleLink3B     3.944     4.689    74.668     4.689
7                            MuscleLink12A     0.000     0.777    73.852     0.777
8                            MuscleLink12B     2.268     3.040    75.554    28.905
9                            MuscleLink23A     0.000     0.777    73.852     0.777
10                           MuscleLink23B     2.268     3.040    75.554     3.449
------------ F ------------
f =     128.9606
------------ AEQ' ------------
------------  BEQ ------------
7.52927    -39.05015     11.36277
------------ AEQ*X-BEQ ------------
0.00000      0.00000      0.00000

Next the states are printed again. If the equilibration was successful all values should be zero (or very small)

Attempting to set muscles and neurons to equilibrium...
Index State                                  Value        Locked
==== ======================================= ============ ===
1                                    base1tx  0.000E+000 *
2                                    base1ty  0.000E+000 *
3                                    base1tz  0.000E+000 *
4                                    base1rx  0.000E+000 *
5                                    base1ry  0.000E+000 *
6                                    base1rz  0.000E+000 *
10                                   base1tx  0.000E+000 *
11                                   base1ty  0.000E+000 *
12                                   base1tz  0.000E+000 *
13                                   base1rx  0.000E+000 *
14                                   base1ry  0.000E+000 *
15                                   base1rz  0.000E+000 *
=====================================================================

Finally the model is linearized at the equilibrium point and a stability analysis is performed. The eigenvalues, doubling (or halving) time, period, damping ratio, and undamped natural frequency are printed for each mode.

________________________________________________________________________________
********************************************************************************
************************** STABILITY ANALYSIS **********************************
********************************************************************************
There are 0 rigid body modes                 (RBM)
0 pure oscillatory pairs        (POM)
2 stable non-oscillatory modes  (SNM)
1 stable oscillatory pairs      (SOM)
2 UNstable non-oscillatory modes(UNM)
0 UNstable oscillatory pairs    (UOM)
________________________________________________________________________________
TYPE    EIGENVALUE          DoublingTime      Period   DampingRatio UndmpdNatFre
________________________________________________________________________________
UNM  4.10E+000 +  0.00E+000i   1.69E-001   Infinity   -1.00E+000   4.10E+000
UNM  7.62E-001 +  0.00E+000i   9.10E-001   Infinity   -1.00E+000   7.62E-001
SNM -4.82E+000 +  0.00E+000i  -1.44E-001   Infinity    1.00E+000   4.82E+000
SOM -6.41E+000 +  1.31E+001i  -1.08E-001   4.80E-001   4.40E-001   1.46E+001
SOM -6.41E+000 + -1.31E+001i  -1.08E-001   4.80E-001   4.40E-001   1.46E+001
SNM -1.94E+002 +  0.00E+000i  -3.57E-003   Infinity    1.00E+000   1.94E+002
________________________________________________________________________________

## Visualizing muscle activation and adjusting the equilibration

Muscle saturation preferences

There is a tool to help visualize the activity level of muscles. You can set an "on" and "off" color of muscles and the shading of the muscles will vary between the two colors 1. Preferences->Rendering Tab->Muscles->MuscleColor and 2. Preferences->Rendering Tab->Muscles->SaturatedForceColor or Preferences->Rendering Tab->Muscles->SaturatedActivationColor. If the SaturatedForceColor is selected then the forces corresponding to the two colors are set in Preferences->Rendering Tab->Muscles->LowForceSaturation and Preferences->Rendering Tab->Muscles->HighForceSaturation.

• Make MuscleColor red
• Make SaturatedForceColor cyan
• Make sure SaturatedForceColor check box is selected
• Make LowForceSaturation 0
• Make HighForceSaturation 40

Notice from the color of the muscles that one of the biarticular muscles has the greatest force based on the last equilibration (Force values of the solution are the last column here)

Muscle ------------------------------------- C ------- LB ------ UB ------ X ------
1                             MuscleLink1A     0.000     0.787    74.800     0.787
2                             MuscleLink1B     3.944     4.689    74.668     4.689
3                             MuscleLink2A     0.000     0.787    74.800    14.606
4                             MuscleLink2B     3.944     4.689    74.668     4.689
5                             MuscleLink3A     0.000     0.787    74.800     0.787
6                             MuscleLink3B     3.944     4.689    74.668     4.689
7                            MuscleLink12A     0.000     0.777    73.852     0.777
8                            MuscleLink12B     2.268     3.040    75.554    28.905
9                            MuscleLink23A     0.000     0.777    73.852     0.777
10                           MuscleLink23B     2.268     3.040    75.554     3.449

Let's try equilibrating again after adjusting some parameters of the quadratic programming problem.

• Make the Linear cost of all the bi-articular muscles 1 (penalizing the activation of those muscles)
• Equilibrate and watch what happens to the muscle colors. I get the following new solution.
  Muscle ------------------------------------- C ------- LB ------ UB ------ X ------
1                             MuscleLink1A     0.000     0.787    74.800     0.787
2                             MuscleLink1B     3.944     4.689    74.668    33.415
3                             MuscleLink2A     0.000     0.787    74.800     0.787
4                             MuscleLink2B     3.944     4.689    74.668    17.614
5                             MuscleLink3A     0.000     0.787    74.800     0.787
6                             MuscleLink3B     3.944     4.689    74.668     5.143
7                            MuscleLink12A    77.739     0.777    73.852     0.777
8                            MuscleLink12B    79.411     3.040    75.554     3.040
9                            MuscleLink23A    77.739     0.777    73.852     0.777
10                            MuscleLink23B    79.411     3.040    75.554     3.040
• Run a simulation. The pendulum is in equilibrium and shouldn't move