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_deact = 50e-3 #young adults
    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)
    
    #-------
    dadt = (u-a)/tau_a # euler
    
    a = a + dadt*dt
    #-------
    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 [ ]: