Incorporating a firing model to induce different stages of sleep

Moderator: wwlytton

Post Reply
pascal
Posts: 84
Joined: Thu Apr 26, 2012 11:51 am

Incorporating a firing model to induce different stages of sleep

Post by pascal » Thu Feb 27, 2020 7:33 pm

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
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

Post by ted » Fri Feb 28, 2020 11:15 am

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
Posts: 195
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

Post by ramcdougal » Fri Feb 28, 2020 11:37 am

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
Posts: 84
Joined: Thu Apr 26, 2012 11:51 am

Re: Incorporating a firing model to induce different stages of sleep

Post by pascal » Fri Feb 28, 2020 4:04 pm

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

ted
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

Post by ted » Fri Feb 28, 2020 5:32 pm

Thanks for posting an excellent question.

pascal
Posts: 84
Joined: Thu Apr 26, 2012 11:51 am

Re: Incorporating a firing model to induce different stages of sleep

Post by pascal » Mon Mar 02, 2020 12:50 pm

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
Posts: 195
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

Post by ramcdougal » Tue Mar 03, 2020 11:27 am

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
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

Post by ted » Tue Mar 03, 2020 11:46 am

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
Posts: 84
Joined: Thu Apr 26, 2012 11:51 am

Re: Incorporating a firing model to induce different stages of sleep

Post by pascal » Tue Mar 03, 2020 1:32 pm

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
Posts: 195
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

Post by ramcdougal » Tue Mar 03, 2020 1:53 pm

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
Posts: 84
Joined: Thu Apr 26, 2012 11:51 am

Re: Incorporating a firing model to induce different stages of sleep

Post by pascal » Wed Mar 04, 2020 12:53 pm

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
Posts: 84
Joined: Thu Apr 26, 2012 11:51 am

Re: Incorporating a firing model to induce different stages of sleep

Post by pascal » Wed Mar 11, 2020 4:27 pm

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
Posts: 195
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

Post by ramcdougal » Sun Mar 15, 2020 11:32 am

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
Posts: 84
Joined: Thu Apr 26, 2012 11:51 am

Re: Incorporating a firing model to induce different stages of sleep

Post by pascal » Sun Mar 15, 2020 7:03 pm

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
 }

Post Reply