This notebook presents a simulation of 5000 ms of 400 descending commands and 800 motoneurons from soleus. The force is prduced by a Hill-type muscle model.


In [11]:
import sys
sys.path.insert(0, '..')
import time
import matplotlib.pyplot as plt
%matplotlib notebook 
from IPython.display import set_matplotlib_formats
set_matplotlib_formats('pdf', 'png')
plt.rcParams['savefig.dpi'] = 75

plt.rcParams['figure.autolayout'] = False
plt.rcParams['figure.figsize'] = 6, 6
plt.rcParams['axes.labelsize'] = 18
plt.rcParams['axes.titlesize'] = 20
plt.rcParams['font.size'] = 16
plt.rcParams['lines.linewidth'] = 2.0
plt.rcParams['lines.markersize'] = 8
plt.rcParams['legend.fontsize'] = 14

plt.rcParams['text.usetex'] = True
plt.rcParams['font.family'] = "serif"
plt.rcParams['font.serif'] = "cm"
plt.rcParams['text.latex.preamble'] = "\usepackage{subdepth}, \usepackage{type1cm}"


import numpy as np

from Configuration import Configuration
from MotorUnitPool import MotorUnitPool
from NeuralTract import NeuralTract
from AfferentPool import AfferentPool
from SynapsesFactory import SynapsesFactory
from jointAnkleForceTask import jointAnkleForceTask

import scipy as sc
from scipy.signal import *

In [12]:
conf = Configuration('confMUProperties.rmto')
conf.simDuration_ms = 700 # Here I change simulation duration without changing the Configuration file.
t = np.arange(0.0, conf.simDuration_ms, conf.timeStep_ms)

In [13]:
pools = dict()
pools[0] = MotorUnitPool(conf, 'TA')


Muscle spindle from muscle TA built.
Motor Unit Pool TA built

In [14]:
aS1 = np.zeros_like(t)
aS3 = np.zeros_like(t)
aFR1 = np.zeros_like(t)
aFR3 = np.zeros_like(t)
aFF1 = np.zeros_like(t)
aFF3 = np.zeros_like(t)
tic = time.time()
for i in xrange(0, len(t)-1):
    for j in xrange(len(pools[0].unit)):
        if t[i] <= 25:
            pools[0].iInjected[2*j+1] = 10+8*j
        else:
            pools[0].iInjected[2*j+1] = 0
    pools[0].atualizeMotorUnitPool(t[i])
    aS1[i] = pools[0].Activation.activation_Sat[0] 
    aS3[i] = pools[0].Activation.activation_Sat[2] 
    aFR1[i] = pools[0].Activation.activation_Sat[3] 
    aFR3[i] = pools[0].Activation.activation_Sat[5]
    aFF1[i] = pools[0].Activation.activation_Sat[6] 
    aFF3[i] = pools[0].Activation.activation_Sat[8] 
toc = time.time()
print str(toc - tic) + ' seconds'


292.335968018 seconds

In [ ]:
tc = np.array([])
for i in xrange(len(pools[0].unit)):
    tc = np.append(tc, pools[0].unit[i].TwitchTc_ms)

In [ ]:
plt.figure()
plt.hist(tc,14)
plt.show()
print len(pools[0].unit)

tc1 = np.empty_like(tc) P0 = tc[0] Pn = tc[-1] n = len(tc) for i in xrange(250): tc1[i] = 110 np.exp(np.log(86.5/110)/250i)

for i in xrange(250,300): tc1[i] = 86.5 np.exp(np.log(55.25/86.5)/50(i-250))

for i in xrange(300,350): tc1[i] = 55.25 np.exp(np.log(25/55.25)/50(i-300))

plt.figure() plt.hist(tc1,14) plt.show()

tc2 = np.empty_like(tc) P0 = tc[0] Pn = tc[-1] n = len(tc) for i in xrange(n): tc2[i] = P0 np.exp(np.log(Pn/P0)/ni)

plt.figure() plt.hist(tc2,14) plt.show()

n = np.linspace(0,349,350) tc3 = 110.0 np.exp(np.log(25/110.0)/350.0n)

plt.figure() plt.hist(tc3,14) plt.show()

tc1 = np.empty_like(tc3)

for i in xrange(250): tc1[i] = 110 np.exp(np.log(86.5/110)/250i)

for i in xrange(250,300): tc1[i] = 86.5 np.exp(np.log(55.25/86.5)/50(i-250))

for i in xrange(300,350): tc1[i] = 55.25 np.exp(np.log(25/55.25)/50(i-300))

plt.figure() plt.plot(n,tc1,'r-') plt.xlim((0,350)) plt.ylim((20,120)) plt.show()

tc4 = np.empty_like(tc1)

for i in xrange(250): tc4[i] = (110-40.0) np.exp(-0.02i)+40.0

for i in xrange(250,350): tc4[i]= (40 - 20.25) (1 - np.exp(1/50.0np.log((35-40.0)/(20.25-40))*(350-(i)))) + 20.25

for i in xrange(300,350):

tc4[i] = 55.25 - 30.25(1 - np.exp(-30.25/500(i-300)))

plt.figure() plt.plot(n,tc4,'r-') plt.xlim((0,350)) plt.ylim((20,120)) plt.show()

plt.figure() plt.hist(tc4,14) plt.show()

plt.plot(n,(250-50)np.exp(-0.01n)+(50)np.exp(-0.01n))

tc5 = np.empty_like(tc1) ''' VS1 = 100.0 VFR1 = 45.0 VFF1 = 35.0 VFFN = 20.25 ''' VS1 = 20.25 VFR1 = 35.0 VFF1 = 45.0 VFFN = 100 Nmu = 350 NmuS = 250 NmuFF = 50

for i in xrange(350): tc5[i] = ((VS1 - VFR1) np.exp(-5.0i/NmuS) + (VFR1 - VFFN) (1.0 - np.exp(1.0/NmuFFnp.log((VFF1 - VFR1)/(VFFN - VFR1)) * (Nmu - i))) + VFFN)

plt.figure() plt.plot(n,tc5,'r-') plt.xlim((0,350)) plt.ylim((20,120)) plt.show()

plt.figure() plt.hist(tc5,14,rwidth=0.6,color='grey') plt.show()


In [ ]:


In [7]:
plt.figure()
plt.plot(t,aS1)
plt.show()



In [8]:
plt.figure()
plt.plot(t,aS3)
plt.show()



In [9]:
2/(1+np.exp(-0.7))-1


Out[9]:
0.33637554433633232

In [10]:
1/2.97


Out[10]:
0.33670033670033667

In [11]:
plt.figure()
plt.plot(t, pools[0].Muscle.force)


Out[11]:
[<matplotlib.lines.Line2D at 0x7f3c4e3f1fd0>]

In [ ]: