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)