Okay, now I'm attempting to incorporate my firing rate model within a toy parallel network model. I started with the parallel ring network outlined here:
viewtopic.php?f=12&t=4270
I then added a dummy section (for the firing rate model) to the 0th host. In general, I want the activity of the firing rate model to influence the ionic conductances of all neurons in the network (on all hosts). I used a callback to implement the same silly example as earlier in this forum post, where if the activity of the firing rate model is above 0.5, then the sodium conductances get shut off throughout the whole network. (Obviously not biophysically realistic, but just a proof of principle demo.) I also used the callback to compute a very crappy "LFP" of the network (really just summed voltage of all segments). In principle, I'd also like the LFP to influence the activity of the firing rate model (though that is not incorporated in this toy model).
The code does exactly what I expect when I run it in serial: the network is shut down initially because the FR model's activity starts at 10. Then when the activity decays below 0.5 (after 50 or 60 ms), the sodium channels are turned on and the network starts firing.
However, when I try to run the program in parallel, it just hangs. I'm not really sure how to debug it because I can't get print statements (of the form
if idhost==0: print("Made it.") ) to work. Here's the program:
Code: Select all
from neuron import h
from matplotlib import pyplot as plt
import json
import numpy as np
h.load_file('stdrun.hoc')
class Pyramidal:
def __init__(self, gid):
self._gid = gid
self._createSections()
self._defineBiophysics()
self._register_netcon()
def __repr__(self):
return "Pyramidal[%d]" % self._gid
def _createSections(self):
self.soma=h.Section(name='soma', cell=self)
self.dend=h.Section(name='dend', cell=self)
self.dend.connect(self.soma(1))
self.all=self.soma.wholetree()
def _defineBiophysics(self):
self.soma.insert('hh')
self.dend.insert('pas')
def _register_netcon(self):
nc = h.NetCon(self.soma(0.5)._ref_v, None, sec=self.soma)
nc.threshold = -20 #because Hodgkin-Huxley spikes do not always reach 0mV
pc.set_gid2node(self._gid, idhost)
pc.cell(self._gid, nc)
del nc # discard netcon, since its only purpose is to detect spikes, and the previous line of code will now make sure that happens
class Network:
def __init__(self,num):
mygids= list(range(idhost,num,nhost))
self.cells = [Pyramidal(i) for i in mygids]
#set up an excitable ExpSyn on each cell's dendrites
self.syns = [h.ExpSyn(cell.dend(0.5)) for cell in self.cells]
for syn in self.syns:
syn.e = 0.0
#connect cell (i-1)%num to cell i
self.ncs = []
for i, syn in zip(mygids,self.syns):
nc = pc.gid_connect((i-1) % num, syn) #note that cell (i-1)%num is in general NOT on this host (when run in parallel)
nc.delay = 1.0
nc.weight[0] = 10.0
self.ncs.append(nc)
#set up spike recording
self.tVec = h.Vector() # spike time of all cells on this processor
self.idVec = h.Vector() # cell ids of spike times
for i in mygids:
pc.spike_record(i, self.tVec, self.idVec) # Record spikes of this cell
#create dummy section to house firing rate model on idhost=0
if idhost==0:
self.dummy_section = h.Section(name='dummy_section')
self.dummy_section.insert('frate')
pc = h.ParallelContext()
pc.set_maxstep(10)
idhost = int(pc.id())
nhost = int(pc.nhost())
n = Network(10)
#drive the 0th cell
if pc.gid_exists(0):
#constant stimulation to 0th cell, to elicit regular firing from it
stim = h.IClamp(pc.gid2cell(0).soma(0.5))
stim.dur = 1e9
stim.amp = 75.0
stim.delay = 10.0
t = h.Vector().record(h._ref_t)
r = h.Vector().record(n.dummy_section(0.5)._ref_f_frate)
v = [h.Vector().record(cell.soma(0.5)._ref_v) for cell in n.cells]
v_rec=[] #to record voltage summed over all segments
def callback():
if idhost==0: print("Made it into callback.")
#add up voltages from all segments
vsum = 0
for cell in n.cells:
for sec in cell.all:
for seg in sec:
vsum = vsum + seg.v
v_rec.append(vsum)
#send frate from host 0 to every other host
pc.barrier()
gather = pc.py_alltoall([n.dummy_section(0.5).f_frate]*nhost if idhost == 0 else [None]*nhost)
pc.barrier()
f=gather[0]
if(f>0.5):
for cell in n.cells:
cell.soma.gnabar_hh = 0.0 #if firing rate gets too high, it shuts off sodium current in HH cell, so cells won't fire; right now, frate.mod starts with activity of 10, so network will be shut off at beginning of simulation, until activity decays below 0.5
else:
for cell in n.cells:
cell.soma.gnabar_hh = 0.12 #standard value in HH model
h.cvode.extra_scatter_gather(0,callback)
h.stdinit()
pc.psolve(100)
#for myv in v:
# plt.plot(t,myv)
#plt.xlabel('Time (ms)')
#plt.ylabel('Voltage (mV)')
#plt.show()
#
#plt.figure()
#plt.plot(t, r, label='firing rate')
"""Gather spikes from all nodes/hosts"""
data = [None]*nhost
data[0] = {'tVec': n.tVec, 'idVec': n.idVec}
pc.barrier()
gather=pc.py_alltoall(data)
pc.barrier()
tVecAll = []
idVecAll = []
if idhost==0:
for d in gather:
tVecAll.extend(list(d['tVec']))
idVecAll.extend(list(d['idVec']))
'''Gather summed voltage from all nodes/hosts'''
data = [None]*nhost #EACH NODE has this list, the i^th element of which will be sent to node i
data[0] = {'v_rec': v_rec} #by making only the zeroth element something other than 'None,' this means each node will be sending data only to node 0
pc.barrier()
gather=pc.py_alltoall(data)
pc.barrier()
if idhost==0:
net_v_rec=np.zeros(len(gather[0]['v_rec'])) #start sum at zeros, and make np array same length as v_rec lists
for d in gather:
net_v_rec += d['v_rec'] #compute summed voltage, summed over contributions from segments on all hosts
#dump raster data and "LFP" to JSON file
if idhost==0:
raster_dat = {'tVecAll': tVecAll, 'idVecAll': idVecAll}
with open("raster.json", "w") as outfile:
json.dump(raster_dat, outfile)
v_dat = {'net_v_rec': net_v_rec.tolist()}
with open("v_rec.json", "w") as outfile:
json.dump(v_dat, outfile)
pc.barrier()
h.quit()
Here is frate.mod:
Code: Select all
COMMENT
This mod file implements a simple firing rate model to represent the firing rate of
a population of cells
ENDCOMMENT
NEURON {
SUFFIX frate
RANGE tau, decay, f, f_of_i_of_t
}
PARAMETER {
tau = 10.0 (ms)
decay = 0.5
}
ASSIGNED {
f_of_i_of_t
}
STATE {
f
}
INITIAL {
f = 10.0 :initial activity of this FR model
f_of_i_of_t = 0.0
}
BREAKPOINT {
SOLVE states METHOD cnexp
}
Here's a little routine to make the raster plot:
Code: Select all
from matplotlib import pyplot as plt
import json
# Opening JSON file
with open('raster.json', 'r') as openfile:
# Reading from json file
raster_dict = json.load(openfile)
stimes = raster_dict['tVecAll']
neurid = raster_dict['idVecAll']
plt.figure()
for i in range(len(stimes)):
plt.vlines(stimes[i], neurid[i]+0.5, neurid[i]+1.5)
plt.show()
And here's a little routine to plot the "LFP":
Code: Select all
from matplotlib import pyplot as plt
import json
import numpy as np
dt = 0.025 # time step for simulation
# Opening JSON file
with open('v_rec.json', 'r') as openfile:
# Reading from json file
v_dict = json.load(openfile)
net_v_rec = v_dict['net_v_rec']
plt.figure()
plt.plot(dt*np.arange(0,len(net_v_rec)),net_v_rec)
plt.show()
I have three questions:
1) Is this even the correct way to approach parallelizing this program? For example, I assume I should make a dummy section on just one host (might as well be host zero) for the FR model, but perhaps I should make identical copies on all hosts? The problem is that I'd like to eventually have the network activity influence the FR model, with some stochasticity, and it sounds like a nightmare trying to get each copy of the FR model to be identical across all hosts.
2) Is it even worth parallelizing this model? Since it requires all the hosts to talk to each other on every time step, I suspect it may not be. Unfortunately, I'm unable to get the program to run in parallel, so I can't compare serial versus parallel performance.
3) If it *is* worth parallelizing, then what is wrong with my code, and/or how do I go about debugging it in parallel?
Thanks for the help!