Python/Neuron model runs on Windows but does not on Linux.(h.finitialize)
Posted: Mon Jul 22, 2024 12:08 pm
I have a python/Neuron model which runs on windows but will not run on our server. No error flags are presented and we are using the multiprocessing.pool command. There seems to be an issue with the h.finitialize command where it has an issue with its children processes. Any help would be greatly appreciated.
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()