Bug in Exp2Syn?
Posted: Wed Jul 12, 2017 3:28 am
When I was playing around with synapses in a neuron model, I figured out that the start of current change in Exp2Syn is delayed by 1 time step (which is not the case in ExpSyn). See the code below for a minimal example. I found this behavior strange and found out that A and B in Exp2Syn already change at the previous time step. So if I would compute the current after the simulation: i_syn = (A-B) * (V-E) then the resulting current trace would be shifted by one time step and start where I would expect it. Is this a bug? Why is this happening?
Another point: In the documentation of Exp2Syn it says: "if tau1 -> 0 then we have just single exponential decay." But setting tau1 = 0 in the code below makes the simulation stop at 100 ms (start of the synaptic event) without telling an error. I looked at the code and found that no care was taken for the special case tau1=0. It basically results in a division by 0 when computing the derivative of A: A' = -A/tau1.
For ExpSyn and Exp2Syn I refer to the code found here: http://neuron.yale.edu/hg/neuron/nrn/fi ... expsyn.mod and http://www.neuron.yale.edu/hg/neuron/nr ... xp2syn.mod
Another point: In the documentation of Exp2Syn it says: "if tau1 -> 0 then we have just single exponential decay." But setting tau1 = 0 in the code below makes the simulation stop at 100 ms (start of the synaptic event) without telling an error. I looked at the code and found that no care was taken for the special case tau1=0. It basically results in a division by 0 when computing the derivative of A: A' = -A/tau1.
For ExpSyn and Exp2Syn I refer to the code found here: http://neuron.yale.edu/hg/neuron/nrn/fi ... expsyn.mod and http://www.neuron.yale.edu/hg/neuron/nr ... xp2syn.mod
Code: Select all
from neuron import h
import numpy as np
import pylab as pl
h.load_file("stdrun.hoc") # load NEURON libraries
h("""cvode.active(0)""") # unvariable time step in NEURON
# create simple cell model
cell = h.Section(name='soma')
cell.cm = 1
cell.Ra = 100
cell.diam = 100
cell.L = 100
cell.nseg = 1
cell.insert('pas')
# insert two synapses
syn1 = h.ExpSyn(.5, sec=cell)
syn1.tau = 1 # ms
syn1.e = 0 # mV
syn2 = h.Exp2Syn(.5, sec=cell)
syn2.tau1 = 1 # ms
syn2.tau2 = 2 # ms
syn2.e = 0 # mV
h.nrn_load_dll("/your/path/to/vecevent/x86_64/.libs/libnrnmech.so") # http://www.neuron.yale.edu/hg/neuron/nrn/file/b30234a6f319/share/examples/nrniv/netcon/vecevent.mod
spiketimes = [100] # ms
spiketimes_vec = h.Vector()
spiketimes_vec.from_python(spiketimes) # convert spiketimes to neuron vector
# make stimulus
stim = h.VecStim(.5, sec=cell)
stim.play(spiketimes_vec)
con1 = h.NetCon(stim, syn1, sec=cell)
con1.weight[0] = 1
con1.delay = 0
con2 = h.NetCon(stim, syn2, sec=cell)
con2.weight[0] = 1
con2.delay = 0
# recordings
i_syn1 = h.Vector()
i_syn1.record(syn1._ref_i)
i_syn2 = h.Vector()
i_syn2.record(syn2._ref_i)
# run simulation
h.tstop = 200
h.dt = 0.1
h.init()
h.run()
# change to arrays
i_syn1 = np.array(i_syn1)
i_syn2 = np.array(i_syn2)
t = np.arange(0, h.tstop+h.dt, h.dt)
# plot
pl.figure()
pl.plot(t, i_syn1, 'r')
pl.plot(t, i_syn2, 'b')
pl.show()