Conditional termination of simulation

[pascal]

I am simulating a single neuron with adaptation currents, and I would like to run a simulation until the most recent inter-spike interval differs from the second-most-recent inter-spike interval by less than some prescribed threshold. I have basic code that records spike times, so I'm wondering about the best approach to get NEURON to monitor the spike times and terminate the simulation using this condition. Thank you.

[ted]

Use events. Specifically,
1. use an FInitializeHandler to set the spike count to 0 at the start of the simulation
2. usethe NetCon class's record() method to call a function every time the cell spikes. This function should do the following:

Code: Select all

increase the spike count by 1
if spike count == 1
  tnow = t
  isinow = 1e9
if spike count > 1
  isiprev = isinow
  isinow = t - tnow
  if (abs(isiprev - isinow) < epsilon) stop the simulation
  tnow = t

[pascal]

Thanks for the help, Ted. I have to admit I don't understand FInitializeHandler very well even after reading its documentation page. In particular, I can't figure out how to initialize a variable 'spike_count' that is associated with the NetCon, and then connect the 'spike_count' variable defined for FInitializeHandler to 'spike_count' in the event callback function. I have tried many different variations of the following code, and I keep getting the error

spike_count = spike_count + 1
UnboundLocalError: local variable 'spike_count' referenced before assignment

Code: Select all

mycell = Cell()

nc = h.NetCon(mycell.soma(0.5)._ref_v, None, sec=mycell.soma)

def set_spike_count():
    spike_count = 0

fih = h.FInitializeHandler(set_spike_count)

tol = 0.001

def spike_callback():
    spike_count = spike_count + 1
    if spike_count == 1:
        tnow = h.t
        isinow = 1e9
    elif spike_count > 1:
        isiprev = isinow
        isinow = h.t - tnow
        if abs( (isinow-isiprev) / isinow) < tol:
            h.stoprun = 1
        tnow = h.t

nc.record(spike_callback)

mycell.createIClamp(1,9000,1000)

h.load_file('stdrun.hoc')
h.finitialize(-65)
h.continuerun(10000)

I am also wondering whether the variables tnow, isinow, and isiprev are saved from event to event, or if I need to somehow manually save them and pass them as inputs to 'spike_callback.'

Thank you for the help.

[pascal]

If I change the function to FInitializeHandler to

def set_spike_count(spike_count=0):
print('set_spike_count: spike_count = %d' % spike_count)

and then replace

nc.record(spike_callback)

with

spike_count = 0
nc.record((spike_callback, spike_count))

then the program runs without error. But it doesn't terminate as I want it to, I think because spike_count is being hard-coded as equal to zero with every call to spike_callback. It seems I need to figure out how to define spike_count as a Python variable and get the callback to change its value.

[pascal]

I got this working, but the code is terribly ugly, and I'd appreciate any insight into how to clean it up. First, I defined a very simple mechanism with the following mod file:

Code: Select all

NEURON {
	SUFFIX sp
	RANGE spike_count, tnow, isinow, isiprev
}

ASSIGNED {
	spike_count
	tnow
	isinow
	isiprev
}

In the Cell template I inserted the 'sp' mechanism into the soma, then ran the following simulation code:

Code: Select all

from neuron import h
from cell_class import Cell
h.load_file('stdrun.hoc')

mycell = Cell()

nc = h.NetCon(mycell.soma(0.5)._ref_v, None, sec=mycell.soma)
nc.threshold = 0

def set_spike_count(sec):
    sec.spike_count_sp = 0
    print('set_spike_count: spike_count = %d' % sec.spike_count_sp)

fih = h.FInitializeHandler((set_spike_count,mycell.soma))

tol = h.dt

def spike_callback(sec, tol):
    sec.spike_count_sp = sec.spike_count_sp + 1
    if sec.spike_count_sp == 1:
        sec.tnow_sp = h.t
        sec.isinow_sp = 1e9
    elif sec.spike_count_sp > 1:
        sec.isiprev_sp = sec.isinow_sp
        sec.isinow_sp = h.t - sec.tnow_sp
        if abs(sec.isinow_sp-sec.isiprev_sp) < tol:
            h.stoprun = 1
        sec.tnow_sp = h.t

nc.record((spike_callback, (mycell.soma, tol)))

mycell.createIClamp(1,9000,1000)

h.finitialize(-65)
h.continuerun(10000)

Two questions:
1) Is it necessary to create a new NMODL mechanism to accomplish this? (It seems awfully clunky to do so.)
2) If so, is there a better implementation?

[ted]

Good try. There is no need to write any NMODL code--everything can be done with hoc or Python or a combination of the two. Here's an example that assumes the existence of a section called soma with 100 um2 surface area and hh membrane properties. Omitted is the "from neuron import h, gui" statement and the code that creates the section and the following items:
a RunControl panel with all default values except for tstop, which is set to 100 ms.
a "voltage axis Graph" (to show soma(0.5).v vs. t)
a PointProcessManager configured as an IClamp with delay 1 ms, duration 1e9 ms, and amplitude 0.01 nA

The very first thing to do, after importing modules, is to define the criterion for ISI convergence.

Code: Select all

# stop execution when |isi - previsi| < TOL
TOL = h.dt # could be < h.dt

where isi and previsi are the most recent interspike interval, and the ISI before that. Note that if TOL is significantly < h.dt, convergence will be to within roundoff error. Also note that, if subsequent statements change h.dt, you might also want to change the value of TOL before running a simulation.

Then include the code that sets up the model cell.

Next to do are (1) create variables that keep track of the number of spikes, the most recent ISI, and the ISI before that, and (2) create a procedure that will set these variables to appropriate initial values during each initialization (so we can use the RunControl and IClamp for interactive testing of our code).

Code: Select all

# set aside storage for these values:
count = 0 # running tally of spikes
isi = 0 # most recent isi
previsi = 0 # isi before that

# every run starts a new tally
def resetstuff():
  global count, isi, previsi
  count = 0
  # assign isi and previsi values that are way too big and differ by more than TOL
  isi = 1e3 # most recent isi
  previsi = 2e3 # isi before that

rc = h.FInitializeHandler(resetstuff)

Yes, I know, global is not as pythonic as some would like, but my aim is to quickly create something simple that works. Note that isi and previsi are initialized to values that are clearly bogus (could even be negative if you like).

The next part is easy. It seemed appropriate to change threshold for spike detection from the default (10 mV) to -10 mV. You might want to try an even more negative value.

Code: Select all

# detect spikes
nc = h.NetCon(soma(0.5)._ref_v, None)
nc.threshold = -10

Last but not least, the code that tests for ISI convergence and stops simulation execution.

Code: Select all

tsp = 0

# act when a spike happens
def spikehappened():
  global tsp, count, isi, previsi
  count += 1
  if count == 1 :
    tsp = h.t
    # leave isi and previsi with their initialized values
    print("tsp, count, isi, previsi")
  elif count > 1 :
    previsi = isi
    isi = h.t - tsp
    tsp = h.t
  print("%7.3f %2d %7.3f %7.3f" % (tsp, count, isi, previsi))
  if abs(isi - previsi) < TOL: h.stoprun = 1
  # or stop a couple ms later so the last spike is captured
  # writing code to do that is left as an exercise to the reader

nc.record(spikehappened)

[pascal]

Thanks, Ted, that worked like a charm! I hadn't thought to use global variables, and that is a perfectly acceptable solution for my purposes. Appreciate the help.

[ted]

Well, "crude but effective" is still effective.

An alternative is to explicitly pass all variables whose values are to be changed as arguments, and return them as results, like so:

Code: Select all

a = somevalue
b = someothervalue
c = yetanothervalue

def foo(a, b, c):
  . . . statements that change a, b, and c . . .
  return a, b, c

a, b, c = foo(a, b, c)

which is arguably less prone to bugs in the future because the values of a, b, and c aren't altered as "side-effects" of calling foo.

[pascal]

I decided I couldn't live with myself if I didn't try to implement the more acceptable approach you suggested (passing all variables as arguments, and then returning them). However, I'm running into some trouble with syntax.

I can change spikehappened to the following:

Code: Select all

def spikehappened(tsp, count, isi, previsi):
  count += 1
  if count == 1 :
    tsp = h.t
    # leave isi and previsi with their initialized values
    print("tsp, count, isi, previsi")
  elif count > 1 :
    previsi = isi
    isi = h.t - tsp
    tsp = h.t
  print("%7.3f %2d %7.3f %7.3f" % (tsp, count, isi, previsi))
  if abs(isi - previsi) < TOL: h.stoprun = 1
  return tsp, count, isi, previsi

but then how do I adapt the line nc.record(spikehappened)? I can specify the input arguments like using nc.record((spikehappened,(tsp, count, isi, previsi))), but how do I specify the outputs?

It also occurred that I don't understand why FInitializeHandler and resetstuff are necessary in the following code to initialize count, isi, and previsi:

Code: Select all

# set aside storage for these values:
count = 0 # running tally of spikes
isi = 0 # most recent isi
previsi = 0 # isi before that

# every run starts a new tally
def resetstuff():
  global count, isi, previsi
  count = 0
  # assign isi and previsi values that are way too big and differ by more than TOL
  isi = 1e3 # most recent isi
  previsi = 2e3 # isi before that

rc = h.FInitializeHandler(resetstuff)

Naively it seems to me that count, isi, and previsi are Python variables, and the first three lines of the above code seem sufficient to initialize them (assuming we used the values in resetstuff). I don't see anything in the documentation for FInitializeHandler that helps me understand why exactly this is necessary.

Thank you!

[ted]

how do I adapt the line nc.record(spikehappened)? I can specify the input arguments like using nc.record((spikehappened,(tsp, count, isi, previsi))), but how do I specify the outputs?

Excellent question. I don't have an answer for it. If there is such an answer, maybe someone more knowledgable about Python can tell us what it is.

I don't understand why FInitializeHandler and resetstuff are necessary . . . to initialize count, isi, and previsi

The custom initialization facilitates interactive usage, which is one of the main benefits of having a GUI. Specifically, it enables this

Code: Select all

REPEAT
  tweak a parameter
  run a simulation
  observe and interpret results
UNTIL some insight has been obtained

In this particular example, the named variables must be properly initialized before each run. Absent a custom initialization, one would have to

Code: Select all

REPEAT
  revise source code
  save to a file
  use NEURON to execute that file
  observe and interpret results
  exit NEURON
UNTIL some insight has been obtained

which breaks the user's focus by introducing several manual steps and imposing delays for model setup and program exiting.

[pascal]

Thanks, Ted, that makes sense that this would be useful for the purposes of a GUI. Maybe Robert can help with the other question...