Code: Select all
import os
from neuron import h
import matplotlib.pyplot as plt
import numpy as np
import multiprocessing
import time
from scipy import interpolate
import h5py
#------------------------------------------------------------------------
#Load Motor Neuron template with incorporated mechanisms. (Kim 2017)
h.load_file('kim.hoc')
#Load standard NEURON run controls.
h.load_file("stdrun.hoc")
#------------------------------------------------------------------------
#Function Definitions----------------------------------------------------
def heav(x):#Needed function to calculate current of stimulation source.
if x<0:
return 0
else:
return 1
def Soma_Diam_Change(NMU):#Determine Motor Pool Motor Neuron size. (Volk 2021)
num_MU = NMU
RR = 50.9
a =np.log(RR)/num_MU
RTE=np.empty(num_MU)
for j in range(1,num_MU+1):
RTE[j-1]=np.exp(a*j)
diam_ori =48.8
diameters = RTE + diam_ori
return diameters#Return list of soma sizes in order of Motor Units.
def MU_List(N):#Generate Python Object from the Motor Neuron template.
MU_list = []
MU=h.Kim()
MU_list.append(MU)
return MU_list#Returns a list with a Motor Unit Py Object
def Soma_List(MU_list):#Isolate the Soma Segment of each Motor Unit Object
soma_list=[]#Initialize list
dend_list=[]
for section in h.allsec(MU_list):#For each Soma section on each Motor Neuron Object add them to a list
if h.issection('.*so.*', sec=section):
soma_list.append(section)
if h.issection('.*dend.*', sec=section):
dend_list.append(section)
return soma_list, dend_list#Return list of Somas
def Change_Dia(mu_list, soma_list,dend_list,N,NMU):#Change Diameters of the Somas of each Motor Neuron Object
NMU=NMU #Number of Motor Neurons desired
diameters = Soma_Diam_Change(NMU) #Calculate Soma Diameters
for i in range(len(mu_list)):#For each Motor Neuron Object in list
for j in range(soma_list[i].n3d()):#For each section of each Soma of each Motor Neuron Object
if j == 0:
mu_list[i].soma.pt3dchange(j, 0, 0, 0, diameters[N])
else:
mu_list[i].soma.pt3dchange(j, diameters[N], 0, 0, diameters[N])
for i in range(len(mu_list)):
mu_list[i].soma.gnafbar_Naf = ((mu_list[i].soma.gnafbar_Naf)-np.random.normal(0.0,0.1))
# mu_list[i].soma.ena = (mu_list[i].soma.ena+np.random.normal(0,5))
# mu_list[i].soma.g_pas = (mu_list[i].soma.g_pas+np.random.normal(0,1e-6))
# mu_list[i].soma.ek= (mu_list[i].soma.ek+np.random.normal(0,5))
mu_list[i].hillock.gnafbar_Naf = ((mu_list[i].hillock.gnafbar_Naf*np.random.normal(0.5,0.1)))
mu_list[i].hillock.gkdrbar_KDr = ((mu_list[i].hillock.gkdrbar_KDr+np.random.normal(0,0.1)))
# for j in range(len(dend_list)):
# mu_list[i].dend[j].g_pas =mu_list[i].dend[j].g_pas+np.random.normal(0,1e-6)
# mu_list[i].dend[j].e_pas =mu_list[i].dend[j].e_pas+np.random.normal(0,5)
# for i in range(len(mu_list)):
# for j in range(len(mu_list[i].iCaL)):
# if mu_list[i].iCaL[j] is None:
# ww=0
# else:
# mu_list[i].iCaL[j].gcalbar =mu_list[i].iCaL[j].gcalbar + np.random.normal(0,0.0008)*i
return mu_list#Return list of updated Motor Neuron Objects
def Record(stimulus, MU_list):#Record simulation characteristics of interest
#Initialize lists
t = []
v = []
amp = []
mgi = []
ina = []
# h_naf = []
for i in range(len(MU_list)):#For each Motor Neuron Object record Time, Soma Voltage, Stimulus Current, and Muscle Unit force
t.append(h.Vector().record(h._ref_t))
v.append(h.Vector().record(MU_list[i].soma(0.5)._ref_v))
amp.append(h.Vector().record(stimulus[i]._ref_amp))
mgi.append(h.Vector().record(MU_list[i].muscle_unit( 0.5 )._ref_mgi))
cell_t, cell_id = Spike_Detector(MU_list)
ina.append(h.Vector().record(MU_list[i].soma(0.5)._ref_ina))
return t, v, amp, mgi, cell_t, cell_id, ina #Return each variable recorded above as a NEURON Vector
def Spike_Detector(mu_list):
cell_t = h.Vector()
cell_id = h.Vector()
for i in range(len(mu_list)):
nc = h.NetCon(mu_list[i].soma(0.5)._ref_v,None, sec=mu_list[i].soma)
nc.record(cell_t, cell_id)
return cell_t, cell_id
def Run_Sim(N):#Simulation parameters
mu_list = MU_List(N)#Make Motor Neuron object list
soma_list, dend_list = Soma_List(mu_list)#Make Soma list
MU_Fin = Change_Dia(mu_list,soma_list,dend_list,N,200)#Updated Motor Neuron object list
#Print Coordinates of cells (Used to check results)
# print(MU_Fin[0].soma.x3d(0),MU_Fin[0].soma.y3d(1),MU_Fin[0].soma.z3d(0),MU_Fin[0].soma.diam3d(0))
# print(MU_Fin[1].soma.x3d(0),MU_Fin[1].soma.y3d(0),MU_Fin[1].soma.z3d(0),MU_Fin[1].soma.diam3d(0))
# print(MU_Fin[2].soma.x3d(0),MU_Fin[2].soma.y3d(0),MU_Fin[2].soma.z3d(0),MU_Fin[2].soma.diam3d(0))
# print(MU_Fin[0].soma.x3d(1),MU_Fin[0].soma.y3d(1),MU_Fin[0].soma.z3d(1),MU_Fin[0].soma.diam3d(1))
# print(MU_Fin[1].soma.x3d(1),MU_Fin[1].soma.y3d(1),MU_Fin[1].soma.z3d(1),MU_Fin[1].soma.diam3d(1))
# print(MU_Fin[2].soma.x3d(1),MU_Fin[2].soma.y3d(1),MU_Fin[2].soma.z3d(1),MU_Fin[2].soma.diam3d(1))
mul=40
WU = 5000
WUf = 5000
time1=50000 + WU+ WUf#Total time of simulation (ms)
dela=0
dur=32500
dura = (dur+WU)*mul
ramp_dur = (17500+WU)*mul
WUd = WU*mul
WUfd = WUf*mul
pkamp=18.1
bias=0
min = 6
std=1
#(nA)
on = 0*mul #ms
#---------------------------------------------------
fr_mean = np.genfromtxt('Trap_signal_T3.txt', delimiter=",")
fr_mean = np.interp(fr_mean, (fr_mean.min(), fr_mean.max()), (min, pkamp))
fr_array = np.ones([len(fr_mean),2])
for i in range(len(fr_mean)):
fr_array[i,:] = [i,fr_mean[i]]
desired_time = 5
time1= int(desired_time)
iv=np.arange(0,time1,0.001) #Time vector for current calculation.
#Make Neuron Vector.
istim_cur = np.empty(len(iv))
ival=0
time_fact = desired_time/np.max(fr_array[:,0])
fr_array[:,0]=fr_array[:,0]*time_fact
f = interpolate.interp1d(fr_array[:,0],fr_array[:,1])
istim_cur = f(iv)
frequency = 1
amplitude=1
sampling_rate= len(istim_cur)/time1
duration=time1
t=np.arange(0, duration, 1/sampling_rate)
sine_wave= amplitude * np.sin(2*np.pi*frequency*t)
istim_cur += sine_wave
for i in range(len(istim_cur)):
istim_cur[i] += np.random.normal(0,1)
stim = [] #Initialize list
iv=iv*1000
ivv = h.Vector().from_python(iv)
for i in range(len(MU_Fin)): #Generate and apply Stimulation object to each Soma of each Motor Neuron object.
stim.append(h.IClamp(MU_Fin[i].soma(0.5)))
istim_cur_1 = istim_cur
for i in range(len(stim)): #Generate stimulation parameters for each Stimulation object.
istim_cur_1 = istim_cur
for j in range(len(istim_cur)):
istim_cur_1[j] = istim_cur[j] + np.random.normal(0,1.5)
i_tv=h.Vector().from_python(istim_cur_1) #Generate Neuron Vector.
stim[i].dur = time1*1000
stim[i].delay = dela
i_tv.play(stim[i]._ref_amp, ivv)
stimulus = stim
t, v, amp, mgi, cell_t, cell_id, ina = Record(stimulus, MU_Fin) #Record desired characteristics.
#Initialize Simulation.
h.finitialize(-70)
h.celsius = 36
h.dt = 0.01
h.continuerun(time1*1000)
return t, v, amp, mgi, cell_t, cell_id, ina #Return desired recrodings.
#------------------------------------------------------------------------
if __name__ == "__main__": #Semantic statement for Multiprocessing.
hdf5_file_path = 'simulation_data.h5'
# Delete the existing HDF5 file if it exists
if os.path.exists(hdf5_file_path):
os.remove(hdf5_file_path)
print(f"Deleted existing {hdf5_file_path}")
print("I Started")
num=200 #Num of desired Neurons
N = [i for i in range(num)] #Generate range for Multiprocessing command.
MU_Model = [] #Model List.
start = time.perf_counter() #For elapsed time counter
with multiprocessing.Pool(6) as pool: #Run simulation with as many available cpus as possible.
MU_Model.append(pool.map(Run_Sim, N))
print('Finished Sim ready to write')
#MU_model[i][j][k][l]
# i: access the list
# j: access desired Motor Neuron
# k: access desired recorded characteristic (Time, Voltage, Current, Force, mgi)
# l: acces desired characteristic vector
holder = np.ones(len(MU_Model[0]))
for i in range(len(MU_Model[0])):
holder[i] = len(MU_Model[0][i][4])
mval = np.argmax(holder)
t_array = np.empty([len(MU_Model[0][0][0][0]),len(MU_Model[0])])
v_array = np.empty([len(MU_Model[0][0][0][0]),len(MU_Model[0])])
amp_array = np.empty([len(MU_Model[0][0][0][0]),len(MU_Model[0])])
mgi_array = np.empty([len(MU_Model[0][0][3][0]),len(MU_Model[0])])
cell_t_array = np.ones([len(MU_Model[0][mval][4]),len(MU_Model[0])])
cell_id_array = np.ones([len(MU_Model[0][0][5]),len(MU_Model[0])])
ina_array = np.ones([len(MU_Model[0][0][6][0]),len(MU_Model[0])])
# Create an HDF5 file to store your data
with h5py.File('simulation_data.h5', 'w') as file:
# Create datasets for each type of data
for i in range(len(MU_Model[0])):
file.create_dataset(f't_array_{i}', data=np.array(MU_Model[0][i][0][0]))
file.create_dataset(f'v_array_{i}', data=np.array(MU_Model[0][i][1][0]))
file.create_dataset(f'amp_array_{i}', data=np.array(MU_Model[0][i][2][0]))
file.create_dataset(f'mgi_array_{i}', data=np.array(MU_Model[0][i][3][:]))
if len(MU_Model[0][i][4])>0:
file.create_dataset(f'cell_t_array_{i}', data=np.array(MU_Model[0][i][4]))
else:
file.create_dataset(f'cell_t_array_{i}', data=[0])
file.create_dataset(f'ina_array_{i}', data=np.array(MU_Model[0][i][6][0]))
stop = time.perf_counter()
print(f"elapsed time: {stop - start} sec")
plt.plot(MU_Model[0][0][0][0],MU_Model[0][0][2][0])
plt.show()