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: 5704
 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: 196
 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: 5704
 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: 196
 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: 5704
 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: 196
 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: 196
 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
}