Bug in Exp2Syn?

[cafischer]

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

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

[ted]

No bug. Advancing the solution from t to t+dt by one step using the default implicit Euler method begins by evaluating the currents at t. This calculation is done by executing the assignment statements in the BREAKPOINT block, and it uses the values of the state variables (including v) at time t. After that, voltage is updated. Consequenty, at the end of fadvance() (the Python name for which is h.fadvance()) the states have values that correspond to t+dt, but the currents have values that are consistent with v and the other state values at t. That's the source of the apparent 1 step delay in the Exp2Syn's current.

By the way, the same thing happens if the synapse is an instance of the ExpSyn class--at the end of each fadvance(), the current is consistent with the membrane potential and synaptic conductance at the previous time step. You just didn't notice it. If you want to check this for yourself, initialize the simulation to the model's resting potential, then h.continuerun to the time at which the first synaptic activation occurs. At that point print h.t, membrane potential, and the conductance and current belonging to one synapse. Then advance by a single step with h.fadvance(), and print out the new h.t, membrane potential, and synaptic conductance and current. Finally, calculate the actual synaptic current at the new time by evaluating the product of synaptic conductance and the driving force for synaptic current current flow.

If it is important for currents, voltages, and states to all be mutually consistent throughout a simulation, use adaptive integration ("cvode"). If you are using adaptive integration but want to record the solution at fixed time intervals, set up vector recording with interpolation. Do that with this syntax in Python
vdest.record(var_reference, Dt)
or
vdest.record(&var, Dt) in hoc. Be sure to read the Programmer's Reference entry on the Vector class's record method.

Another point: In the documentation of Exp2Syn it says: "if tau1 -> 0 then we have just single exponential decay."

Sorry about that. We should have written "as tau1 -> 0" with the intention that it would be read "as tau1 approaches 0". A correct and unambiguous statement would be "In the limit as tau1 becomes much smaller than tau2, we have single exponential decay."

no care was taken for the special case tau1=0

I suppose we could bulletproof this by doing something like this in the initial block

Code: Select all

if (tau1<1e-9) {
  tau1 = 1e-9
}

but that would deprive users of the beneficial exercise of debugging their own model specification code in situations where they set tau1<=0. Where's the fun in that?

[ted]

A revised exp2syn.mod has been pushed to github. The substantive changes are
1. The "If tau2-tau1" comments have been changed to

If tau2-tau1 is very small compared to tau1, this is an alphasynapse with time constant tau2.
If tau1/tau2 is very small, this is single exponential decay with time constant tau2.

2. In the INITIAL block, the following test has been inserted immediately after the test for tau1/tau2 > 0.9999:

Code: Select all

if (tau1/tau2 < 1e-9) {
  tau1 = tau2*1e-9
}