# Task 14 - Motor Control¶

### Introduction to modeling and simulation of human movement¶

https://github.com/BMClab/bmc/blob/master/courses/ModSim2018.md

Implement a simulation of the ankle joint model using the parameters from Thelen (2003) and Elias (2014)

In [1]:
import numpy as np
import pandas as pd
%matplotlib notebook
import matplotlib.pyplot as plt
import math


### Muscle properties¶

Parameters obtained from Table 2 of Thelen (2003).

In [2]:
Lslack = 2.4*0.09 # tendon slack length
Lce_o  = 0.09 # optimal muscle fiber length
Fmax   = 1400 #maximal isometric DF force
alpha  = 7*math.pi/180 # DF muscle fiber pennation angle


Parameters from Nigg & Herzog (2006).

In [3]:
Umax  = 0.04 # SEE strain at Fmax
width = 0.63 # Max relative length change of CE


Activation dynamics parameters

In [4]:
a = 1 #inital conditional for ativation
u = 1 #Initial conditional for Brain's activation
#b = .25*10#*Lce_o


## Subject's anthropometrics¶

Parameters obtained experimentally or from Winter's book.

In [5]:
M      = 75 #total body mass (kg)
Lseg   = 0.26 #segment length (m)
m      = 0.0145*M #foot mass (kg)
g      = 9.81 #acceleration of gravity (m/s2)
Rcm    = Lseg*0.5 #distance from ankle joint to center of mass (m)
I      = m*(0.69*Lseg)**2 #moment of inertia
legAng = math.pi/2 #angle of the leg with horizontal (90 deg)


### Initial conditions¶

In [6]:
phi  = 0 #start as neutral ankle angle (0 degrees)
phid = 0 #zero velocity
Lm0  = 0.31 #initial total lenght of the muscle
Lnorm_ce = .087/Lce_o #norm
t0 = 0 #Initial time
tf = 0.5 #Final Time
h  = 1e-4 #integration step size and step counter

In [7]:
t = np.arange(t0,tf,h) # time array
# preallocating
F          = np.empty(t.shape)
phivec     = np.empty(t.shape)
Fkpe       = np.empty(t.shape)
FiberLen   = np.empty(t.shape)
TendonLen  = np.empty(t.shape)
a_dynamics = np.empty(t.shape)
Moment     = np.empty(t.shape)


## Simulation - Series¶

In [8]:
def momentArmDF(phi):
'''
Calculate the tibialis anterior moment arm according to Elias et al (2014)
Input:
phi: Ankle joint angle in radians
Output:
Rarm: TA moment arm
'''
# Consider neutral ankle position as zero degrees
phi = phi*180/np.pi # converting to degrees
Rf = 4.3 + 1.66E-2*phi + -3.89E-4*phi**2 + -4.45E-6*phi**3 + -4.34E-8*phi**4
Rf = Rf/100 # converting to meters
return Rf

In [9]:
def ComputeTotalLengthSizeTA(phi):
'''
Calculate TA MTU length size according to Elias et al (2014)
Input:
phi: ankle angle
'''
phi = phi*180/math.pi # converting to degrees
Lm = 30.6 + -7.44E-2*phi + -1.41E-4*phi**2 + 2.42E-6*phi**3 + 1.5E-8*phi**4
Lm = Lm/100
return Lm

In [10]:
def tendonLength(Lm, Lce_o, Lnorm_ce, alpha):
'''
Compute tendon length

Inputs:
Lm =
Lce_o = optimal length of the fiber
Lnorm_ce = normalized contractile element length

Output:
Lnorm_see = normalized tendon length
'''
Lnorm_see = Lm/Lce_o - Lnorm_ce*np.cos(alpha)

return Lnorm_see

In [11]:
def TendonForce(Lnorm_see, Lslack, Lce_o):
'''
Compute tendon force

Inputs:
Lnorm_see = normalized tendon length
Lslack = slack length of the tendon (non-normalized)
Lce_o = optimal length of the fiber

Output:
Fnorm_tendon = normalized tendon force

'''
Umax = 0.04

if Lnorm_see<Lslack/Lce_o:
Fnorm_tendon = 0
else:
Fnorm_tendon = ((Lnorm_see-Lslack/Lce_o)/(Umax*Lslack/Lce_o))**2

return Fnorm_tendon

In [12]:
def ParallelElementForce(Lnorm_ce):
'''
Compute parallel element force

Inputs:
Lnorm_ce = normalized contractile element length

Output:
Fnorm_kpe = normalized parallel element force

'''
Umax = 1

if Lnorm_ce< 1:
Fnorm_kpe = 0
else:
Fnorm_kpe = ((Lnorm_ce-1)/(Umax*1))**2

return Fnorm_kpe

In [13]:
def ForceLengthCurve(Lnorm_ce, width):
F0 = max([0, (1-((Lnorm_ce-1)/width)**2)])
return F0

In [14]:
def ContractileElementDot(F0, Fnorm_CE, a):

'''
Compute Contractile Element Derivative

Inputs:
F0 = Force-Length Curve
Fce = Contractile element force

Output:
Lnorm_cedot = normalized contractile element length derivative

'''

FMlen = 1.4 # young adults
Vmax = 10  # young adults
Af = 0.25  #force-velocity shape factor

Fnorm_CE = min(FMlen*a*F0 - 0.001, Fnorm_CE)

if Fnorm_CE > a*F0:

b = ((2 + 2/Af)*(a*F0*FMlen - Fnorm_CE))/(FMlen-1)

elif Fnorm_CE <= a*F0:

b = a*F0 + Fnorm_CE/Af

Lnorm_cedot = (.25 + .75*a)*Vmax*((Fnorm_CE - a*F0)/b)

return Lnorm_cedot

In [15]:
def ContractileElementForce(Fnorm_tendon, Fnorm_kpe, alpha):
'''
Compute Contractile Element force

Inputs:
Fnorm_tendon = normalized tendon force
Fnorm_kpe = normalized parallel element force

Output:
Fnorm_CE = normalized contractile element force
'''
Fnorm_CE = Fnorm_tendon/np.cos(alpha) - Fnorm_kpe
return Fnorm_CE

In [16]:
def activation(a, u, dt):
'''
Compute activation

Inputs:
u = idealized muscle excitation signal, 0 <= u <= 1
a = muscular activation
dt = time step

Output:
a = muscular activation
'''

tau_act = 15e-3

if u>a:
tau_a = tau_act*(0.5+1.5*a)
elif u <=a:
tau_a = tau_deact/(0.5+1.5*a)

#-------

#-------
return a

In [17]:
def ComputeTotalLenghtSize(Lm0, phi, Rf, Rcm):
'''
Inputs:
Lm0 = initial lenght of the muscle
Phi = degree flexion of the joint
RF = Moment arm
Lce_o = optimal size of the muscle
Output:
Lm = total muscle lenght
'''
Lm = Lm0 - (phi-(math.pi/2))*Rf #total muscle-tendon length from joint angle
return Lm

In [18]:
def ComputeMomentJoint(Rf, Fnorm_tendon, Fmax, m, g, phi):
'''
Inputs:
RF = Moment arm
Fnorm_tendon = Normalized tendon force
m = Segment Mass
g = Acelleration of gravity
Fmax= maximal isometric force
Output:
M = Total moment with respect to joint
'''
M = Rf*Fnorm_tendon*Fmax - m*g*Rcm*np.sin(math.pi/2 - phi)
return M

In [19]:
def ComputeAngularAcelerationJoint(M, I):
'''
Inputs:
M = Total moment with respect to joint
I = Moment of Inertia
Output:
phidd= angular aceleration of the joint
'''
phidd = M/I
return phidd


## Simulation - Parallel¶

Checar os parâmetros se estão ok aqui!

In [20]:
#Normalizing

for i in range (len(t)):
#ramp
#if t[i]<=1:
#Lm = 0.31
#elif t[i]>1 and t[i]<2:
#Lm = .31 - .04*(t[i]-1)
#print(Lm)

#shortening at 4cm/s
#u = 0.7 + 0.2*np.sin(np.pi*t[i])

Lm = ComputeTotalLengthSizeTA(phi)

Rf = momentArmDF(phi)

Lnorm_see = tendonLength(Lm,Lce_o,Lnorm_ce,alpha)

Fnorm_tendon = TendonForce(Lnorm_see,Lslack, Lce_o)

Fnorm_kpe = ParallelElementForce (Lnorm_ce)

#isometric force at Lce from CE force length relationship
F0 = ForceLengthCurve (Lnorm_ce,width)

Fnorm_CE = ContractileElementForce(Fnorm_tendon,Fnorm_kpe, alpha) #Fnorm_CE = ~Fm

#computing activation
a = activation(a,u,h)

#calculate CE velocity from Hill's equation
Lnorm_cedot = ContractileElementDot(F0, Fnorm_CE,a)

#Compute MomentJoint
M = ComputeMomentJoint(Rf,Fnorm_tendon,Fmax,m,g,phi)

#Compute Angular Aceleration Joint
phidd = ComputeAngularAcelerationJoint (M,I)

# Euler integration steps
Lnorm_ce = Lnorm_ce + h*Lnorm_cedot
phid= phid + h*phidd
phi  = phi  + h*phid
phideg= (phi*180)/math.pi #convert joint angle from radians to degree
# Store variables in vectors
F[i] = Fnorm_tendon*Fmax
Fkpe[i] = Fnorm_kpe*Fmax
FiberLen[i] = Lnorm_ce*Lce_o
TendonLen[i] = Lnorm_see*Lce_o
a_dynamics[i] = a
phivec[i] = phideg
Moment[i] = M


## Plots¶

In [26]:
fig, ax = plt.subplots(1, 1, figsize=(6,4))

ax.plot(t,a_dynamics,c='magenta')
plt.grid()
plt.xlabel('time (s)')
plt.ylabel('Activation dynamics')
plt.show()

In [27]:
fig, ax = plt.subplots(1, 1, figsize=(6,4))

ax.plot(t, Moment)
plt.grid()
plt.xlabel('time (s)')
plt.ylabel('joint moment')
plt.show()

In [28]:
fig, ax = plt.subplots(1, 1, figsize=(6,4))

ax.plot(t, F, c='red')
plt.grid()
plt.xlabel('time (s)')
plt.ylabel('Force (N)')
plt.show()

In [29]:
fig, ax = plt.subplots(1, 1, figsize=(6,4))

ax.plot(t,phivec,c='red')
plt.grid()
plt.xlabel('time (s)')
plt.ylabel('Joint angle (deg)')
plt.show()

In [30]:
fig, ax = plt.subplots(1, 1, figsize=(6,4))

ax.plot(t,FiberLen, label = 'fiber')
ax.plot(t,TendonLen, label = 'tendon')
plt.grid()
plt.xlabel('time (s)')
plt.ylabel('Length (m)')
ax.legend(loc='best')

fig, ax = plt.subplots(1, 3, figsize=(9,4), sharex=True, sharey=True)
ax[0].plot(t,FiberLen, label = 'fiber')
ax[1].plot(t,TendonLen, label = 'tendon')
ax[2].plot(t,FiberLen + TendonLen, label = 'muscle (tendon + fiber)')

ax[1].set_xlabel('time (s)')
ax[0].set_ylabel('Length (m)')
ax[0].legend(loc='best')
ax[1].legend(loc='best')
ax[2].legend(loc='best')
plt.show()

In [ ]: