Alpha-Synapse evaluation time in Backward Euler scheme

[danarwed]

Dear Neuron-Experts,

I have a simple passive model (see the attached code below) consisting of a uniform cylinder with a alpha synapse located in the middle of the cyliner (therefore using an odd number of nseg).
The synapse is activated after 1ms with a duration of 1 ms. I try to solve the model with the Backward Euler (BE) method with dt=0.005ms and as I wanted to check the outcome of the alpha-function in the synapse
and the induced synaptic current I put a line like

Code: Select all

printf("t: %.5g  tau: %.3g gmax: %.3g  alpha: %.16g  v: %.16g\n", t, tau, gmax, alpha(_threadargscomma_ (t-onset)/ tau), v);

in the function _nrn_current (in syn,c) which I identified to handle the synapse activation (I am sure there are better ways using the high end user functionality but I am still a very beginner of Neuron).
What I do not understand is the output, I'll place some lines around t = 1ms here:

Code: Select all

0.99 0.002827275848336286
t: 0.9925  tau: 1 gmax: 10  alpha: 0  v: 0.003827275848336286
t: 0.9925  tau: 1 gmax: 10  alpha: 0  v: 0.002827275848336286
0.995 0.002823746165629249
t: 0.9975  tau: 1 gmax: 10  alpha: 0  v: 0.003823746165629249
t: 0.9975  tau: 1 gmax: 10  alpha: 0  v: 0.002823746165629249
1 0.002820220889517352
t: 1.0025  tau: 1 gmax: 10  alpha: 0.006778736528582666  v: 0.003820220889517352
t: 1.0025  tau: 1 gmax: 10  alpha: 0.006778736528582666  v: 0.002820220889517352
1.005 0.00490483032996971
t: 1.0075  tau: 1 gmax: 10  alpha: 0.02023478231735222  v: 0.4240529365368073
t: 1.0075  tau: 1 gmax: 10  alpha: 0.02023478231735222  v: 0.4230529365368073
1.01 0.01304325012706313
t: 1.0125  tau: 1 gmax: 10  alpha: 0.0335564348658731  v: 1.644988063163831
t: 1.0125  tau: 1 gmax: 10  alpha: 0.0335564348658731  v: 1.643988063163831

My question is:
It seems to me the synapse is evaluated at t = t_old + 0.5 * dt. Is this really true? I thought that for BE, one has to solve (1-dt*A) v_new = v_old + dt * f(v_new) where A is the spatial discretization matrix,
v_new the potential at t_new=t+dt and v_old at t=t_old. f(v+new) is the synaptical current at t=t_new, that is i(v_new) = g_syn(t_new) [v_new-erev]. So the synapse should be evaluated at t_new = t +dt but not
at t + 0.5 *dt. Hopefully I am right in my assumptions and haven't done a stupid mistake in my considerations.

Thank you for your help, it is very much appreciated.
Dan

P.S. Here the model to be solved:

Code: Select all

secondorder = 0   // Backward Euler, default
ndend=1              // number of dendrites
nsegments=1      // Nr. of segments, should be odd

create soma
create dend[ndend]

/* specify anatomical and biophysical properties */
soma {
    nseg = nsegments    // number of segments
    diam(0:1) = 4:4         // [µm], dendritic diameter tapers along its length
    diam = 4                  // [µm], dendritic diameter tapers along its length
    L = 158.113885        // [µm], length
    Ra = 200.                // [Ohm*cm], axial resistivity, default: 35.4
    // passive current
    insert pas
    g_pas = 0.0005      // [S/cm^2], conductive
    e_pas = 0              // [mv], equilibrium potential
    cm = 2.                 // [uf/cm^2]
}

dend[0] {
    nseg = nsegments    // number of segments
    diam(0:1) = 4:4         // [µm], dendritic diameter tapers along its length
    diam = 4                  // [µm], dendritic diameter tapers along its length
    L    = 158.113885     // [µm], length
    Ra   = 200.              // [Ohm*cm], axial resistivity, default: 35.4
    // passive current
    insert pas
    g_pas = 0.0005       // [S/cm^2], conductive
    e_pas = 0               // [mv], equilibrium potential
    cm    = 2.               // [uf/cm^2]
}

/* connect soma and dends, check topology */
connect dend[0](0), soma(1)
topology()

/* synapse */
objectvar syn
dend[0] syn  = new AlphaSynapse(0.5) // put in the middle of the dendrite
syn.tau        = 1.0                               // [ms] t_peak
syn.onset    = 1.0                               // [ms] t_start
syn.e          = 50                                // [mV] synaptical reversal potential
syn.gmax    = 10                                // [nS] maximal synaptical conductivity

/* initial condition and other model parameters */
dt    = 0.005                // [ms], time step for use with fadvance
tstop = 10                  // [ms], end time
finitialize(0.00362067) // [mv], init membrane potential

/* solve */
proc integrate() {
    print soma.nseg, dend[0].nseg
    printf("%g %.16g\n", t, soma.v(0.5))
    while(t < tstop) {
        fadvance()
        printf("%g %.16g\n", t, soma.v(0.5))
    }
}

[ted]

as I wanted to check the outcome of the alpha-function in the synapse
and the induced synaptic current I put a line like

Code: Select all

printf("t: %.5g  tau: %.3g gmax: %.3g  alpha: %.16g  v: %.16g\n", t, tau, gmax, alpha(_threadargscomma_ (t-onset)/ tau), v);

in the function _nrn_current (in syn,c) which I identified to handle the synapse activation

Bad idea. You sound like the fictional C programmer who was in a hurry to write a program in Python. Whenever he came to something that he didn't know how to do in Python, he just hacked the C source code so the Python interpreter would do what he wanted.

Most of what you need to do with NEURON--especially setting up models, controlling simulations, and analyzing results--can be done by code written in its interpreted languages hoc and/or Python. NEURON's compiled language NMODL is helpful for adding new functions (including new biophysical mechanisms) to NEURON (which are callable from hoc or Python), and occasionally for very special purposes it can be useful to insert C code into a VERBATIM block inside an NMODL file.

But for the most part you only need hoc and/or Python. Until you know the limits of what can be done with these tools, and understand much more about how NEURON is organized and how it works, hacking its source code is unlikely to be useful.

Take the time to read some documentation--one or two of the "Key papers about NEURON" listed on the WWW site's Documentation page http://www.neuron.yale.edu/neuron/docs, maybe chapter 5 and 6 in The NEURON Book. Examine some of the example programs that are included with NEURON, while looking up keywords in the Programmer's Reference. Browse through the FAQ list. Check out the Forum. Ask questions if you don't find the answers.

Here's the Programmer's Reference entry about the AlphaSynapse class.
http://www.neuron.yale.edu/neuron/stati ... phaSynapse
If syn is an objref for an AlphaSynapse, syn.i returns the current delivered by that synapse.

What I do not understand is the output

NEURON has a complex run time system that controls the sequence of computations that are performed in advancing the solution in time. These computations can involve multiple evaluations of internal functions that produce intermediate results; your printf is giving you a peek at some of these results. For details, read chapter 7 of The NEURON Book.

Stick with hoc or Python, which get the final results of each fadvance().