From output of #5 write the actual library

• Single value activation
• Activation array
• Motor neuron equations
• Motor neuron fatigue
• Fiber equations
• Fiber fatigue
• Default parameters
• Configurable parameters

This figure appears to show a relationship between the normalized firing rate and the normalized force (as a percentage of maximum normalized force) for each motor unit. The figure caption states that this was modeled to be the same for each motor unit.

However when I plot this relationship for any motor unit in my code the range below 0.4 is almost totally absent.

The values input to the equation to normalize the firing rates jump from 0 to 8 imp/s. This means that a large portion of normalized force range is missing for most units.

e.g. For the first motor unit the minimum normalized firing rate is 0.72, well above the 0.4 threshold in the normalized force function.

Implement equation from "Fatigue-related changes in motor unit force capacity" section of Potvin 2017

• Does this make sense for the whole library API?
• Should we use validators for complex input?
• Should we validate the config yaml?

• What is the intent of the API
• What are the goals?
• How would someone use the library?
• How would they contribute?
• Versioning Plan

What are the run times of equivalent simulations in:

• Unvectorized Matlab Code

Run time: 1372.633170 seconds (~22 minutes)

• Vectorized Matlab Code

19.214544 seconds (~70x faster)

• Library Code

3.785 seconds (~360x faster)

All numbers are from a 50% excitation level 10 "second" run. (fthtime = 10)

• Bianca
• Potvin
• Will?

Design Guidelines

1. Look at the API from the eyes of a naive user

The most important thing is to have empathy with your primary target audience and constantly ask yourself what makes most sense for them in their perspective. This applies to everything, from the first moments reading the README.md examples to the corner-case use cases of the API. The Principle of Least Surprise applies strongly here.

2. Define the purpose and level of abstraction of your library

What problem does it solve?

It enables a more accurate and dynamic model of the relationship between control signals and force output. It manages the stateful, time-dependent relationship between control signals and changes to motor units. This can be coupled with a simulation of the muscle body to capture the full control signal to force output relationship.

What is the level of abstraction?

The collection of motor units which form a muscle.

What problems does it not solve?

• Simulating the muscle body.
• Proprioception / feedback

3. Tradeoffs always exist, but unimportant disadvantages also exist

What are the major tradeoffs?

Accuracy vs Speed.

The goal is to capture the complexity of the problem that biological creatures are faced with, while being fast enough that useful interesting research can be done quickly.

Pros

More accurate.
More dynamic.

Cons

Slower than a purely linear relationship.
Complexity might not be needed.

4. Make the correct way easy, make the wrong way possible

Keep the library user comfortable in using your library in the right way, while making him cautious in using the library in the wrong way.

4.1. Short names

4.2. Common words

4.3. Long names for 'incorrect' way of doing things.

e.g. dangerouslySetInnerHTML

5. Be always mindful of the types

• Mock out types in initial design

So, write the type signatures.

What kind of composability strategy you are using?

OOP

6. Cooperate with the host language

a.k.a. stay Pythonic

7. Be open to feedback

Guidelines by Andre Statlz

Output data from the library should match the below exactly for the same experimental parameters.

• Pipfiles
• Pip install
• pipenv install
• setup.py install
• Releases
• Stake out name on pypi
• Upload wheel
• Verify minimal-chart-output.png is real

Sources / Summaries

Thelen 2003 - Complete OpenSim Implementation Description

https://simtk-confluence.stanford.edu/download/attachments/2624181/CompleteDescriptionOfTheThelen2003MuscleModel.pdf?version=1&modificationDate=1319838594036&api=v2

Seow

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3840917/

"data also provide justification for the use of a hyperbolic equation not as a mere empirical description but as a meaningful explanation for force–velocity behavior based on cross-bridge kinetics"

"During an isotonic quick release, the muscle is suddenly released from its isometric force to a lower and constant force (i.e., isotonic load). In response to the sudden change in load, the muscle shortens in a characteristic fashion"

Characteristic Force - Velocity Curve

NOTE Immediately upon decreasing the force there is a small (ms) window of time where the muscle is acting like a damped (not quite critically damped) spring. Then contraction begins taking up the slack

The standard model of force velocity relationships is not perfect

"Careful measurements of force–velocity properties of single skeletal muscle fibers have revealed that at low loads (less than ∼5% Fmax), the measured velocities exceed those predicted by the Hill hyperbola. At the extrapolated zero load, Hill’s equation underestimates the value of Vmax by ∼6–7% (Edman et al., 1976; Julian et al., 1986). For whole muscle preparations with mixed fiber types, the underestimation of Vmax by the Hill hyperbola is found to be much greater (Claflin and Faulkner, 1989)"

Summary

The article generally supports the use of the standard Hill equation for the force velocity relationship. It goes on to discuss how this relationship might arise from the properties of cross-bridge creation, detachment, and force production under varying shortening velocities. Empirical data from Piazzesi et al. (2007) are used to back up this model.

Fatigue - The curvature of the force-velocity relationship increases with fatigue. "there must be a substantial decrease in the rate constant for attachment in the Huxley (1957) model to account for the observed decrease in power and increase in force–velocity curvature in fatigued muscle."

Clay Anderson (and others) BME 599

http://rrg.utk.edu/resources/BME599/assignments/BME599_lab_1_dynamic_simulation_jumping.pdf

MuscleModelEquations.pdf

Jovanović et al

Desmos comparison of Jovanović and Anderson

http://www.doiserbia.nb.rs/img/doi/1451-4869/2015/1451-48691501053J.pdf

Rosario et al.

Rosario et al. - 2016 - Muscle-spring dynamics in time-limited, elastic movements.pdf

"muscle contractile force is transmitted to skeletal structures through elastic structures, inextricably coupling muscle and spring dynamics." - GOOD QUOTE

"muscle force declines with contraction velocity" - References Hill 1938

"in situations where rapid spring-loading is beneficial (e.g. escape jumps and predatory ambushes), organisms may not have enough time to fully load their springs before the onset of movement. Although these organisms are not generating maximal muscle force, it is possible that their muscle – spring prop- erties maximize elastic energy storage for submaximal force production."

"We simulated dynamics within muscle–spring systems by con- necting, in series, a model of a muscle to a model of a spring ... Specifically, we connected a Hill-type muscle to a Hookean spring [1,5,16]."

Doesn't a hill-type model already have an in-series spring to account for tendon forces?

Equations used and a discussion of the parameters

Desmos graph including the supplied parameter values

The graph at lengths greater than 1 muscle-length never falls back to zero. Why?

Contrast with this graph or this graph which show force of the contractile portion going to zero as it nears 2 muscle-lengths

Summary

The regime we're targeting is much closer to the fast acting frog regime than slow-force-building grasshopper regime. Once satisfactory equations and paramaters for the force length and force velocity relationships have been established, the standard charts that take into account passive component force generation at muscle lengths >1 should be sufficient for our purposes.

Jones 2010

"There are three factors contributing to the loss of power in mammalian muscle at physiological temperatures, a decrease in isometric force, which mainly indicates a reduction in the number of active cross bridges, a slowing of the maximum velocity of unloaded shortening and an increased curvature of the force–velocity relationship."

The recovery of force in muscle takes minutes and follows an asymptotic return to normal

"is that it is the rate constant for attachment, f, rather than g1, which decreases with fatigue, giving rise to the change in curvature and being a major factor in the loss of power."

READ: Changes in the force–velocity relationship of fatigued muscle: implications for power production and possible causes

Summary

The paper discusses the metabolic and actin-myosin level changes that may occur during muscle fatigue. While interesting these factors are likely too specific and low level to be incorporated into a PyMuscle model.

In vitro force measurements for single fibers

https://www.jove.com/video/52695/measurement-maximum-isometric-force-generated-permeabilized-skeletal

Physics of the Body (older textbook)

https://www.medicalphysics.org/documents/WebPOTB.pdf

Muscles contract 15-20% of their length.

Curling a 44lb weight requires ~330 lb of force.

"A trained individual could probably curl about 200 N (~44 lb) requiring
the biceps to provide 1500 N (~330 lb) force."

Why Animals Have Different Muscle Fiber Types - Rome 1988

Kami Link

Design Guidelines

Why a second draft?

The initial API was not well aligned with the API for actuators in simulation libraries. The MuJoCo muscles were recently released are partialy defined by their peak force. In Box2d (the 2d physics library that Gym uses) joint motors are defined by their max torque and by additional limits. A PyMunk simple motor also allows for a max force parameter.

While a biologist or neuroscientist might think of a muscle in terms of its composition (the number of motor units) for the target audience of this library it is the functional aspect of muscles (force production) which is primary.

The second draft of the API will focus on revising the instantiation of muscles in terms of their peak force output.

StandardMuscle

The base Muscle class will continue to be in terms of motor units as this is well aligned with the physiology literature. A new top level class StandardMuscle will be implemented to provide the new force-centric API.

What kind of information will the user provide?

Setup

The primary (optional) paramater for a StandardMuscle will be max_force. This value will be interpreted as Newtons. From this value we will calculate how many motor units a muscle would have if it could generate this force. That calculation is essentially scaling up or down from the First Dorsal Interossei which is the muscle modeled by Potvin and Fuglevand.

Runtime

The primary method on a Muscle instance is step(). The main argument should be one of the following.

1. A single value between 0.0 and 1.0 representing a percentage of maximum excitation.

This item is new and is intended as a simplification of the API.

Range Note: 1.0 - will be translated into the excitation frequency which would stimulate a fully rested muscle to produce it's maximum force output. This value will be calculated internally at the time the muscle is instantiated.

2. A single float representing a stimulus input frequency which will be passed to each motor unit
3. An array of floats representing stimulus input frequencies for each motor unit. This array must have length N where N is the number of motor units in the model.

What kind of information will the API provide?

Primary - Motor unit output at a given timestep.
Secondary - Motor unit capacity at a given timestep.

What are the logical divisions in the functionality?

Outer class

OpenAI Gym compatability
    Action space definition NOTE: Maybe we shouldn't do this? Make it generic?
    Observation space definition
Input validation
State maintenance
Motor Neuron Model
Input -> Output behavior of the motor neurons
Fatigue of the motor neurons
Recovery of motor neurons
Fiber Model
Input -> Output behavior of the muscle fibers
Fatigue of the muscle fibers
Recovery of the muscle fibers

How does a user write a plugin for this library? If they can't how do they extend it?

...

As requested by Prof Potvin.

Update API design based on #3 and #4

• 'Simple' single value activations
• 'Long but possible' array of activations

Answer the Following Questions

How should users deal with multiple / many muscles?*

For 0.1 I'm going to focus on the Muscle API and not worry about a higher level wrapper.

Should they have to send input to each one with .step()?

See above

What about inspecting each muscle for output and state?

See above

Try to do this without gym dependency.

• describe biological input - output relationship (blocked by this question
• explain units of input / output
• scope - positive and negative examples

This will help verify the fatigue implementations of both the neuron pool and the muscle fibers.

• Add ability to disable central fatigue while maintaining peripheral
• Disable central fatigue in a new muscle class that extends PotvinAndFuglevandMuscle
• Implement basic recovery in the fibers model
• Update examples to use new extended muscle class

What is the best way to map one to the other?

What experimental values can inform this relationship?

Use muscles in cartpole with Gym

Investigate and document any slow points. Collect ideas for speed ups including

• parallelization
• GPU support
• caching

• Other muscle sim libraries?

None that do exactly what we need. The fatigable muscle model by Yong and the Motor Unit Stack by Lin come close. Thelen 2003 and Millard 2012 also provide useful reference models. Most recently Kim 2015, 2018 provide very detailed models of everything from motorneuron to muscle fiber force output.

• What methods do they have?

No great insights here.

• How do they break down subtasks?

Both Motor Unit Stack and PyMus break down the chain of events as you would expect. First the motor neuron train is calculated, then the fiber trains.

• When do we need #! ?
• Author and date

After attempting #22 and #23 I discovered the model produces unexpected results. Fatigue is far too fast. I spent several hours trying to debug my code implementing the equations from the paper to no avail.

So now I need to go back to the original code and compare outputs at each stage very carefully to detect bugs.

Travis or similar.

Document Religiously

• Document every class, interface, method, constructor, parameter, and exception
  • Class: what an instance represents
  • Method: contract between method and its client
  • Preconditions, postconditions, side-effects
  • Parameter: indicate units, form, ownership
• Document state space very carefully

MuJoCO

• Get MuJoCo installed and running
• Run through an Ant sample from elsewhere
• Tweak the model to learn the format
• Read MuJoCo docs
• Add one constraint
• Add many constraints (I'm guessing it will take a minimum of 16)

Notes

Tendons

These are minimum length paths which can wrap around geometry. They are defined using the following tag hierarchy

<tendon>
  <spatial>
    <site/>
    <site/>
    ...
    <site/>
  </spatial>
</tendon>

Each <site/> tag defines a connection point for the path. If a <geom/> tag is in the list the path will wrap around that geometry.

Both sites and geometries for tendons must be defined and named in the bodies to which they are associated. The names are then used when defining a spatial tendon to calculate the path.

Generating force

Tendons can be associated with actuators directly, can transfer force along their lengths, or generate forces as springs.

A tendon transferring force along its length to multiple connection points from a weight.

As springs

To enable the spring-like properties of a tendon you need to set the stiffness attribute of the tag. In my experiments to produce noticeable force this needs to be quite high (at least 1000). In addition you need to set the damping attribute otherwise the spring will not settle. A value of 2.0 seems to have a noticeable effect.

As actuators

To use a tendon as an actuator you need to specify an actuator as follows

<actuator>
  <general tendon="{tendon name}" />
  <general tendon="{tendon name}" />
  ...
  <general tendon="{tendon name}" />
</actuator> 

By specifying a tendon the general actuator is associated with you make it into a tendon actuator. These actuators take control values and translate them into changes in length of the tendon and (I believe) output forces. Usually the length is directly proportional to the input control value but this can be changed through the gear (ratio) attribute to make it a multiple or fraction of the input. There is also a gainprm attribute which multiplies the force output. In the sample balanceBeam.xml this gain is set to 50 which seems quite high.

A tendon actuator being sent increasingly negative control values causes the ball to rise. Here -19 with a gainprm of 50.0.

Note that applying force using a tendon actuator does not change the resting length of the tendon if it is also a spring.

Applying sufficient force with a tendon actuator can cause the length of the tendon to go to zero.

The Force Value Chain in MuscledAnt sim

A StandardMuscle is instantiated with a max_force value. This is the largest force (in Newtons) that this muscle can theoretically produce. By default this value is 100N.

In reality the force it can produce is less than this value because the strength of a muscle and the Maximum Voluntary Force are different.

A StandardMuscle step() method takes a value between 0 and 1 (inclusive)

This is then passed through the current state of the muscle simulation and an output in arbitrary force units (NOT newtons) is returned.

To use this value to control the actuators in the ant sim you need to know how the muscle actuators are defined. By default they are constructed as General actuators attached to Tendons. They have a control range of -1 to 0. They have a gain parameter of 100.

If you pass a control value of -1 to one of these actuators there will be a contractile force of 100 newtons produced at that attachment sites of the tendons.

Note that the gainprm value in the General actuator and the max_force value in the StandardMuscle should be the same.

To translate from the arbitrary units of a muscle output to the control range expected by the general actuators this calculation is performed

ctrl_value = -1 * (current_output / muscle.max_output)

This is pretty confusing. Need to simplify.

Discussion

It's not clear what combination of the three factors should be included as part of the muscle control.

• Muscle length usually decreases when muscle activation increases (but not always)
• Muscle stiffness usually increases when muscle activation increases
• Force output (almost by definition) increases when activation increase outside of fatigue

Thanks for the package :)

At the current version, the example pymuk-gym-example doesn't run because of a small type in pymuk-arm.py.

In line 120 and line 108 the "pymunk.constraint" should be replaced by "pymunk.constraints"