Incorporating a firing model to induce different stages of sleep
Moderator: wwlytton
Incorporating a firing model to induce different stages of sleep
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 momenttomoment 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!
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 momenttomoment 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!

 Site Admin
 Posts: 5869
 Joined: Wed May 18, 2005 4:50 pm
 Location: Yale University School of Medicine
 Contact:
Re: Incorporating a firing model to induce different stages of sleep
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.
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.
The contrary would surprise me.Surprisingly, I haven't been able to find much on implementing firing rate models in NEURON.
Pointers and fictive concentrations come to mind. You might glean some ideas from the implementations of dopamine effects in https://modeldb.yale.edu/399492) How do I then connect . . .

 Posts: 223
 Joined: Fri Nov 28, 2008 3:38 pm
 Location: Yale School of Public Health
Re: Incorporating a firing model to induce different stages of sleep
You can also do this using the LinearMechanism class.
LinearMechanism solves vector ODEs of the form:
In particular, let's consider the simple firing rate model:
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:
LinearMechanism solves vector ODEs of the form:
Code: Select all
c y' + gy = b
Code: Select all
tau * y' + decay * y = f(I(t))
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()
Re: Incorporating a firing model to induce different stages of sleep
Thank you both very much! That will certainly get me started.

 Site Admin
 Posts: 5869
 Joined: Wed May 18, 2005 4:50 pm
 Location: Yale University School of Medicine
 Contact:
Re: Incorporating a firing model to induce different stages of sleep
Thanks for posting an excellent question.
Re: Incorporating a firing model to induce different stages of sleep
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 differentialalgebraic 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 timestep 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?
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 differentialalgebraic 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 timestep 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?

 Posts: 223
 Joined: Fri Nov 28, 2008 3:38 pm
 Location: Yale School of Public Health
Re: Incorporating a firing model to induce different stages of sleep
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:
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.
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:
In particular, the first equation in the example in the documentation is y[1] = 0... but this isn't solved as a standalone 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.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.
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.

 Site Admin
 Posts: 5869
 Joined: Wed May 18, 2005 4:50 pm
 Location: Yale University School of Medicine
 Contact:
Re: Incorporating a firing model to induce different stages of sleep
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 angleit's still a bottle of ink. When the need to add my own ODEs arises, I resort to NMODL.
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 angleit's still a bottle of ink. When the need to add my own ODEs arises, I resort to NMODL.
Re: Incorporating a firing model to induce different stages of sleep
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:
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...
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()

 Posts: 223
 Joined: Fri Nov 28, 2008 3:38 pm
 Location: Yale School of Public Health
Re: Incorporating a firing model to induce different stages of sleep
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.
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.
Re: Incorporating a firing model to induce different stages of sleep
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:
Thanks for the help!
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()
Re: Incorporating a firing model to induce different stages of sleep
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:
And here is the mod file:
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()
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
}

 Posts: 223
 Joined: Fri Nov 28, 2008 3:38 pm
 Location: Yale School of Public Health
Re: Incorporating a firing model to induce different stages of sleep
Dereference a NEURON pointer using [0]. e.g.
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.
Code: Select all
>>> from neuron import h
>>> h.t = 17
>>> print(h._ref_t[0])
17.0
Re: Incorporating a firing model to induce different stages of sleep
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
}
Re: Incorporating a firing model to induce different stages of sleep
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:
Here is frate.mod:
Here's a little routine to make the raster plot:
And here's a little routine to plot the "LFP":
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!
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 HodgkinHuxley 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 (i1)%num to cell i
self.ncs = []
for i, syn in zip(mygids,self.syns):
nc = pc.gid_connect((i1) % num, syn) #note that cell (i1)%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()
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
}
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()
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()
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!