sl6609 wrote:Say I have a neuron with a synaptic point process. I also have a Vector with recorded voltage levels that I would like to input into the synapse. Each element of the Vector contains a voltage value separated by a given dt.
If your model involves several different presynaptic waveforms that have to "drive" several synaptic mechanisms, the most direct and efficient approach is to first use a program that does the following:
1. analyzes a recorded voltage waveform to find the times at which it crosses your assumed threshold for spike-triggered synaptic transmission
2. saves these times to a file
Then write your simulation code so that it
1. sets up the model cell or network
2. reads the spike time file into a Vector
3. uses that Vector and an instance of the VecStim class to deliver events at those times to your synaptic mechanism
Source code for the VecStim class is in nrn-x.x/examples/nrniv/netcon/ where you will find vecevent.mod plus a hoc and ses file that demonstrates how it works. Here is a program that reads and analyzes voltage data, detects threshold crossing times and appends them to a Vector, and confirms its performance by plotting v vs. t and marking the times at which v rises to or above a user-specified threshold:
Code: Select all
/*
dattospikes.hoc
reads data in the format
first line
number_of_data_pairs
t0 v0
t1 v1
t2 v2
etc.
*/
load_file("nrngui.hoc")
THRESH = -45 // threshold for "spike detection"
strdef tstr
objref dfil, tdata, vdata
tdata = new Vector()
vdata = new Vector()
dfil = new File()
dfil.ropen("v.dat")
// first line is label:str
dfil.gets(tstr)
// second line is number of data pairs
numpairs = 0
numpairs = dfil.scanvar()
// subsequent lines are tval vval pairs
tdata.resize(numpairs)
vdata.resize(numpairs)
for i=0,numpairs-1 {
tdata.x[i] = dfil.scanvar()
vdata.x[i] = dfil.scanvar()
}
///// find where vdata rises above thresh
thresh = THRESH
objref tvdata, ttdata, stvec
ttdata = tdata.c
tvdata = vdata.c
stvec = new Vector() // will hold the times of threshold crossings
/*
The algorithm:
repeat
find the first datum >= threshold
record the corresponding time
discard data up to that point
look for the first datum < threshold
discard data up to that point
until there are no more data
Before entering this loop, one must be at the first datum < threshold.
*/
// $o1 and $o2 are the input time and data vectors
// returns a vector of "spike times"
obfunc spiketimes() { local i localobj tmpvec, ttmp, vtmp
ttmp = $o1.c
vtmp = $o2.c
tmpvec = new Vector()
i = vtmp.indwhere("<", thresh) // find first value < thresh
// print "below threshold at ", ttmp.x[i]
while (i>=0) { // look for the next spike time
// ith datum is subthreshold
vtmp.remove(0,i) // discard previous data
ttmp.remove(0,i)
i = vtmp.indwhere(">=", thresh) // find first value >= thresh
// print "crosses threshold at ", ttmp.x[i]
if (i>=0) {
// ith datum is at or above threshold
tmpvec.append(ttmp.x[i])
vtmp.remove(0,i) // discard previous data
ttmp.remove(0,i)
i = vtmp.indwhere("<", thresh) // find first value < thresh
// ith datum is below threshold
// print "below threshold at ", ttmp.x[i]
}
}
return tmpvec
}
stvec = spiketimes(tdata,vdata)
print "Spike Times"
stvec.printf
// seeing is believing
objref g
g = new Graph()
vdata.plot(g, tdata) // show the data
// mark the spike times
objref yval
yval = new Vector(stvec.size())
yval.fill(thresh)
yval.mark(g, stvec, "|", 12, 2, 1) // red vertical line
objref xx,yy
xx = new Vector(2)
yy = new Vector(2)
xx.x[0] = 0
xx.x[1] = 100
yy.fill(thresh)
yy.plot(g, xx, 2, 1) // red horizontal line
g.exec_menu("View = plot")
All you have to do is add a procedure that writes the spike times to a file.
If you are only going to have a single synapse that is driven by just one presynaptic voltage, it is quicker and easier to use Vector.play to drive your voltage waveform into a range variable that belongs to a point process that has a WATCH statement that generates an event whenever the range variable crosses threshold. Here is the source code for such a point process
Code: Select all
COMMENT
watchstim.mod
Play a Vector into WatchStim.x, and this mechanism will
launch an event when x rises above thresh.
Note: WATCH statements don't work with ARTIFICIAL_CELL.
ARTIFICIAL_CELLs are not associated with a cvode instance,
and this makes it problematic to associate them with the
proper cvode instance and thread.
ENDCOMMENT
NEURON {
POINT_PROCESS WatchStim
RANGE vv, thresh
}
PARAMETER { thresh }
ASSIGNED { vv } : the WATCHed variable must not be called "x"
INITIAL {
net_send(0, 1) : launch a self-event with flag==1
:to execute the WATCH statement
}
NET_RECEIVE(w) {
if (flag == 1) {
: if vv crosses thresh in a positive-going direction,
: send a self-event with flag==2
WATCH (vv > thresh) 2
} else if (flag == 2) {
: send an event to all NetCons that have
: this point process as their source
net_event(t)
}
}
and here is a demonstration program that drives it with a ramp waveform that crosses the "spike threshold" halfway through the simulation run
Code: Select all
load_file("nrngui.hoc")
create dummy // WatchStim is a POINT_PROCESS so needs a "host" section
objref ws
dummy ws = new WatchStim(0.5)
// drive ws.vv with a Vector's contents
objref vvec, tvec
tvec = new Vector(2)
tvec.x[0] = 0
tvec.x[1] = tstop
vvec = tvec.c
vvec.play(&ws.vv, tvec, 1) // "continuous" Vector.play
// vvec is treated as a continuous ramp that rises from 0 at t==0
// to tstop at t==tstop
ws.thresh = tstop/2 // thresh crossing halfway through the simulation
// report the threshold crossings
proc reportx() {
print "crossed threshold at t = ", t, " v = ", ws.vv
}
objref ncws, nil
ncws = new NetCon(ws, nil)
ncws.record("reportx()")
run()