Muscle simulation¶

Marcos Duarte

Let's simulate the 3-component Hill-type muscle model we described in Muscle modeling and illustrated below:

Pennation angle¶

The pennation angle will vary during muscle activation; for instance, Kawakami et al. (1998) showed that the pennation angle of the medial gastrocnemius muscle can vary from 22$^o$ to 67$^o$ during activation. The most used approach is to assume that the muscle width (defined as the length of the perpendicular line between the lines of the muscle origin and insertion) remains constant (Scott & Winter, 1991):

$$w = L_{M,0} \sin\alpha_0$$

The pennation angle as a function of time will be given by:

$$\alpha = \sin^{-1} \left(\frac{w}{L_M}\right)$$

The cosine of the pennation angle can be given by (if $L_M$ is known):

$$\cos \alpha = \frac{\sqrt{L_M^2-w^2}}{L_M} = \sqrt{1-\left(\frac{w}{L_M}\right)^2}$$

or (if $L_M$ is not known):

$$\cos \alpha = \frac{L_{MT}-L_T}{L_M} = \frac{1}{\sqrt{1 + \left(\frac{w}{L_{MT}-L_T}\right)^2}}$$

Muscle force¶

In general, the dependence of the force of the contractile element with its length and velocity and with the activation level are assumed independent of each other:

$$F_{CE}(a, L_{CE}, \dot{L}_{CE}) = a \: f_l(L_{CE}) \: f_v(\dot{L}_{CE}) \: F_{M0}$$

where $f_l(L_M)$ and $f_v(\dot{L}_M)$ are mathematical functions describing the force-length and force-velocity relationships of the contractile element (typically these functions are normalized by $F_{M0}$, the maximum isometric (at zero velocity) muscle force, so we have to multiply the right side of the equation by $F_{M0}$).

And for the muscle force:

$$F_{M}(a, L_M, \dot{L}_M) = \left[a \: f_l(L_M)f_v(\dot{L}_M) + F_{PE}(L_M)\right]F_{M0}$$

This equation for the muscle force, with $a$, $L_{M}$, and $\dot{L}_{M}$ as state variables, can be used to simulate the dynamics of a muscle given an excitation and determine the muscle force and length. We can rearrange the equation, invert the expression for $f_v$, and integrate the resulting first-order ordinary differential equation (ODE) to obatin $L_M$:

$$\dot{L}_M = f_v^{-1}\left(\frac{F_{SE}(L_{MT}-L_M\cos\alpha)/\cos\alpha - F_{PE}(L_M)}{a f_l(L_M)}\right)$$

This approach is the most commonly employed in the literature (see for example, OpenSim; McLean, Su, van den Bogert, 2003; Thelen, 2003; Nigg and Herzog, 2007).

Although the equation for the muscle force doesn't have numerical singularities, the differential equation for muscle velocity has four (OpenSim Millard 2012 Muscle Models):
When $a \rightarrow 0$; when $f_l(L_M) \rightarrow 0$; when $\alpha \rightarrow \pi/2$; and when $\partial f_v/\partial v \rightarrow 0$.
The following solutions can be employed to avoid the numerical singularities (OpenSim Millard 2012 Muscle Models):
Adopt a minimum value for $a$; e.g., $a_{min}=0.01$; adopt a minimum value for $f_l(L_M)$; e.g., $f_l(0.1)$; adopt a maximum value for pennation angle; e.g., constrain $\alpha$ to $\cos\alpha > 0.1 \; (\alpha < 84.26^o)$; and make the slope of $f_V$ at and beyond maximum velocity different than zero (for both concentric and excentric activations). We will adopt these solutions to avoid singularities in the simulation of muscle mechanics. A probem of imposing values to variables as described above is that we can make the ordinary differential equation numerically stiff, which will increase the computational cost of the numerical integration. A better solution would be to modify the model to not have these singularities (see OpenSim Millard 2012 Muscle Models).

Simulation¶

Let's simulate muscle dynamics using the Thelen2003Muscle model we defined in Muscle modeling. For the simulation of the Thelen2003Muscle, we simply have to integrate the equation:

$$V_M = (0.25+0.75a)\,V_{Mmax}\frac{\bar{F}_M-a\bar{f}_{l,CE}}{b}$$

where

$$b = \left\{ \begin{array}{l l l} a\bar{f}_{l,CE} + \bar{F}_M/A_f \quad & \text{if} \quad \bar{F}_M \leq a\bar{f}_{l,CE} & \text{(shortening)} \\ \\ \frac{(2+2/A_f)(a\bar{f}_{l,CE}\bar{f}_{CEmax} - \bar{F}_M)}{\bar{f}_{CEmax}-1} \quad & \text{if} \quad \bar{F}_M > a\bar{f}_{l,CE} & \text{(lengthening)} \end{array} \right.$$

The equation above already contains the terms for actvation, $a$, and force-length dependence, $\bar{f}_{l,CE}$. The equation is too complicated for solving analytically, we will solve it by numerical integration using the scipy.integrate.ode class of numeric integrators, particularly the dopri5, an explicit runge-kutta method of order (4)5 due to Dormand and Prince (a.k.a. ode45 in Matlab). We could run a simulation using OpenSim; it would be faster, but for fun, let's program in Python. All the necessary functions for the Thelen2003Muscle model described in Muscle modeling were grouped in one file (module), muscles.py. Besides these functions, the module muscles.py contains a function for the muscle velocity, vm_eq, which will be called by the function that specifies the numerical integration, lm_sol; a standard way of performing numerical integration in scientific computing:

def vm_eq(self, t, lm, lm0, lmt0, lmopt, ltslack, alpha0, vmmax, fm0):
"""Equation for muscle velocity."""
if lm < 0.1*lmopt:
lm = 0.1*lmopt
a     = self.activation(t)
lmt   = self.lmt_eq(t, lmt0)
alpha = self.penn_ang(lmt=lmt, lm=lm, lm0=lm0, alpha0=alpha0)
lt    = lmt - lm*np.cos(alpha)
fse   = self.force_se(lt=lt, ltslack=ltslack)
fpe   = self.force_pe(lm=lm/lmopt)
fl    = self.force_l(lm=lm/lmopt)
fce_t = fse/np.cos(alpha) - fpe
vm    = self.velo_fm(fm=fce_t, a=a, fl=fl)
return vm

def lm_sol(self, fun, t0, t1, lm0, lmt0, ltslack, lmopt, alpha0, vmmax, fm0, show, axs):
"""Runge-Kutta (4)5 ODE solver for muscle length."""
if fun is None:
fun = self.vm_eq
f = ode(fun).set_integrator('dopri5', nsteps=1, max_step=0.005, atol=1e-8)
f.set_initial_value(lm0, t0).set_f_params(lm0, lmt0, lmopt, ltslack, alpha0, vmmax, fm0)
# suppress Fortran warning
warnings.filterwarnings("ignore", category=UserWarning)
data = []
while f.t < t1:
f.integrate(t1, step=True)
d = self.calc_data(f.t, f.y, lm0, lmt0, ltslack, lmopt, alpha0, fm0)
data.append(d)
warnings.resetwarnings()
data = np.array(data)
self.lm_data = data
if show:
self.lm_plot(data, axs)
return data


muscles.py also contains some auxiliary functions for entering data and for plotting the results. Let's import the necessary Python libraries and customize the environment in order to run some simulations using muscles.py:

In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
#%matplotlib nbagg
import matplotlib
matplotlib.rcParams['lines.linewidth'] = 3
matplotlib.rcParams['font.size'] = 13
matplotlib.rcParams['lines.markersize'] = 5
matplotlib.rc('axes', grid=False, labelsize=14, titlesize=16, ymargin=0.05)
matplotlib.rc('legend', numpoints=1, fontsize=11)
# import the muscles.py module
import sys
sys.path.insert(1, r'./../functions')
import muscles


The muscles.py module contains the class Thelen2003() which has the functions we want to use. For such, we need to create an instance of this class:

In [2]:
ms = muscles.Thelen2003()


Now, we need to enter the parameters and states for the simulation: we can load files with these values or enter as input parameters when calling the function (method) 'set_parameters()' and 'set_states()'. If nothing if inputed, these methods assume that the parameters and states are stored in the files 'muscle_parameter.txt' and 'muscle_state.txt' inside the directory './../data/'. Let's use some of the parameters and states from an exercise of the chapter 4 of Nigg and Herzog (2006).

In [3]:
ms.set_parameters()
ms.set_states()

The parameters were successfully loaded and are stored in the variable P.
"lt0" value "" was replaced by NaN.
The states were successfully loaded and are stored in the variable S.


We can see the parameters and states:

In [4]:
print('Parameters:\n', ms.P)
print('States:\n', ms.S)

Parameters:
{'u_min': 0.01, 'ltslack': 0.223, 'id': '', 'lmopt': 0.093, 't_act': 0.015, 'epst0': 0.04, 'af': 0.25, 'gammal': 0.45, 'fmlen': 1.4, 'kttoe': 3.0, 'epsm0': 0.6, 'name': '', 'kpe': 5.0, 'alpha0': 0.0, 'u_max': 1.0, 'fm0': 7400.0, 't_deact': 0.05, 'vmmax': 10.0}
States:
{'lm0': 0.09, 'lmt0': 0.313, 'lt0': nan, 'id': '', 'name': ''}


We can plot the muscle-tendon forces considering these parameters and initial states:

In [5]:
ms.muscle_plot()


Let's simulate an isometric activation (and since we didn't enter an activation level, $a=1$ will be used):

In [6]:
def lmt_eq(t, lmt0):
# isometric activation

lmt = lmt0

return lmt

ms.lmt_eq = lmt_eq

In [7]:
data = ms.lm_sol()


We can input a prescribed muscle-tendon length for the simulation:

In [8]:
def lmt_eq(t, lmt0):
# prescribed change in the muscle-tendon length

if t < 1:
lmt = lmt0
if 1 <= t < 2:
lmt = lmt0 - 0.04*(t - 1)
if t >= 2:
lmt = lmt0 - 0.04

return lmt

ms.lmt_eq = lmt_eq

In [9]:
data = ms.lm_sol()


Let's simulate a pennated muscle with an angle of $30^o$. We don't need to enter all parameters again, we can change only the parameter alpha0:

In [10]:
ms.P['alpha0'] = 30*np.pi/180
print('New initial pennation angle:', ms.P['alpha0'])

New initial pennation angle: 0.5235987755982988


Because the muscle length is now shortenned by $\cos(30^o)$, we will also have to change the initial muscle-tendon length if we want to start with the tendon at its slack length:

In [11]:
ms.S['lmt0'] = ms.S['lmt0'] - ms.S['lm0'] + ms.S['lm0']*np.cos(ms.P['alpha0'])
print('New initial muscle-tendon length:', ms.S['lmt0'])

New initial muscle-tendon length: 0.300942286341

In [12]:
data = ms.lm_sol()


Here is a plot of the simulated pennation angle:

In [13]:
plt.plot(data[:, 0], data[:, 9]*180/np.pi)
plt.xlabel('Time (s)')
plt.ylabel('Pennation angle $(^o)$')
plt.show()


Change back to the old values:

In [14]:
ms.P['alpha0'] = 0
ms.S['lmt0']   = 0.313


We can change the initial states to show the role of the passive parallel element:

In [15]:
ms.S = {'id': '', 'lt0': np.nan, 'lmt0': 0.323, 'lm0': 0.10, 'name': ''}

In [16]:
ms.muscle_plot()


Let's also change the excitation signal:

In [17]:
def excitation(t, u_max=1, u_min=0.01, t0=1, t1=2):
"""Excitation signal, a hat signal."""
u = u_min
if t >= t0 and t <= t1:
u = u_max
return u

ms.excitation = excitation

In [18]:
act = ms.activation_sol()


And let's simulate an isometric contraction:

In [19]:
def lmt_eq(t, lmt0):
# isometric activation

lmt = lmt0

return lmt

ms.lmt_eq = lmt_eq

In [20]:
data = ms.lm_sol()


Let's use as excitation a train of pulses:

In [21]:
def excitation(t, u_max=.5, u_min=0.01, t0=.2, t1=2):
"""Excitation signal, a train of square pulses."""

u = u_min
ts = np.arange(1, 2.0, .1)
#ts = np.delete(ts, np.arange(2, ts.size, 3))
if t >= ts[0] and t <= ts[1]:
u = u_max
elif t >= ts[2] and t <= ts[3]:
u = u_max
elif t >= ts[4] and t <= ts[5]:
u = u_max
elif t >= ts[6] and t <= ts[7]:
u = u_max
elif t >= ts[8] and t <= ts[9]:
u = u_max

return u

ms.excitation = excitation

In [22]:
act = ms.activation_sol()

In [23]:
data = ms.lm_sol()


Module muscles.py¶

In [ ]:
# %load ./../functions/muscles.py
"""Muscle modeling and simulation."""

from __future__ import division, print_function
import numpy as np
from scipy.integrate import ode
import warnings
import configparser

__author__ = 'Marcos Duarte, https://github.com/demotu/BMC'
__version__ = 'muscles.py v.1 2015/03/01'

class Thelen2003():
""" Thelen (2003) muscle model.
"""

def __init__(self, parameters=None, states=None):
if parameters is not None:
self.set_parameters(parameters)
if states is not None:
self.set_states(states)

self.lm_data = []
self.act_data = []

def set_parameters(self, var=None):
"""Load and set parameters for the muscle model.
"""

if var is None:
var = './../data/muscle_parameter.txt'
if isinstance(var, str):
self.P = self.config_parser(var, 'parameters')
elif isinstance(var, dict):
self.P = var
else:
raise ValueError('Wrong parameters!')

print('The parameters were successfully loaded ' +
'and are stored in the variable P.')

def set_states(self, var=None):
"""Load and set states for the muscle model.
"""

if var is None:
var = './../data/muscle_state.txt'
if isinstance(var, str):
self.S = self.config_parser(var, 'states')
elif isinstance(var, dict):
self.S = var
else:
raise ValueError('Wrong states!')

print('The states were successfully loaded ' +
'and are stored in the variable S.')

def config_parser(self, filename, var):

parser = configparser.ConfigParser()
parser.optionxform = str  # make option names case sensitive
if not parser:
#if not 'Muscle' in parser.sections()[0]:
#    raise ValueError('Wrong %s file!' %var)
var = {}
for key, value in parser.items(parser.sections()[0]):
if key.lower() in ['name', 'id']:
var.update({key: value})
else:
try:
value = float(value)
except ValueError:
print('"%s" value "%s" was replaced by NaN.' %(key, value))
value = np.nan
var.update({key: value})

return var

def force_l(self, lm, gammal=None):
"""Thelen (2003) force of the contractile element vs. muscle length.

Parameters
----------
lm : float
normalized muscle fiber length
gammal : float, optional (default from parameter file)
shape factor

Returns
-------
fl : float
normalized force of the muscle contractile element
"""

if gammal is None: gammal = self.P['gammal']

fl = np.exp(-(lm-1)**2/gammal)

return fl

def force_pe(self, lm, kpe=None, epsm0=None):
"""Thelen (2003) force of the muscle parallel element vs. muscle length.

Parameters
----------
lm : float
normalized muscle fiber length
kpe : float, optional (default from parameter file)
exponential shape factor
epsm0 : float, optional (default from parameter file)
passive muscle strain due to maximum isometric force

Returns
-------
fpe : float
normalized force of the muscle parallel (passive) element
"""

if kpe is None: kpe = self.P['kpe']
if epsm0 is None: epsm0 = self.P['epsm0']

if lm <= 1:
fpe = 0
else:
fpe = (np.exp(kpe*(lm-1)/epsm0)-1)/(np.exp(kpe)-1)

return fpe

def force_se(self, lt, ltslack=None, epst0=None, kttoe=None):
"""Thelen (2003) force-length relationship of tendon vs. tendon length.

Parameters
----------
lt : float
tendon length (normalized or not)
ltslack : float, optional (default from parameter file)
tendon slack length (normalized or not)
epst0 : float, optional (default from parameter file)
tendon strain at the maximal isometric muscle force
kttoe : float, optional (default from parameter file)
linear scale factor

Returns
-------
fse : float
normalized force of the tendon series element
"""

if ltslack is None: ltslack = self.P['ltslack']
if epst0 is None: epst0 = self.P['epst0']
if kttoe is None: kttoe = self.P['kttoe']

epst = (lt-ltslack)/ltslack
fttoe = .33
# values from OpenSim Thelen2003Muscle
epsttoe =  .99*epst0*np.e**3/(1.66*np.e**3 - .67)
ktlin =  .67/(epst0 - epsttoe)
#
if epst <= 0:
fse = 0
elif epst <= epsttoe:
fse = fttoe/(np.exp(kttoe)-1)*(np.exp(kttoe*epst/epsttoe)-1)
else:
fse = ktlin*(epst-epsttoe) + fttoe

return fse

def velo_fm(self, fm, a, fl, lmopt=None, vmmax=None, fmlen=None, af=None):
"""Thelen (2003) velocity of the force-velocity relationship vs. CE force.

Parameters
----------
fm : float
normalized muscle force
a : float
muscle activation level
fl : float
normalized muscle force due to the force-length relationship
lmopt : float, optional (default from parameter file)
optimal muscle fiber length
vmmax : float, optional (default from parameter file)
normalized maximum muscle velocity for concentric activation
fmlen : float, optional (default from parameter file)
normalized maximum force generated at the lengthening phase
af : float, optional (default from parameter file)
shape factor

Returns
-------
vm : float
velocity of the muscle
"""

if lmopt is None: lmopt = self.P['lmopt']
if vmmax is None: vmmax = self.P['vmmax']
if fmlen is None: fmlen = self.P['fmlen']
if af is None: af = self.P['af']

if fm <= a*fl:  # isometric and concentric activation
if fm > 0:
b = a*fl + fm/af
else:
b = a*fl
else:           # eccentric activation
asyE_thresh = 0.95  # from OpenSim Thelen2003Muscle
if fm < a*fl*fmlen*asyE_thresh:
b = (2 + 2/af)*(a*fl*fmlen - fm)/(fmlen - 1)
else:
fm0 = a*fl*fmlen*asyE_thresh
b = (2 + 2/af)*(a*fl*fmlen - fm0)/(fmlen - 1)

vm = (0.25  + 0.75*a)*1*(fm - a*fl)/b
vm = vm*vmmax*lmopt

return vm

def force_vm(self, vm, a, fl, lmopt=None, vmmax=None, fmlen=None, af=None):
"""Thelen (2003) force of the contractile element vs. muscle velocity.

Parameters
----------
vm : float
muscle velocity
a : float
muscle activation level
fl : float
normalized muscle force due to the force-length relationship
lmopt : float, optional (default from parameter file)
optimal muscle fiber length
vmmax : float, optional (default from parameter file)
normalized maximum muscle velocity for concentric activation
fmlen : float, optional (default from parameter file)
normalized normalized maximum force generated at the lengthening phase
af : float, optional (default from parameter file)
shape factor

Returns
-------
fvm : float
normalized force of the muscle contractile element
"""

if lmopt is None: lmopt = self.P['lmopt']
if vmmax is None: vmmax = self.P['vmmax']
if fmlen is None: fmlen = self.P['fmlen']
if af is None: af = self.P['af']

vmmax = vmmax*lmopt
if vm <= 0:  # isometric and concentric activation
fvm = af*a*fl*(4*vm + vmmax*(3*a + 1))/(-4*vm + vmmax*af*(3*a + 1))
else:        # eccentric activation
fvm = a*fl*(af*vmmax*(3*a*fmlen - 3*a + fmlen - 1) + \
8*vm*fmlen*(af + 1)) / \
(af*vmmax*(3*a*fmlen - 3*a + fmlen - 1) + 8*vm*(af + 1))

return fvm

def lmt_eq(self, t, lmt0=None):
"""Equation for muscle-tendon length."""

if lmt0 is None:
lmt0 = self.S['lmt0']

return lmt0

def vm_eq(self, t, lm, lm0, lmt0, lmopt, ltslack, alpha0, vmmax, fm0):
"""Equation for muscle velocity."""

if lm < 0.1*lmopt:
lm = 0.1*lmopt
#lt0 = lmt0 - lm0*np.cos(alpha0)
a = self.activation(t)
lmt = self.lmt_eq(t, lmt0)
alpha = self.penn_ang(lmt=lmt, lm=lm, lm0=lm0, alpha0=alpha0)
lt = lmt - lm*np.cos(alpha)
fse = self.force_se(lt=lt, ltslack=ltslack)
fpe = self.force_pe(lm=lm/lmopt)
fl = self.force_l(lm=lm/lmopt)
fce_t = fse/np.cos(alpha) - fpe
#if fce_t < 0: fce_t=0
vm = self.velo_fm(fm=fce_t, a=a, fl=fl)

return vm

def lm_sol(self, fun=None, t0=0, t1=3, lm0=None, lmt0=None, ltslack=None, lmopt=None,
alpha0=None, vmmax=None, fm0=None,  show=True, axs=None):
"""Runge-Kutta (4)5 ODE solver for muscle length."""

if lm0 is None: lm0 = self.S['lm0']
if lmt0 is None: lmt0 = self.S['lmt0']
if ltslack is None: ltslack = self.P['ltslack']
if alpha0 is None: alpha0 = self.P['alpha0']
if lmopt is None: lmopt = self.P['lmopt']
if vmmax is None: vmmax = self.P['vmmax']
if fm0 is None: fm0 = self.P['fm0']

if fun is None:
fun = self.vm_eq
f = ode(fun).set_integrator('dopri5', nsteps=1, max_step=0.005, atol=1e-8)
f.set_initial_value(lm0, t0).set_f_params(lm0, lmt0, lmopt, ltslack, alpha0, vmmax, fm0)
# suppress Fortran warning
warnings.filterwarnings("ignore", category=UserWarning)
data = []
while f.t < t1:
f.integrate(t1, step=True)
d = self.calc_data(f.t, np.max([f.y, 0.1*lmopt]), lm0, lmt0,
ltslack, lmopt, alpha0, fm0)
data.append(d)

warnings.resetwarnings()
data = np.array(data)
self.lm_data = data
if show:
self.lm_plot(data, axs)

return data

def calc_data(self, t, lm, lm0, lmt0, ltslack, lmopt, alpha0, fm0):
"""Calculus of muscle-tendon variables."""

a = self.activation(t)
lmt = self.lmt_eq(t, lmt0=lmt0)
alpha = self.penn_ang(lmt=lmt, lm=lm, lm0=lm0, alpha0=alpha0)
lt = lmt - lm*np.cos(alpha)
fl = self.force_l(lm=lm/lmopt)
fpe = self.force_pe(lm=lm/lmopt)
fse = self.force_se(lt=lt, ltslack=ltslack)
fce_t = fse/np.cos(alpha) - fpe
vm = self.velo_fm(fm=fce_t, a=a, fl=fl, lmopt=lmopt)
fm = self.force_vm(vm=vm, fl=fl, lmopt=lmopt, a=a) + fpe
data = [t, lmt, lm, lt, vm, fm*fm0, fse*fm0, a*fl*fm0, fpe*fm0, alpha]

return data

def muscle_plot(self, a=1, axs=None):
"""Plot muscle-tendon relationships with length and velocity."""

try:
import matplotlib.pyplot as plt
except ImportError:
print('matplotlib is not available.')
return

if axs is None:
_, axs = plt.subplots(nrows=1, ncols=3, figsize=(9, 4))

lmopt   = self.P['lmopt']
ltslack = self.P['ltslack']
vmmax   = self.P['vmmax']
alpha0  = self.P['alpha0']
fm0     = self.P['fm0']
lm0     = self.S['lm0']
lmt0    = self.S['lmt0']
lt0     = self.S['lt0']
if np.isnan(lt0):
lt0 = lmt0 - lm0*np.cos(alpha0)

lm  = np.linspace(0, 2, 101)
lt  = np.linspace(0, 1, 101)*0.05 + 1
vm  = np.linspace(-1, 1, 101)*vmmax*lmopt
fl  = np.zeros(lm.size)
fpe = np.zeros(lm.size)
fse = np.zeros(lt.size)
fvm = np.zeros(vm.size)

fl_lm0  = self.force_l(lm0/lmopt)
fpe_lm0 = self.force_pe(lm0/lmopt)
fm_lm0  = fl_lm0 + fpe_lm0
ft_lt0  = self.force_se(lt0, ltslack)*fm0

for i in range(101):
fl[i]  = self.force_l(lm[i])
fpe[i] = self.force_pe(lm[i])
fse[i] = self.force_se(lt[i], ltslack=1)
fvm[i] = self.force_vm(vm[i], a=a, fl=fl_lm0)

lm  = lm*lmopt
lt  = lt*ltslack
fl  = fl
fpe = fpe
fse = fse*fm0
fvm = fvm*fm0

xlim = self.margins(lm, margin=.05, minmargin=False)
axs[0].set_xlim(xlim)
ylim = self.margins([0, 2], margin=.05)
axs[0].set_ylim(ylim)
axs[0].plot(lm, fl, 'b', label='Active')
axs[0].plot(lm, fpe, 'b--', label='Passive')
axs[0].plot(lm, fl+fpe, 'b:', label='')
axs[0].plot([lm0, lm0], [ylim[0], fm_lm0], 'k:', lw=2, label='')
axs[0].plot([xlim[0], lm0], [fm_lm0, fm_lm0], 'k:', lw=2, label='')
axs[0].plot(lm0, fm_lm0, 'o', ms=6, mfc='r', mec='r', mew=2, label='fl(LM0)')
axs[0].legend(loc='best', frameon=True, framealpha=.5)
axs[0].set_xlabel('Length [m]')
axs[0].set_ylabel('Scale factor')
axs[0].xaxis.set_major_locator(plt.MaxNLocator(4))
axs[0].yaxis.set_major_locator(plt.MaxNLocator(4))
axs[0].set_title('Muscle F-L (a=1)')

xlim = self.margins([0, np.min(vm), np.max(vm)], margin=.05, minmargin=False)
axs[1].set_xlim(xlim)
ylim = self.margins([0, fm0*1.2, np.max(fvm)*1.5], margin=.025)
axs[1].set_ylim(ylim)
axs[1].plot(vm, fvm, label='')
axs[1].set_xlabel('$\mathbf{^{CON}}\;$ Velocity [m/s] $\;\mathbf{^{EXC}}$')
axs[1].plot([0, 0], [ylim[0], fvm[50]], 'k:', lw=2, label='')
axs[1].plot([xlim[0], 0], [fvm[50], fvm[50]], 'k:', lw=2, label='')
axs[1].plot(0, fvm[50], 'o', ms=6, mfc='r', mec='r', mew=2, label='FM0(LM0)')
axs[1].plot(xlim[0], fm0, '+', ms=10, mfc='r', mec='r', mew=2, label='')
axs[1].text(vm[0], fm0, 'FM0')
axs[1].legend(loc='upper right', frameon=True, framealpha=.5)
axs[1].set_ylabel('Force [N]')
axs[1].xaxis.set_major_locator(plt.MaxNLocator(4))
axs[1].yaxis.set_major_locator(plt.MaxNLocator(4))
axs[1].set_title('Muscle F-V (a=1)')

xlim = self.margins([lt0, ltslack, np.min(lt), np.max(lt)], margin=.05,
minmargin=False)
axs[2].set_xlim(xlim)
ylim = self.margins([ft_lt0, 0, np.max(fse)], margin=.05)
axs[2].set_ylim(ylim)
axs[2].plot(lt, fse, label='')
axs[2].set_xlabel('Length [m]')
axs[2].plot([lt0, lt0], [ylim[0], ft_lt0], 'k:', lw=2, label='')
axs[2].plot([xlim[0], lt0], [ft_lt0, ft_lt0], 'k:', lw=2, label='')
axs[2].plot(lt0, ft_lt0, 'o', ms=6, mfc='r', mec='r', mew=2, label='FT(LT0)')
axs[2].legend(loc='upper left', frameon=True, framealpha=.5)
axs[2].set_ylabel('Force [N]')
axs[2].xaxis.set_major_locator(plt.MaxNLocator(4))
axs[2].yaxis.set_major_locator(plt.MaxNLocator(4))
axs[2].set_title('Tendon')
plt.suptitle('Muscle-tendon mechanics', fontsize=18, y=1.03)
plt.show()

def lm_plot(self, x, axs):
"""Plot results of actdyn_ode45 function.
data = [t, lmt, lm, lt, vm, fm*fm0, fse*fm0, fl*fm0, fpe*fm0, alpha]
"""

try:
import matplotlib.pyplot as plt
except ImportError:
print('matplotlib is not available.')
return

if axs is None:
_, axs = plt.subplots(nrows=3, ncols=2, sharex=True, figsize=(10, 6))

axs[0, 0].plot(x[:, 0], x[:, 1], 'b', label='LMT')
lmt = x[:, 2]*np.cos(x[:, 9]) + x[:, 3]
if np.sum(x[:, 9]) > 0:
axs[0, 0].plot(x[:, 0], lmt, 'g--', label=r'$LM \cos \alpha + LT$')
else:
axs[0, 0].plot(x[:, 0], lmt, 'g--', label=r'LM+LT')
ylim = self.margins(x[:, 1], margin=.1)
axs[0, 0].set_ylim(ylim)
axs[0, 0].legend(framealpha=.5, loc='best')

axs[0, 1].plot(x[:, 0], x[:, 3], 'b')
#axs[0, 1].plot(x[:, 0], lt0*np.ones(len(x)), 'r')
ylim = self.margins(x[:, 3], margin=.1)
axs[0, 1].set_ylim(ylim)

axs[1, 0].plot(x[:, 0], x[:, 2], 'b')
#axs[1, 0].plot(x[:, 0], lmopt*np.ones(len(x)), 'r')
ylim = self.margins(x[:, 2], margin=.1)
axs[1, 0].set_ylim(ylim)

axs[1, 1].plot(x[:, 0], x[:, 4], 'b')
ylim = self.margins(x[:, 4], margin=.1)
axs[1, 1].set_ylim(ylim)

axs[2, 0].plot(x[:, 0], x[:, 5], 'b', label='Muscle')
axs[2, 0].plot(x[:, 0], x[:, 6], 'g--', label='Tendon')
ylim = self.margins(x[:, [5, 6]], margin=.1)
axs[2, 0].set_ylim(ylim)
axs[2, 0].set_xlabel('Time (s)')
axs[2, 0].legend(framealpha=.5, loc='best')

axs[2, 1].plot(x[:, 0], x[:, 8], 'b', label='PE')
ylim = self.margins(x[:, 8], margin=.1)
axs[2, 1].set_ylim(ylim)
axs[2, 1].set_xlabel('Time (s)')
axs[2, 1].legend(framealpha=.5, loc='best')

axs = axs.flatten()
ylabel = ['$L_{MT}\,(m)$', '$L_{T}\,(m)$', '$L_{M}\,(m)$',
'$V_{CE}\,(m/s)$', '$Force\,(N)$', '$Force\,(N)$']
for i, axi in enumerate(axs):
axi.set_ylabel(ylabel[i], fontsize=14)
axi.yaxis.set_major_locator(plt.MaxNLocator(4))
axi.yaxis.set_label_coords(-.2, 0.5)

plt.suptitle('Simulation of muscle-tendon mechanics', fontsize=18,
y=1.03)
plt.tight_layout()
plt.show()

def penn_ang(self, lmt, lm, lt=None, lm0=None, alpha0=None):
"""Pennation angle.

Parameters
----------
lmt : float
muscle-tendon length
lt : float, optional (default=None)
tendon length
lm : float, optional (default=None)
muscle fiber length
lm0 : float, optional (default from states file)
initial muscle fiber length
alpha0 : float, optional (default from parameter file)
initial pennation angle

Returns
-------
alpha : float
pennation angle
"""

if lm0 is None: lm0 = self.S['lm0']
if alpha0 is None: alpha0 = self.P['alpha0']

alpha = alpha0
if alpha0 != 0:
w = lm0*np.sin(alpha0)
if lm is not None:
cosalpha = np.sqrt(1-(w/lm)**2)
elif lmt is not None and lt is not None:
cosalpha = 1/(np.sqrt(1 + (w/(lmt-lt))**2))
alpha = np.arccos(cosalpha)

if alpha > 1.4706289:  # np.arccos(0.1), 84.2608 degrees
alpha = 1.4706289

return alpha

def excitation(self, t, u_max=None, u_min=None, t0=0, t1=5):
"""Excitation signal, a square wave.

Parameters
----------
t : float
time instant [s]
u_max : float (0 < u_max <= 1), optional (default from parameter file)
maximum value for muscle excitation
u_min : float (0 < u_min < 1), optional (default from parameter file)
minimum value for muscle excitation
t0 : float, optional (default=0)
initial time instant for muscle excitation equals to u_max [s]
t1 : float, optional (default=5)
final time instant for muscle excitation equals to u_max [s]

Returns
-------
u : float (0 < u <= 1)
excitation signal
"""

if u_max is None: u_max = self.P['u_max']
if u_min is None: u_min = self.P['u_min']

u = u_min
if t >= t0 and t <= t1:
u = u_max

return u

def activation_dyn(self, t, a, t_act=None, t_deact=None):
"""Thelen (2003) activation dynamics, the derivative of a at t.

Parameters
----------
t : float
time instant [s]
a : float (0 <= a <= 1)
muscle activation
t_act : float, optional (default from parameter file)
activation time constant [s]
t_deact : float, optional (default from parameter file)
deactivation time constant [s]

Returns
-------
derivative of a at t
"""

if t_act is None: t_act = self.P['t_act']
if t_deact is None: t_deact = self.P['t_deact']

u = self.excitation(t)
if u > a:
adot = (u - a)/(t_act*(0.5 + 1.5*a))
else:
adot = (u - a)/(t_deact/(0.5 + 1.5*a))

def activation_sol(self, fun=None, t0=0, t1=3, a0=0, u_min=None,
t_act=None, t_deact=None, show=True, axs=None):
"""Runge-Kutta (4)5 ODE solver for activation dynamics.

Parameters
----------
fun : function object, optional (default is None and actdyn is used)
function with ODE to be solved
t0 : float, optional (default=0)
initial time instant for the simulation [s]
t1 : float, optional (default=0)
final time instant for the simulation [s]
a0 : float, optional (default=0)
initial muscle activation
u_max : float (0 < u_max <= 1), optional (default from parameter file)
maximum value for muscle excitation
u_min : float (0 < u_min < 1), optional (default from parameter file)
minimum value for muscle excitation
t_act : float, optional (default from parameter file)
activation time constant [s]
t_deact : float, optional (default from parameter file)
deactivation time constant [s]
show : bool, optional (default = True)
if True (1), plot data in matplotlib figure
axs : a matplotlib.axes.Axes instance, optional (default = None)

Returns
-------
data : 2-d array
array with columns [time, excitation, activation]

"""

if u_min is None: u_min = self.P['u_min']
if t_act is None: t_act = self.P['t_act']
if t_deact is None: t_deact = self.P['t_deact']

if fun is None:
fun = self.activation_dyn
f = ode(fun).set_integrator('dopri5', nsteps=1, max_step=0.005, atol=1e-8)
f.set_initial_value(a0, t0).set_f_params(t_act, t_deact)
# suppress Fortran warning
warnings.filterwarnings("ignore", category=UserWarning)
data = []
while f.t < t1:
f.integrate(t1, step=True)
data.append([f.t, self.excitation(f.t), np.max([f.y, u_min])])
warnings.resetwarnings()
data = np.array(data)
if show:
self.actvation_plot(data, axs)

self.act_data = data

return data

def activation(self, t=None):
"""Activation signal."""

data = self.act_data
if t is not None and len(data):
if t <= self.act_data[0, 0]:
a = self.act_data[0, 2]
elif t >= self.act_data[-1, 0]:
a = self.act_data[-1, 2]
else:
a = np.interp(t, self.act_data[:, 0], self.act_data[:, 2])
else:
a = 1

return a

def actvation_plot(self, data, axs):
"""Plot results of actdyn_ode45 function."""

try:
import matplotlib.pyplot as plt
except ImportError:
print('matplotlib is not available.')
return

if axs is None:
_, axs = plt.subplots(nrows=1, ncols=1, figsize=(6, 4))

axs.plot(data[:, 0], data[:, 1], color=[1, 0, 0, .6], label='Excitation')
axs.plot(data[:, 0], data[:, 2], color=[0, 0, 1, .6], label='Activation')
axs.set_xlabel('Time [s]')
axs.set_ylabel('Level')
axs.legend()
plt.title('Activation dynamics')
plt.tight_layout()
plt.show()

def margins(self, x, margin=0.01, minmargin=True):
"""Calculate plot limits with extra margins.
"""
rang = np.nanmax(x) - np.nanmin(x)
if rang < 0.001 and minmargin:
rang = 0.001*np.nanmean(x)/margin
if rang < 1:
rang = 1
lim = [np.nanmin(x) - rang*margin, np.nanmax(x) + rang*margin]

return lim