# Lecture 9 - Motor Control¶

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

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

• In class:
In [115]:
import numpy as np
#import pandas as pd
#import pylab as pl
import matplotlib.pyplot as plt
import math

%matplotlib notebook


### Muscle properties¶

In [116]:
Lslack = .223
Lce_o = .093 #optmal l

Fmax = 3000


### Initial conditions¶

In [117]:
LceNorm = .087/Lce_o
Lm = .31
act = 0
t0 = 0
tf = 2.99
h = 1e-3

In [118]:
t = np.arange(t0,tf,h)
F = np.empty(t.shape)
Fkpe = np.empty(t.shape)
fiberLength = np.empty(t.shape)
tendonLength = np.empty(t.shape)
activation = np.zeros(t.shape)
neural = np.zeros(t.shape)

fiberLength[0] = LceNorm * Lce_o
tendonLength[0] = Lm - fiberLength[0]

In [119]:
def computeTendonForce(LseeNorm, Lslack, Lce_o):
'''
Compute Tendon Force

Imputs:

LseeNorm - Normalized Tendon length

Lslack - slack length of the tendon (not-normalized)

Lce_o - Optimal length of the fiber

Output:

FTendonNorm - Normalized tendon force
'''

Umax = .04
if LseeNorm<Lslack/Lce_o:
FTendonNorm = 0
else:
FTendonNorm = ((LseeNorm-Lslack/Lce_o)/(Umax*Lslack/Lce_o))**2

return FTendonNorm

In [120]:
def computeParallelElementForce(LceNorm):
Umax = 1

if LceNorm<1:
FkpeNorm = 0
else:
FkpeNorm = ((LceNorm-1)/Umax)**2

return FkpeNorm

In [121]:
def computeForceLengthCurve(LceNorm):

width = .63
F0 = max([0, (1-((LceNorm-1)/width)**2)])

return F0

In [122]:
def computeContractilElementDerivative (act, F0, FCE):

A = 0.25
Fm_len = 1.4
v_max = 10

if  FCE <= act*F0:
b = act*F0 + FCE/A

else:
b = (2 + 2/A)*(act*F0*Fm_len + FCE)/(Fm_len - 1)

LceNormdot = (0.25+0.75*act)*v_max*(FCE-act*F0)/b

return LceNormdot

In [123]:
def computeContractilElementForce (FTendonNorm, FkpeNorm):
FCE = FTendonNorm - FkpeNorm
return FCE

In [124]:
def computeTendonLength(Lm, Lce_o, LceNorm):
LseeNorm = Lm/Lce_o - LceNorm
return LseeNorm

In [125]:
def  computeActivation (act, n, dt):

tau = 0.12

dactdt = (n-act)/tau
act = act + dt*dactdt

return act


## Simulation - Parallel¶

In [126]:
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)

if t[i]<=1 and np.random.rand(1)< 0.3:
n = 0.5
elif t[i]>1.5 and t[i]<2.5 and np.random.rand(1)< 0.5:
n = 0.5
elif t[i]>=2.5 and np.random.rand(1)< 0.01:
n = 0.5
else:
n = 0.5

##################################################
LseeNorm = computeTendonLength(Lm, Lce_o, LceNorm)

FTendonNorm = computeTendonForce(LseeNorm, Lslack, Lce_o)

FkpeNorm = computeParallelElementForce(LceNorm)

F0 = computeForceLengthCurve (LceNorm)

FCE = computeContractilElementForce (FTendonNorm, FkpeNorm)

act = computeActivation(act, n, h)

LceNormdot = computeContractilElementForce (act, F0, FCE)

# --- Euler integration step
LceNorm = LceNorm + h*LceNormdot

F[i] = FTendonNorm * Fmax
Fkpe[i] = FkpeNorm * Fmax
fiberLength[i] = LceNorm * Lce_o
tendonLength[i] = LseeNorm * Lce_o
neural[i] = n


---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-126-8b186d5c87f6> in <module>()
29     act = computeActivation(act, n, h)
30
---> 31     LceNormdot = computeContractilElementForce (act, F0, FCE)
32
33

TypeError: computeContractilElementForce() takes 2 positional arguments but 3 were given

## Plot¶

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

ax.plot(t,fiberLength, label = 'Fiber')
ax.plot(t,tendonLength, label = 'Tendon')
ax.plot(t,fiberLength+tendonLength, label = 'Fiber+Tendon')
plt.grid()
plt.legend(loc = 'best')
plt.xlabel('time (s)')
plt.ylabel('Length (m)')
plt. tight_layout()

In [ ]: