iaf_cond_exp_sfa_rr

class brainpy.state.iaf_cond_exp_sfa_rr(in_size, E_L=Quantity(-70., "mV"), C_m=Quantity(289.5, "pF"), t_ref=Quantity(0.5, "ms"), V_th=Quantity(-57., "mV"), V_reset=Quantity(-70., "mV"), E_ex=Quantity(0., "mV"), E_in=Quantity(-75., "mV"), g_L=Quantity(28.95, "nS"), tau_syn_ex=Quantity(1.5, "ms"), tau_syn_in=Quantity(10., "ms"), tau_sfa=Quantity(110., "ms"), tau_rr=Quantity(1.97, "ms"), E_sfa=Quantity(-70., "mV"), E_rr=Quantity(-70., "mV"), q_sfa=Quantity(14.48, "nS"), q_rr=Quantity(3214., "nS"), I_e=Quantity(0., "pA"), gsl_error_tol=1e-06, V_initializer=Constant(value=-70. mV), g_ex_initializer=Constant(value=0. nS), g_in_initializer=Constant(value=0. nS), g_sfa_initializer=Constant(value=0. nS), g_rr_initializer=Constant(value=0. nS), spk_fun=ReluGrad(alpha=0.3, width=1.0), spk_reset='hard', ref_var=False, name=None)

NEST-compatible conductance-based LIF neuron with spike-frequency adaptation and relative refractory mechanisms.

This model implements a conductance-based leaky integrate-and-fire neuron with exponential synaptic conductances, spike-frequency adaptation (SFA), and a relative refractory (RR) conductance mechanism. It follows the NEST iaf_cond_exp_sfa_rr model dynamics and update ordering exactly.

Mathematical Description

The model evolves five state variables:

1. Synaptic conductances (exponential decay):

   \[\frac{dg_{\mathrm{ex}}}{dt} = -\frac{g_{\mathrm{ex}}}{\tau_{\mathrm{syn,ex}}}, \qquad \frac{dg_{\mathrm{in}}}{dt} = -\frac{g_{\mathrm{in}}}{\tau_{\mathrm{syn,in}}}\]

2. Adaptation and relative refractory conductances (exponential decay):

   \[\frac{dg_{\mathrm{sfa}}}{dt} = -\frac{g_{\mathrm{sfa}}}{\tau_{\mathrm{sfa}}}, \qquad \frac{dg_{\mathrm{rr}}}{dt} = -\frac{g_{\mathrm{rr}}}{\tau_{\mathrm{rr}}}\]

3. Membrane potential:

   \[\frac{dV}{dt} = \frac{-I_{\mathrm{L}} + I_e + I_{\mathrm{stim}} - I_{\mathrm{syn,ex}} - I_{\mathrm{syn,in}} - I_{\mathrm{sfa}} - I_{\mathrm{rr}}}{C_m}\]

   where the individual currents are computed as:

   \[\begin{split}\begin{aligned} I_{\mathrm{L}} &= g_{\mathrm{L}} (V_{\mathrm{eff}} - E_{\mathrm{L}}) \\ I_{\mathrm{syn,ex}} &= g_{\mathrm{ex}} (V_{\mathrm{eff}} - E_{\mathrm{ex}}) \\ I_{\mathrm{syn,in}} &= g_{\mathrm{in}} (V_{\mathrm{eff}} - E_{\mathrm{in}}) \\ I_{\mathrm{sfa}} &= g_{\mathrm{sfa}} (V_{\mathrm{eff}} - E_{\mathrm{sfa}}) \\ I_{\mathrm{rr}} &= g_{\mathrm{rr}} (V_{\mathrm{eff}} - E_{\mathrm{rr}}) \end{aligned}\end{split}\]

   The effective voltage \(V_{\mathrm{eff}}\) implements NEST voltage clamping:

   • During absolute refractory period: \(V_{\mathrm{eff}} = V_{\mathrm{reset}}\)

   • Otherwise: \(V_{\mathrm{eff}} = \min(V, V_{\mathrm{th}})\)

   During absolute refractory period, \(dV/dt = 0\) while all conductances continue decaying.

Spike Dynamics

When \(V \geq V_{\mathrm{th}}\) and the neuron is not refractory:

1. A spike is emitted

2. \(V \leftarrow V_{\mathrm{reset}}\)

3. Absolute refractory period begins (duration \(t_{\mathrm{ref}}\))

4. Adaptation and RR conductances are incremented:

   \[g_{\mathrm{sfa}} \leftarrow g_{\mathrm{sfa}} + q_{\mathrm{sfa}}, \qquad g_{\mathrm{rr}} \leftarrow g_{\mathrm{rr}} + q_{\mathrm{rr}}\]

Numerical Integration

The ODEs are integrated using adaptive Runge-Kutta-Fehlberg 4(5) (RKF45) with absolute error tolerance gsl_error_tol. Each neuron maintains its own adaptive time step size (stored in integration_step), which is adjusted based on local error estimates. The minimum step size is _MIN_H = 1e-8 ms and maximum iterations per simulation step is _MAX_ITERS = 100000.

Update Ordering (NEST Semantics)

Per simulation step at time t:

1. Integrate ODEs over \((t, t+dt]\) using RKF45

2. Apply spike inputs: Add incoming delta inputs to g_ex and g_in

3. Refractory countdown: Decrement refractory counter if neuron is refractory

4. Threshold test: If \(V \geq V_{\mathrm{th}}\) and not refractory, emit spike

5. Reset and adaptation: On spike, reset voltage and increment g_sfa and g_rr

6. Buffer current: Store current input x into I_stim (one-step delay)

The one-step delayed current input mirrors NEST’s ring-buffer semantics.

Biological Interpretation

• g_sfa: Models spike-frequency adaptation through a slow potassium-like current that accumulates with repeated spiking and gradually decays. This causes firing rate to decrease during sustained input.

• g_rr: Models relative refractoriness through a transient hyperpolarizing current that makes the neuron harder to excite immediately after a spike, beyond the absolute refractory period. This provides a smooth transition back to normal excitability.

Parameters:

• in_size (int, tuple of int) – Population shape. Can be an integer for 1D population or tuple for multi-dimensional.

• E_L (ArrayLike, default: -70 mV) – Leak reversal potential. Must be less than V_th.

• C_m (ArrayLike, default: 289.5 pF) – Membrane capacitance. Must be strictly positive.

• t_ref (ArrayLike, default: 0.5 ms) – Absolute refractory period duration. Must be non-negative. During this period, voltage is clamped to V_reset and no spikes can occur.

• V_th (ArrayLike, default: -57 mV) – Spike threshold potential. Must be greater than V_reset.

• V_reset (ArrayLike, default: -70 mV) – Reset potential after spike. Must be less than V_th.

• E_ex (ArrayLike, default: 0 mV) – Excitatory synaptic reversal potential. Typically set to 0 mV (depolarizing).

• E_in (ArrayLike, default: -75 mV) – Inhibitory synaptic reversal potential. Typically set below E_L (hyperpolarizing).

• g_L (ArrayLike, default: 28.95 nS) – Leak conductance. Determines membrane time constant \(\tau_m = C_m / g_L\).

• tau_syn_ex (ArrayLike, default: 1.5 ms) – Excitatory synaptic conductance decay time constant. Must be strictly positive. Fast excitatory synapses (AMPA-like).

• tau_syn_in (ArrayLike, default: 10.0 ms) – Inhibitory synaptic conductance decay time constant. Must be strictly positive. Slower inhibitory synapses (GABA-A-like).

• tau_sfa (ArrayLike, default: 110.0 ms) – Spike-frequency adaptation conductance decay time constant. Must be strictly positive. Long timescale for slow adaptation (calcium-activated potassium currents).

• tau_rr (ArrayLike, default: 1.97 ms) – Relative refractory conductance decay time constant. Must be strictly positive. Short timescale for post-spike transient refractoriness.

• E_sfa (ArrayLike, default: -70 mV) – Adaptation reversal potential. Typically set to or below E_L for hyperpolarizing effect.

• E_rr (ArrayLike, default: -70 mV) – Relative refractory reversal potential. Typically set to or below E_L for hyperpolarizing effect.

• q_sfa (ArrayLike, default: 14.48 nS) – Spike-triggered adaptation conductance increment. Added to g_sfa on each spike. Controls adaptation strength.

• q_rr (ArrayLike, default: 3214.0 nS) – Spike-triggered relative refractory conductance increment. Added to g_rr on each spike. Controls relative refractoriness strength. Large value creates strong transient hyperpolarization.

• I_e (ArrayLike, default: 0 pA) – Constant external current injection. Positive values are depolarizing.

• gsl_error_tol (ArrayLike) – Unitless local RKF45 error tolerance, broadcastable and strictly positive.

• V_initializer (Callable, default: Constant(-70 mV)) – Initializer for membrane potential. Called as V_initializer(shape).

• g_ex_initializer (Callable, default: Constant(0 nS)) – Initializer for excitatory conductance. Called as g_ex_initializer(shape).

• g_in_initializer (Callable, default: Constant(0 nS)) – Initializer for inhibitory conductance. Called as g_in_initializer(shape).

• g_sfa_initializer (Callable, default: Constant(0 nS)) – Initializer for adaptation conductance. Called as g_sfa_initializer(shape).

• g_rr_initializer (Callable, default: Constant(0 nS)) – Initializer for relative refractory conductance. Called as g_rr_initializer(shape).

• spk_fun (Callable, default: ReluGrad()) – Surrogate gradient function for differentiable spike generation. Maps scaled voltage difference \((V - V_{\mathrm{th}}) / (V_{\mathrm{th}} - V_{\mathrm{reset}})\) to spike probability in \([0, 1]\).

• spk_reset (str, default: 'hard') – Spike reset mode. 'hard' uses stop_gradient (matches NEST), 'soft' allows gradient flow through reset.

• ref_var (bool, default: False) – If True, expose boolean state variable refractory indicating whether neuron is in absolute refractory period.

• name (str, optional) – Name of the neuron group. If None, auto-generated.

Parameter Mapping

Parameter	Default	Math equivalent
in_size	(required)	—
E_L	-70 mV	\(E_\mathrm{L}\)
C_m	289.5 pF	\(C_\mathrm{m}\)
t_ref	0.5 ms	\(t_\mathrm{ref}\)
V_th	-57 mV	\(V_\mathrm{th}\)
V_reset	-70 mV	\(V_\mathrm{reset}\)
E_ex	0 mV	\(E_\mathrm{ex}\)
E_in	-75 mV	\(E_\mathrm{in}\)
g_L	28.95 nS	\(g_\mathrm{L}\)
tau_syn_ex	1.5 ms	\(\tau_{\mathrm{syn,ex}}\)
tau_syn_in	10.0 ms	\(\tau_{\mathrm{syn,in}}\)
tau_sfa	110.0 ms	\(\tau_{\mathrm{sfa}}\)
tau_rr	1.97 ms	\(\tau_{\mathrm{rr}}\)
E_sfa	-70 mV	\(E_\mathrm{sfa}\)
E_rr	-70 mV	\(E_\mathrm{rr}\)
q_sfa	14.48 nS	\(q_\mathrm{sfa}\)
q_rr	3214.0 nS	\(q_\mathrm{rr}\)
I_e	0 pA	\(I_\mathrm{e}\)

State Variables

• V: HiddenState (float, shape in_size) — Membrane potential in mV

• g_ex: HiddenState (float, shape in_size) — Excitatory conductance in nS

• g_in: HiddenState (float, shape in_size) — Inhibitory conductance in nS

• g_sfa: HiddenState (float, shape in_size) — Adaptation conductance in nS

• g_rr: HiddenState (float, shape in_size) — Relative refractory conductance in nS

• refractory_step_count: ShortTermState (int32, shape in_size) — Remaining refractory steps

• integration_step: ShortTermState (float, shape in_size) — Adaptive RKF45 step size in ms

• I_stim: ShortTermState (float, shape in_size) — One-step delayed current buffer in pA

• last_spike_time: ShortTermState (float, shape in_size) — Last spike time in ms

• refractory: ShortTermState (bool, shape in_size) — Boolean refractory indicator (if ref_var=True)

Raises:

• ValueError – If V_reset >= V_th (reset must be below threshold)

• ValueError – If C_m <= 0 (capacitance must be positive)

• ValueError – If t_ref < 0 (refractory period cannot be negative)

• ValueError – If any time constant (tau_syn_ex, tau_syn_in, tau_sfa, tau_rr) is non-positive

See also

iaf_cond_exp

Simpler conductance-based LIF without adaptation or relative refractoriness

iaf_cond_alpha

Conductance-based LIF with alpha-function synaptic conductances

aeif_cond_exp

Exponential integrate-and-fire with conductance-based synapses

References

[1]

Meffin H, Burkitt AN, Grayden DB (2004). An analytical model for the large, fluctuating synaptic conductance state typical of neocortical neurons in vivo. Journal of Computational Neuroscience, 16:159-175. DOI: https://doi.org/10.1023/B:JCNS.0000014108.03012.81

[2]

Dayan P, Abbott LF (2001). Theoretical Neuroscience: Computational and Mathematical Modeling of Neural Systems. MIT Press.

[3]

NEST Simulator. iaf_cond_exp_sfa_rr documentation. https://nest-simulator.readthedocs.io/

Examples

Create a single neuron with default parameters:

>>> import brainstate as bst
>>> import saiunit as u
>>> import brainpy.state as bp
>>> neuron = bp.iaf_cond_exp_sfa_rr(in_size=1)
>>> with bst.environ.context(dt=0.1 * u.ms):
...     neuron.init_all_states()
...     print(neuron.V.value)
[-70.] mV

Simulate a population with constant current injection:

>>> import matplotlib.pyplot as plt
>>> neuron = bp.iaf_cond_exp_sfa_rr(in_size=10, I_e=500 * u.pA)
>>> with bst.environ.context(dt=0.1 * u.ms):
...     neuron.init_all_states()
...     voltages = []
...     for _ in range(1000):
...         spike = neuron.update()
...         voltages.append(neuron.V.value[0])
>>> # plt.plot(voltages)  # Shows adaptation: decreasing firing rate

Demonstrate spike-frequency adaptation:

>>> # Strong adaptation (large q_sfa)
>>> neuron_adapt = bp.iaf_cond_exp_sfa_rr(in_size=1, I_e=600*u.pA, q_sfa=50*u.nS)
>>> # Weak adaptation (small q_sfa)
>>> neuron_weak = bp.iaf_cond_exp_sfa_rr(in_size=1, I_e=600*u.pA, q_sfa=5*u.nS)
>>> # neuron_adapt will show stronger decrease in firing rate over time

Notes

• Computational cost: This model uses adaptive RKF45 integration, which is more expensive than fixed-step exponential Euler used in simpler models like iaf_cond_exp. However, it provides better accuracy for stiff dynamics.

• NEST compatibility: This implementation exactly reproduces NEST behavior including voltage clamping, update ordering, and one-step delayed current semantics.

• Gradient flow: The 'hard' reset mode (default) uses stop_gradient on reset, which is necessary for NEST compatibility but prevents gradient flow through spike reset. Use spk_reset='soft' for better gradient-based learning, at the cost of deviating from NEST semantics.

• Parameter tuning: The default q_sfa and q_rr values are taken from NEST defaults and produce moderate adaptation. Increase q_sfa for stronger adaptation (more pronounced firing rate decrease). Increase q_rr for stronger post-spike hyperpolarization.

__init__(in_size, E_L=Quantity(-70., "mV"), C_m=Quantity(289.5, "pF"), t_ref=Quantity(0.5, "ms"), V_th=Quantity(-57., "mV"), V_reset=Quantity(-70., "mV"), E_ex=Quantity(0., "mV"), E_in=Quantity(-75., "mV"), g_L=Quantity(28.95, "nS"), tau_syn_ex=Quantity(1.5, "ms"), tau_syn_in=Quantity(10., "ms"), tau_sfa=Quantity(110., "ms"), tau_rr=Quantity(1.97, "ms"), E_sfa=Quantity(-70., "mV"), E_rr=Quantity(-70., "mV"), q_sfa=Quantity(14.48, "nS"), q_rr=Quantity(3214., "nS"), I_e=Quantity(0., "pA"), gsl_error_tol=1e-06, V_initializer=Constant(value=-70. mV), g_ex_initializer=Constant(value=0. nS), g_in_initializer=Constant(value=0. nS), g_sfa_initializer=Constant(value=0. nS), g_rr_initializer=Constant(value=0. nS), spk_fun=ReluGrad(alpha=0.3, width=1.0), spk_reset='hard', ref_var=False, name=None)

Initialize the iaf_cond_exp_sfa_rr neuron model.

All parameters are validated to ensure physical consistency. Parameters can be scalars (broadcast to all neurons) or arrays matching in_size.

get_spike(V=None)

Compute differentiable spike output using surrogate gradient.

Applies the surrogate gradient function to the scaled voltage difference. The voltage is scaled to \([0, 1]\) range where 0 corresponds to V_reset and 1 corresponds to V_th. Values above 1 (threshold crossing) produce spike output near 1.

Parameters:

V (ArrayLike, optional) – Membrane potential in mV. If None, uses current self.V.value.

Returns:

Differentiable spike indicator in [0, 1], shape matching V. Values near 1 indicate spike, near 0 indicate no spike.

Return type:

ArrayLike

Notes

• The surrogate function is specified by spk_fun parameter

• The scaling ensures consistent behavior across different voltage ranges

• During backpropagation, gradients flow through the surrogate function

init_state(**kwargs)

Initialize persistent and short-term state variables.

Parameters:

**kwargs – Unused compatibility parameters accepted by the base-state API.

Raises:

• ValueError – If an initializer cannot be broadcast to requested shape.

• TypeError – If initializer outputs have incompatible units/dtypes for the corresponding state variables.

update(x=Quantity(0., 'pA'))

Advance the neuron state by one simulation time step.

Implements the complete NEST update cycle:

1. ODE Integration: Integrate voltage and conductances over [t, t+dt] using RKF45

2. Spike Input: Apply delta inputs from incoming spikes to g_ex and g_in

3. Refractory Logic: Handle absolute refractory countdown and voltage clamping

4. Threshold Test: Detect threshold crossing and emit spike if not refractory

5. Reset and Adaptation: On spike, reset voltage and increment g_sfa and g_rr

6. Current Buffering: Store current input x for next time step (one-step delay)

Parameters:

x (ArrayLike, default: 0 pA) – External current input for the next time step in pA, shape matching in_size or broadcastable. This input is buffered and applied with one-step delay, mirroring NEST ring-buffer semantics.

Returns:

Binary spike tensor with dtype jnp.float64 and shape self.V.value.shape. A value of 1.0 indicates at least one internal spike event occurred during the integrated interval \((t, t+dt]\).

Return type:

jax.Array

Raises:

ValueError – If RKF45 integration enters a guarded unstable regime (V < -1e3 mV), indicating divergent dynamics for the current parameter/input regime.

Notes

Integration is performed with an adaptive vectorized RKF45 loop, including in-loop spike/reset/adaptation events and optional multiple spikes per step. All arithmetic is unit-aware via saiunit.math.