Incorporating a firing model to induce different stages of sleep

[pascal]

I am developing a thalamocortical network model of sleep. Currently, various cortical ionic conductances (e.g. potassium currents) are altered by a "hand of God" that sets them to different values to induce different stages of sleep. I would like to develop the model so that the ionic conductances are determined by neuromodulator concentrations (such as acetylcholine), which are in turn determined by a fluctuating firing rate model (representing, for example, the activity of cholinergic neurons in the basal forebrain). The "hand of God" would then directly change the mean firing rate of (for example) the cholinergic population, which would then change ACh concentration, which would affect cortical ionic conductances. The firing rates would have some stochasticity, so that neuromodulator concentrations and ionic conductances would also fluctuate over time, even within the same sleep state.

So I am wondering the following:
1) How do I incorporate a firing rate model within NEURON? I want to specify one differential equation that represents the activity of an entire population of neurons (for the sake of computational efficiency). Surprisingly, I haven't been able to find much on implementing firing rate models in NEURON.

2) How do I then connect these firing rates to determine moment-to-moment ionic conductances? (with the neuromodulator concentration just being an intermediate calculation involved in that connection)

I'm just looking for a push in the right direction. If there are any examples of something similar on ModelDB you could point me to, that would be excellent. Thanks!

[ted]

ODEs can be implemented in mod files. A POINT_PROCESS can be used to add a bag of equations to NEURON. State variables are automatically RANGE variables, visible to hoc and Python. ASSIGNED variables need to be specified as RANGE or GLOBAL to be visible to hoc and Python; RANGE is the more flexible choice and avoids potential conflicts with multithreaded or MPI parallel simulation.

Events can be used to construct state machines. The FInitializeHandler class, NetStim and other artificial spiking cell classes, NetCon.event and cvode.event methods, and NMODL's NET_RECEIVE blocks can be helpful.

Surprisingly, I haven't been able to find much on implementing firing rate models in NEURON.

The contrary would surprise me.

2) How do I then connect . . .

Pointers and fictive concentrations come to mind. You might glean some ideas from the implementations of dopamine effects in https://modeldb.yale.edu/39949

[ramcdougal]

You can also do this using the LinearMechanism class.

LinearMechanism solves vector ODEs of the form:

Code: Select all

  c y' + gy = b

In particular, let's consider the simple firing rate model:

Code: Select all

    tau * y' + decay * y = f(I(t))

For this example, we suppose that I(t) and hence f(I(t)) are known in advance; this allows us to precompute them and just play into the vector b.

But these can be updated at runtime based on computed states by adding a callback as the first argument to LinearMechanism that updates the vector b. See the docs (linked above) for more. Alternatively, you could combine this approach and Ted's approach using a MOD file to update the vector b based on simulated values.

Anyways, here's working code for the above firing rate model:

Code: Select all

from neuron import h
from neuron.units import ms
import numpy as np
import matplotlib.pyplot as plt
h.load_file('stdrun.hoc')

tau = 10
decay = 0.5

# define storage, matrices, and initial values
y0 = h.Vector([10])
y = h.Vector([0])
b = h.Vector([0])

c = h.Matrix(1, 1)
g = h.Matrix(1, 1)
g.setval(0, 0, decay)
c.setval(0, 0, tau)

# calculate how the right hand side b should evolve
sample_t = h.Vector(np.linspace(0, 50, 500))
f_of_i_of_t = h.Vector(4 + 4 * np.sin(sample_t))

# link b[0] to the above time series; True means interpolate as needed
f_of_i_of_t.play(b._ref_x[0], sample_t, True)

# now we declare the mechanism we have defined
mech = h.LinearMechanism(c, g, y, y0, b)

# for reasons (has to do with threads), we need at least one section
ignored_section = h.Section(name='ignored_section')

# setup Vectors to store the results
t = h.Vector().record(h._ref_t)
rate = h.Vector().record(y._ref_x[0])

# run the simulation
h.finitialize()
h.continuerun(50 * ms)

plt.plot(t, rate, label='firing rate')
plt.plot(sample_t, f_of_i_of_t, label='f(I(t))')
plt.xlabel('t (ms)')
plt.legend()
plt.show()

[pascal]

Thank you both very much! That will certainly get me started.

[ted]

Thanks for posting an excellent question.

[pascal]

The LinearMechanism class looks very appealing, but I want to make sure I thoroughly understand it. A few questions:

1. I noticed in the documentation for LinearMechanism, there is the following warning:
"Does not work with the CVODE integrator but does work with the differential-algebraic solver IDA. Note that if the standard run system is loaded, h.cvode_active(True) will automatically choose the correct variable step integrator."
I plan to implement this in a parallelized network model (with fixed time-step integration). Does the above warning imply that the entire simulation will have to switch over to the IDA solver, even if I only use LinearMechanism to implement a simple firing rate model that doesn't involve any neural tree structure?

2. I am trying to understand how Example #1 in the documentation implements a voltage clamp. The 'c' matrix is zero'd out, the 'g' matrix is essentially np.array([[0,-1], [1,0]]), and 'b' is the vector [0,10]'. This gives the two equations -y[1]=0 and y[0]=10. I get that y[0] is interpreted to be the voltage, so it makes sense you'd set it equal to 10 if you want to clamp the voltage to 10mV. But if y[1] is the injected current, I don't understand how it ever assumes a value other zero, since it appears the equations set it equal to zero. I must be missing something really basic.

3. Related to number question #2, how do you know y[1] is the injected current? If you had three dependent variables, what physical entity would y[2] correspond to?

[ramcdougal]

For (1):

That's a statement about the variable step solver and does not directly apply to using the fixed step solver. However: above the example, there's also a statement about the fixed step solver, which LinearMechanism forces to use a matrix solver by Kundert (because adding arbitrary equations loses the special matrix structure that allows the fast solver to work). My guess is that in an MPI context, this only affects whatever nodes contain the LinearMechanism, but that's only a guess.

For (2) and (3):

The relevant portion of the documentation is:

The scalar form, x, with the specified section means that the first equation is added to the current balance equation at this location and the first dependent variable is a copy of the membrane potential. If the system is coupled to more than one location, then sl must be a SectionList and xvec a Vector of relative positions (0 … 1) specifying the locations. In this case, the first xvec.size equations are added to the corresponding current balance equations and the first xvec.size dependent y variables are copies of the membrane potentials at this location.

In particular, the first equation in the example in the documentation is -y[1] = 0... but this isn't solved as a stand-alone equation; this is instead added to the current balance equations for soma(0.5)... i.e. (all other currents) - y[1] = 0. This is mathematically necessary as otherwise voltage couldn't be fixed... there needs to be some play to account for the current difference; here we call that current y[1]. In the general case, it couldn't be any of the first len(sl) variables, but it could be defined to be anything after that.

If you're only interested in firing rate models and aren't trying to directly couple to a voltage variable, then you just don't give an x and section or a sl and xvec... and then you have full control of what every variable means.

[ted]

My own two bits:

Regarding your first question, don't worry about it. The comment is a comment, not a warning. Adaptive integration is generally not useful in either serial or parallel simulations of networks, so you won't be using the CVODE integrator. And NEURON won't use IDA unless (1) its analysis of the system equations determines that it must use IDA, or (2) for whatever reason of your own, you deliberately execute statements that force it to (I'm not going to tell you what they are).

With regard to questions 2 and 3, maybe someone else can help you. To me, the documentation of LinearMechanism seems as transparent as a bottle of ink, and reading ramcdougal's comments about it is like looking at the bottle from a slightly different angle--it's still a bottle of ink. When the need to add my own ODEs arises, I resort to NMODL.

[pascal]

Okay, thanks for the replies. I understand better what's going on with LinearMechanism now. I took the toy model from a few replies ago, and I'm trying to tweak it by adding an HH soma and having a dialogue between it and the firing rate model. I added a callback to LinearMechanism that gives the FR model a kick when the HH soma's voltage is above zero, and that works fine. I also added a callback (using extra_scatter_gather) to try to have the FR model shut off the HH cell's sodium current above some firing rate threshold. (These mechanisms obviously are not biologically plausible, but I'm just trying to figure out how to get the FR model to influence a neuron and vice versa.)

Here's the code:

Code: Select all

from neuron import h
from neuron.units import ms
import numpy as np
import matplotlib.pyplot as plt
h.load_file('stdrun.hoc')

#HH neuron's membrane potential will influence firing rate model's activity
soma = h.Section()
soma.insert('hh')
stim = h.IClamp(soma(0.5))
stim.amp = 10
stim.delay = 5
stim.dur = 1000

tau = 10
decay = 0.5

# define storage, matrices, and initial values
y0 = h.Vector([0]) #initial firing rate is zero
y = h.Vector([0])
b = h.Vector([0])

c = h.Matrix(1, 1)
g = h.Matrix(1, 1)
g.setval(0, 0, decay)
c.setval(0, 0, tau)

def lm_callback():
    #while neuron is spiking, firing rate model gets a kick
    if soma(0.5).v > 0:
        b[0] = 10.0
    else: 
        b[0]=0

# now we declare the mechanism we have defined
mech = h.LinearMechanism(lm_callback, c, g, y, y0, b)

# setup Vectors to store the results
t = h.Vector().record(h._ref_t)
rate = h.Vector().record(y._ref_x[0])
soma_v = h.Vector().record(soma(0.5)._ref_v)

def gnabar_callback(frate):
    print("frate=",frate)
    if(frate>0.5):
        soma.gnabar_hh = 0.0 #if firing rate gets too high, it shuts off sodium current in HH cell
    else:
        soma.gnabar_hh = 0.12 #standard value in HH model

runtime_callback = (gnabar_callback,y[0])

# run the simulation
h.cvode.extra_scatter_gather(0,runtime_callback)
h.finitialize()
h.continuerun(100 * ms)

plt.plot(t, rate, label='firing rate')
#plt.plot(sample_t, f_of_i_of_t, label='f(I(t))')
plt.plot(t,soma_v,label='Neural voltage')
plt.xlabel('t (ms)')
plt.legend()
plt.show()

The problem is with the callback with extra_scatter_gather. Somehow the input to gnabar_callback, which is supposed to be y[0] (the firing rate of the FR model), is always equal to zero. Maybe this is some sort of scope issue? If so I'm not sure how to fix it...

[ramcdougal]

Without looking closely at this or directly answering your question, my immediate thought is to wonder if lm_callback gets a meaningful value of y?

You can avoid doing callbacks at every time step (and thus have a faster simulation) by using a StateTransitionEvent to do an action whenever a variable crosses a threshold.

[pascal]

Ah, got it. The solution is to pass the vector 'y' to gnabar_callback, then access the 0th element within the function itself. The way I was doing it before, I was essentially just passing a constant value to gnabar_callback. Here's the updated code:

Code: Select all

from neuron import h
from neuron.units import ms
import numpy as np
import matplotlib.pyplot as plt
h.load_file('stdrun.hoc')

#HH neuron's membrane potential will influence firing rate model's activity
soma = h.Section()
soma.insert('hh')
stim = h.IClamp(soma(0.5))
stim.amp = 10
stim.delay = 5
stim.dur = 1000

tau = 10
decay = 0.5

# define storage, matrices, and initial values
y0 = h.Vector([0]) #initial firing rate is zero
y = h.Vector([0])
b = h.Vector([0])

c = h.Matrix(1, 1)
g = h.Matrix(1, 1)
g.setval(0, 0, decay)
c.setval(0, 0, tau)

# calculate how the right hand side b should evolve
#sample_t = h.Vector(np.linspace(0, 50, 500))
#f_of_i_of_t = h.Vector(4 + 4 * np.sin(sample_t))

# link b[0] to the above time series; True means interpolate as needed
#f_of_i_of_t.play(b._ref_x[0], sample_t, True)


def lm_callback():
    #while neuron is spiking, firing rate model gets a kick
    if soma(0.5).v > 0:
        b[0] = 10.0
    else: 
        b[0]=0

# now we declare the mechanism we have defined
mech = h.LinearMechanism(lm_callback, c, g, y, y0, b)

# for reasons (has to do with threads), we need at least one section
#ignored_section = h.Section(name='ignored_section')

# setup Vectors to store the results
t = h.Vector().record(h._ref_t)
rate = h.Vector().record(y._ref_x[0])
soma_v = h.Vector().record(soma(0.5)._ref_v)

def gnabar_callback(frate): #teh way I've coded this, the input to this function is a vector with one element. So I need to access that element using "frate[0]"
    print("frate=",frate[0])
    if(frate[0]>0.5):
        soma.gnabar_hh = 0.0 #if firing rate gets too high, it shuts off sodium current in HH cell
    else:
        soma.gnabar_hh = 0.12 #standard value in HH model

runtime_callback = (gnabar_callback,y)
# run the simulation
h.cvode.extra_scatter_gather(0,runtime_callback)
h.finitialize()
h.continuerun(100 * ms)

plt.plot(t, rate, label='firing rate')
#plt.plot(sample_t, f_of_i_of_t, label='f(I(t))')
plt.plot(t,soma_v,label='Neural voltage')
plt.xlabel('t (ms)')
plt.legend()
plt.show()

Thanks for the help!

[pascal]

I decided to try to implement the program I last posted using mod files, rather than LinearMechanism. I almost have it, except for an issue with the callback function similar to what I encountered a few posts ago.

When I define the runtime_callback (to be used by extra_scatter_gather), I want to send it the current firing rate from the firing rate model. This value is defined by the variable dummy_section(0.5).f_frate. However, if I pass the variable itself, the callback function only ever sees its values as being zero, even though it is clearly not zero (as evidenced by the plot generated at the end of the simulation). I tried instead passing the pointer, dummy_section(0.5)._ref_f_frate. This is essentially how I solved the problem in the LinearMechanism implementation.

Here, however, I'm stuck because I don't know how to dereference a pointer in Python+NEURON (and I could not find any forum post on this issue). Is there an equivalent to C++'s *ptr ?

Here's the Python code:

Code: Select all

from neuron import h
from neuron.units import ms
import numpy as np
import matplotlib.pyplot as plt
h.load_file('stdrun.hoc')

#HH neuron's membrane potential will influence firing rate model's activity
soma = h.Section()
soma.insert('hh')
stim = h.IClamp(soma(0.5))
stim.amp = 10
stim.delay = 5
stim.dur = 1000

tau = 10
decay = 0.5

# dummy section to insert firing rate mechanism
dummy_section = h.Section(name='dummy_section')

dummy_section.insert('frate')
dummy_section.tau_frate = tau
dummy_section.decay_frate = decay

# setup Vectors to store the results
t = h.Vector().record(h._ref_t)
r = h.Vector().record(dummy_section(0.5)._ref_f_frate)
soma_v = h.Vector().record(soma(0.5)._ref_v)

def callback(frate): 
    #right now 'frate' is a pointer, and I need to figure out how to dereference it
    #print("frate=",frate)
    if(frate>0.5):
        soma.gnabar_hh = 0.0 #if firing rate gets too high, it shuts off sodium current in HH cell
    else:
        soma.gnabar_hh = 0.12 #standard value in HH model
    
    #while neuron is spiking, firing rate model gets a kick
    if soma(0.5).v > 0:
        dummy_section(0.5).f_of_i_of_t_frate = 10.0
    else: 
        dummy_section(0.5).f_of_i_of_t_frate = 0.0

#problem is here: if I pass 'dummy_section(0.5).f_frate', the input to the function is only ever zero (for reasons I don't understand); if I pass 'dummy_section(0.5)._ref_f_frate', I do not know how to dereference this pointer
runtime_callback = (callback,dummy_section(0.5)._ref_f_frate) 

# run the simulation
h.cvode.extra_scatter_gather(0,runtime_callback)
h.finitialize()
h.continuerun(100 * ms)

plt.plot(t, r, label='firing rate')
plt.plot(t,soma_v,label='Neural voltage')
plt.xlabel('t (ms)')
plt.legend()
plt.show()

And here is the mod file:

Code: Select all

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
	f_of_i_of_t = 0.0
}

BREAKPOINT {
	SOLVE states METHOD cnexp
}
 
 DERIVATIVE states {
	f' = (f_of_i_of_t - decay*f)/tau
 }
 

[ramcdougal]

Dereference a NEURON pointer using [0]. e.g.

Code: Select all

>>> from neuron import h
>>> h.t = 17
>>> print(h._ref_t[0])
17.0

This also works in C/C++; the way to think of it is that referring to an array is the same as referring to a block of memory.

[pascal]

Perfect, thanks Robert. I see that your solution is exactly what I used before in the LinearMechanism implementation, but I thought it worked then only because I was passing a vector as the argument (i.e., I didn't realize that was the general way to dereference a pointer). Thanks again, and here's the working code for anyone who might find it useful:

Code: Select all

from neuron import h
from neuron.units import ms
import numpy as np
import matplotlib.pyplot as plt
h.load_file('stdrun.hoc')

#HH neuron's membrane potential will influence firing rate model's activity
soma = h.Section()
soma.insert('hh')
stim = h.IClamp(soma(0.5))
stim.amp = 10
stim.delay = 5
stim.dur = 1000

tau = 10
decay = 0.5

dummy_section = h.Section(name='dummy_section')

dummy_section.insert('frate')
dummy_section.tau_frate = tau
dummy_section.decay_frate = decay

# setup Vectors to store the results
t = h.Vector().record(h._ref_t)
r = h.Vector().record(dummy_section(0.5)._ref_f_frate)
soma_v = h.Vector().record(soma(0.5)._ref_v)

def callback(frate): #teh way I've coded this, the input to this function is a vector with one element. So I need to access that element using "frate[0]"
    print("frate=",frate[0])
    if(frate[0]>0.5):
        soma.gnabar_hh = 0.0 #if firing rate gets too high, it shuts off sodium current in HH cell
    else:
        soma.gnabar_hh = 0.12 #standard value in HH model
    
    #while neuron is spiking, firing rate model gets a kick
    if soma(0.5).v > 0:
        dummy_section(0.5).f_of_i_of_t_frate = 10.0
    else: 
        dummy_section(0.5).f_of_i_of_t_frate = 0.0


runtime_callback = (callback,dummy_section(0.5)._ref_f_frate)
# run the simulation
h.cvode.extra_scatter_gather(0,runtime_callback)
h.finitialize()
h.continuerun(100 * ms)

plt.plot(t, r, label='firing rate')
plt.plot(t,soma_v,label='Neural voltage')
plt.xlabel('t (ms)')
plt.legend()
plt.show()

Code: Select all

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 = 0.0
	f_of_i_of_t = 0.0
}

BREAKPOINT {
	SOLVE states METHOD cnexp
}
 
 DERIVATIVE states {
	f' = (f_of_i_of_t - decay*f)/tau
 }

[pascal]

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!