Add point-neuron and current-source support to the Arbor backend

pyNN/arbor/nmodl/adexp.mod

@@ -0,0 +1,79 @@
+: Adaptive-exponential (Brette-Gerstner) dynamics for building AdExp / EIF point
+: neurons as Arbor cable cells. Like lif.mod this is a POINT_PROCESS placed at the
+: soma whose reset is driven by the cell's threshold_detector via POST_EVENT; it
+: additionally contributes the exponential spike-generating current and the
+: adaptation current w. The sub-threshold leak -GL*(v - EL) is a separate `pas`
+: density mechanism (with the same GL), exactly as for lif.mod.
+:
+: Membrane dynamics realised (with C dv/dt = -sum of mechanism currents):
+:     C dv/dt = -GL(v-EL) + GL*delta*exp((v-vthresh)/delta) - w + I
+:     tau_w dw/dt = a*(v - EL) - w
+: so this mechanism contributes i = GL*delta term (as iexp, negative/inward) + w.
+: On a spike (detector crosses v_spike) POST_EVENT starts a refractory countdown and
+: increments w by b; while refractory a strong clamp holds v at v_reset (and the
+: exponential/adaptation currents are gated off), matching lif.mod. The clamp is
+: gated by a voltage-independent 0/1 multiplier so modcc derives the conductance
+: correctly (see lif.mod).
+NEURON {
+    POINT_PROCESS adexp
+    RANGE v_reset, t_ref, g_reset, GL, delta, vthresh, a, b, tau_w, EL
+    NONSPECIFIC_CURRENT i
+}
+
+UNITS {
+    (mV) = (millivolt)
+    (nA) = (nanoamp)
+    (uS) = (microsiemens)
+    (ms) = (millisecond)
+}
+
+PARAMETER {
+    v_reset = -70.6  (mV)
+    t_ref   = 0.1    (ms)
+    g_reset = 1000   (uS)
+    GL      = 0.03   (uS)   : leak conductance (= C_m/tau_m; matches the pas leak)
+    delta   = 2      (mV)   : steepness of the exponential
+    vthresh = -50.4  (mV)   : exponential (soft) threshold V_T
+    a       = 0.004  (uS)   : subthreshold adaptation conductance
+    b       = 0.0805 (nA)   : spike-triggered adaptation increment
+    tau_w   = 144    (ms)   : adaptation time constant
+    EL      = -70.6  (mV)   : leak reversal (= v_rest)
+}
+
+STATE {
+    refrac (ms)
+    w      (nA)
+}
+
+ASSIGNED { clamp  iexp (nA) }
+
+INITIAL {
+    refrac = 0
+    w = 0
+    clamp = 0
+}
+
+BREAKPOINT {
+    SOLVE states METHOD cnexp
+    clamp = 0
+    if (refrac > 0) { clamp = 1 }
+    iexp = exp_current(v)
+    i = clamp * g_reset * (v - v_reset) + (1 - clamp) * (iexp + w)
+}
+
+DERIVATIVE states {
+    refrac' = -1
+    w' = (a * (v - EL) - w) / tau_w
+}
+
+FUNCTION exp_current(v (mV)) (nA) {
+    LOCAL arg
+    arg = (v - vthresh) / delta
+    if (arg > 50) { arg = 50 }   : guard against overflow before the reset clamps v
+    exp_current = -GL * delta * exp(arg)
+}
+
+POST_EVENT(time) {
+    refrac = t_ref
+    w = w + b
+}

pyNN/arbor/nmodl/alphasyn.mod

@@ -0,0 +1,55 @@
+: Conductance-based synapse with an alpha-function time course, for IF_cond_alpha
+: and conductance-based point neurons. An incoming event of weight `w` (uS) produces
+: a conductance g(t) = w * (t/tau) * exp(1 - t/tau), peaking at `w` a time `tau`
+: after the event; the membrane current is i = g*(v - e).
+:
+: The alpha shape is the impulse response of the coupled linear system
+:     g' = z - g/tau
+:     z' = -z/tau
+: with each event adding w*exp(1)/tau to z (so that the peak of g equals w). Unlike
+: the NEURON backend's alphasyn.mod this needs no spike queue. The system is a
+: non-diagonalisable Jordan block, so `cnexp` cannot solve it; `sparse` (Arbor's
+: implicit/backward-Euler solver) is used instead, which is first-order accurate --
+: the alpha peak is slightly under-estimated at large dt (a few % of the synaptic
+: current for dt ~ tau), but the downstream membrane response matches the exact
+: alpha to <1% for dt <= 0.05 ms.
+NEURON {
+    POINT_PROCESS alphasyn
+    RANGE tau, e
+    NONSPECIFIC_CURRENT i
+}
+
+UNITS {
+    (mV) = (millivolt)
+    (uS) = (microsiemens)
+    (ms) = (millisecond)
+}
+
+PARAMETER {
+    tau = 5.0 (ms)
+    e   = 0   (mV)
+}
+
+STATE {
+    g (uS)
+    z (uS/ms)
+}
+
+INITIAL {
+    g = 0
+    z = 0
+}
+
+BREAKPOINT {
+    SOLVE state METHOD sparse
+    i = g * (v - e)
+}
+
+DERIVATIVE state {
+    g' = z - g/tau
+    z' = -z/tau
+}
+
+NET_RECEIVE(weight (uS)) {
+    z = z + weight * exp(1) / tau
+}

pyNN/arbor/nmodl/alphasyn_curr.mod

@@ -0,0 +1,46 @@
+: Current-based synapse with an alpha-function time course, for IF_curr_alpha and
+: current-based point neurons. An incoming event of weight `w` (nA) produces a
+: synaptic current isyn(t) = w * (t/tau) * exp(1 - t/tau), peaking at `w` a time
+: `tau` after the event. A positive weight injects depolarising (inward) current,
+: matching PyNN's convention, so the contributed membrane current is i = -isyn.
+:
+: Same coupled-ODE realisation of the alpha shape as alphasyn.mod (no spike queue);
+: solved with `sparse` (implicit) because the Jordan-block system is not diagonal.
+NEURON {
+    POINT_PROCESS alphasyn_curr
+    RANGE tau
+    NONSPECIFIC_CURRENT i
+}
+
+UNITS {
+    (nA) = (nanoamp)
+    (ms) = (millisecond)
+}
+
+PARAMETER {
+    tau = 5.0 (ms)
+}
+
+STATE {
+    isyn (nA)
+    z    (nA/ms)
+}
+
+INITIAL {
+    isyn = 0
+    z = 0
+}
+
+BREAKPOINT {
+    SOLVE state METHOD sparse
+    i = -isyn
+}
+
+DERIVATIVE state {
+    isyn' = z - isyn/tau
+    z' = -z/tau
+}
+
+NET_RECEIVE(weight (nA)) {
+    z = z + weight * exp(1) / tau
+}

pyNN/arbor/nmodl/expsyn_curr.mod

@@ -0,0 +1,40 @@
+: Current-based synapse with exponential decay, for IF_curr_exp / current-based
+: point neurons. An incoming event steps the synaptic current by `weight` (nA),
+: which then decays exponentially with time constant `tau`. A positive weight
+: injects depolarising (inward) current, matching PyNN's isyn convention, so the
+: contributed membrane current is i = -isyn.
+NEURON {
+    POINT_PROCESS expsyn_curr
+    RANGE tau
+    NONSPECIFIC_CURRENT i
+}
+
+UNITS {
+    (nA) = (nanoamp)
+    (ms) = (millisecond)
+}
+
+PARAMETER {
+    tau = 5.0 (ms)
+}
+
+STATE {
+    isyn (nA)
+}
+
+INITIAL {
+    isyn = 0
+}
+
+BREAKPOINT {
+    SOLVE state METHOD cnexp
+    i = -isyn
+}
+
+DERIVATIVE state {
+    isyn' = -isyn/tau
+}
+
+NET_RECEIVE(weight) {
+    isyn = isyn + weight
+}

pyNN/arbor/nmodl/hh_traub.mod

@@ -0,0 +1,64 @@
+: Traub-modified Hodgkin-Huxley channels, for building HH_cond_exp as an Arbor
+: single-compartment cable cell. A near-verbatim port of the NEURON backend's
+: hh_traub.mod: the TABLE statement is dropped (Arbor's modcc has no TABLE; the
+: rates are computed exactly each step) and the rate "trap" is written with Arbor's
+: exprelr (vtrap(x,y) = y*exprelr(x/y) = x/(exp(x/y)-1), singularity-safe). The
+: NEST-compatible initial state m = h = n = 0 is kept. Reversal potentials for Na/K
+: come from the ions (set per cell); the leak is a NONSPECIFIC_CURRENT.
+NEURON {
+    SUFFIX hh_traub
+    USEION na READ ena WRITE ina
+    USEION k READ ek WRITE ik
+    NONSPECIFIC_CURRENT il
+    RANGE gnabar, gkbar, gl, el, vT
+}
+
+UNITS {
+    (mV) = (millivolt)
+    (S)  = (siemens)
+}
+
+PARAMETER {
+    gnabar = 0.02    (S/cm2)
+    gkbar  = 0.006   (S/cm2)
+    gl     = 0.00001 (S/cm2)
+    el     = -60.0   (mV)
+    vT     = -63.0   (mV)
+}
+
+STATE { m h n }
+
+BREAKPOINT {
+    SOLVE states METHOD cnexp
+    ina = gnabar * m * m * m * h * (v - ena)
+    ik  = gkbar * n * n * n * n * (v - ek)
+    il  = gl * (v - el)
+}
+
+INITIAL {
+    : for compatibility with NEST, start from m = h = n = 0
+    m = 0
+    h = 0
+    n = 0
+}
+
+DERIVATIVE states {
+    LOCAL u, alpha, beta
+    u = v - vT
+    : sodium activation
+    alpha = 0.32 * vtrap(13 - u, 4)
+    beta  = 0.28 * vtrap(u - 40, 5)
+    m' = alpha - m * (alpha + beta)
+    : sodium inactivation
+    alpha = 0.128 * exp((17 - u) / 18)
+    beta  = 4 / (exp((40 - u) / 5) + 1)
+    h' = alpha - h * (alpha + beta)
+    : potassium activation
+    alpha = 0.032 * vtrap(15 - u, 5)
+    beta  = 0.5 * exp((10 - u) / 40)
+    n' = alpha - n * (alpha + beta)
+}
+
+FUNCTION vtrap(x, y) {
+    vtrap = y * exprelr(x / y)
+}

pyNN/arbor/nmodl/izhikevich.mod

@@ -0,0 +1,75 @@
+: Izhikevich (2003) quadratic integrate-and-fire neuron, for building the PyNN
+: Izhikevich cell type as an Arbor cable cell:
+:     dv/dt = 0.04 v^2 + 5 v + 140 - u + I
+:     du/dt = a (b v - u)
+: with reset, when v reaches vthresh: v -> c, u -> u + d.
+:
+: v is integrated by Arbor's cable solver: with the cell's absolute capacitance set
+: equal to Cm, this POINT_PROCESS supplies i = -Cm(0.04 v^2 + 5 v + 140 - u) so that
+: Cm dv/dt = -i reproduces the Izhikevich dv/dt (the injected current I enters via a
+: separate iclamp for i_offset, contributing I/Cm as in the NEURON backend). There
+: is no leak mechanism. u is a STATE solved by cnexp.
+:
+: Arbor mechanisms cannot write v, so the reset (v -> c) is done as in lif.mod: the
+: threshold_detector (set to vthresh) delivers a POST_EVENT that starts a very short
+: countdown during which a strong clamp holds v at c; the countdown t_reset is tiny
+: (<= one timestep) so there is no artificial refractory period, only a one-step
+: reset. u += d is applied in the POST_EVENT. The clamp is gated by the
+: voltage-independent 0/1 multiplier `clamp` so modcc derives the conductance
+: correctly (see lif.mod).
+NEURON {
+    POINT_PROCESS izhikevich
+    RANGE a, b, c, d, vthresh, Cm, uinit, t_reset, g_reset
+    NONSPECIFIC_CURRENT i
+}
+
+UNITS {
+    (mV) = (millivolt)
+    (nA) = (nanoamp)
+    (nF) = (nanofarad)
+    (uS) = (microsiemens)
+    (ms) = (millisecond)
+}
+
+PARAMETER {
+    a       = 0.02   (/ms)
+    b       = 0.2    (/ms)
+    c       = -65    (mV)      : reset potential
+    d       = 2      (mV/ms)   : reset increment of u
+    vthresh = 30     (mV)      : spike / reset threshold
+    Cm      = 0.001  (nF)      : capacitance (matches the cell's absolute cm)
+    uinit   = -14    (mV/ms)
+    t_reset = 0.001  (ms)      : clamp duration (<= dt: a one-step reset)
+    g_reset = 1000   (uS)
+}
+
+STATE {
+    u      (mV/ms)
+    refrac (ms)
+}
+
+ASSIGNED { clamp }
+
+INITIAL {
+    u = uinit
+    refrac = 0
+    clamp = 0
+}
+
+BREAKPOINT {
+    SOLVE states METHOD cnexp
+    clamp = 0
+    if (refrac > 0) { clamp = 1 }
+    i = clamp * g_reset * (v - c)
+        + (1 - clamp) * (-Cm * (0.04 * v * v + 5 * v + 140 - u))
+}
+
+DERIVATIVE states {
+    u' = a * (b * v - u)
+    refrac' = -1
+}
+
+POST_EVENT(time) {
+    refrac = t_reset
+    u = u + d
+}