aeif_psc_alpha

aeif_psc_alpha#

class brainpy.state.aeif_psc_alpha(in_size, V_peak=Quantity(0., "mV"), V_reset=Quantity(-60., "mV"), t_ref=Quantity(0., "ms"), g_L=Quantity(30., "nS"), C_m=Quantity(281., "pF"), E_L=Quantity(-70.6, "mV"), Delta_T=Quantity(2., "mV"), tau_w=Quantity(144., "ms"), a=Quantity(4., "nS"), b=Quantity(80.5, "pA"), V_th=Quantity(-50.4, "mV"), tau_syn_ex=Quantity(0.2, "ms"), tau_syn_in=Quantity(2., "ms"), I_e=Quantity(0., "pA"), gsl_error_tol=1e-06, V_initializer=Constant(value=-70.6 mV), I_ex_initializer=Constant(value=0. pA), I_in_initializer=Constant(value=0. pA), w_initializer=Constant(value=0. pA), spk_fun=ReluGrad(alpha=0.3, width=1.0), spk_reset='hard', ref_var=False, name=None)#

NEST-compatible adaptive exponential integrate-and-fire neuron with alpha-shaped postsynaptic currents.

This model implements the adaptive exponential integrate-and-fire (AdEx) neuron with current-based synapses following alpha-function kinetics. It replicates the behavior of NEST’s aeif_psc_alpha model, including adaptive Runge-Kutta-Fehlberg (RKF45) numerical integration, in-loop spike detection and reset, and NEST-compatible refractory handling.

1. Mathematical Model

Membrane Dynamics

The subthreshold membrane potential \(V\) evolves according to:

\[C_m \frac{dV}{dt} = -g_L (V - E_L) + g_L \Delta_T \exp\left(\frac{V - V_{th}}{\Delta_T}\right) - w + I_{ex} - I_{in} + I_e + I_{stim}\]

where:

  • \(C_m\) – membrane capacitance

  • \(g_L\) – leak conductance

  • \(E_L\) – leak reversal potential

  • \(\Delta_T\) – exponential slope factor (spike sharpness)

  • \(V_{th}\) – spike initiation threshold

  • \(w\) – adaptation current

  • \(I_{ex}\), \(I_{in}\) – excitatory and inhibitory synaptic currents

  • \(I_e\) – constant external current

  • \(I_{stim}\) – time-varying external input (one-step delayed)

The exponential term \(g_L \Delta_T \exp((V - V_{th})/\Delta_T)\) causes rapid voltage acceleration near \(V_{th}\), producing spike initiation. Setting \(\Delta_T = 0\) recovers the leaky integrate-and-fire (LIF) limit.

Adaptation Dynamics

The adaptation current \(w\) provides spike-frequency adaptation and subthreshold coupling:

\[\tau_w \frac{dw}{dt} = a(V - E_L) - w\]

  • Subthreshold adaptation: parameter \(a\) couples \(w\) to membrane potential

  • Spike-triggered adaptation: at each spike, \(w \leftarrow w + b\)

Alpha-Function Synaptic Currents

Excitatory and inhibitory currents are modeled as alpha functions, each requiring two state variables:

\[\frac{d\,dI_{ex}}{dt} = -\frac{dI_{ex}}{\tau_{syn,ex}}, \quad \frac{dI_{ex}}{dt} = dI_{ex} - \frac{I_{ex}}{\tau_{syn,ex}}\]
\[\frac{d\,dI_{in}}{dt} = -\frac{dI_{in}}{\tau_{syn,in}}, \quad \frac{dI_{in}}{dt} = dI_{in} - \frac{I_{in}}{\tau_{syn,in}}\]

Incoming spike weights \(w_j\) (in pA) are split by sign and delivered as:

\[dI_{ex} \leftarrow dI_{ex} + \frac{e}{\tau_{syn,ex}} \max(w_j, 0)\]
\[dI_{in} \leftarrow dI_{in} + \frac{e}{\tau_{syn,in}} \max(-w_j, 0)\]

where \(e = \exp(1)\) provides the alpha-function normalization.

2. Spike Detection and Reset

Threshold Crossing

Spike detection threshold depends on \(\Delta_T\):

  • If \(\Delta_T > 0\): spike when \(V \geq V_{peak}\)

  • If \(\Delta_T = 0\): spike when \(V \geq V_{th}\) (LIF-like)

Reset Mechanism

Upon spike detection:

  1. \(V \leftarrow V_{reset}\)

  2. \(w \leftarrow w + b\) (spike-triggered adaptation)

  3. Refractory counter set to \(\lceil t_{ref}/dt \rceil + 1\) (if \(t_{ref} > 0\))

Spike detection and reset occur inside the RKF45 integration substeps, allowing multiple spikes per simulation time step when \(t_{ref} = 0\).

3. Refractory Period Handling

During the refractory period (\(r > 0\) steps remaining):

  • Membrane potential clamped: \(V_{eff} = V_{reset}\)

  • Voltage derivative forced: \(dV/dt = 0\)

  • Alpha currents and adaptation continue evolving normally

After each time step, the refractory counter is decremented: \(r \leftarrow r - 1\).

4. Numerical Integration

The model uses adaptive Runge-Kutta-Fehlberg (RKF45) with local error control:

  • Order: 5th-order accurate solution with 4th-order error estimate

  • Error tolerance: controlled by gsl_error_tol (default \(10^{-6}\))

  • Step size adaptation: \(h_{new} = h \cdot \min(5, \max(0.2, 0.9 (\epsilon/\text{err})^{0.2}))\)

  • Minimum step: \(h_{min} = 10^{-8}\) ms to prevent stalling

  • Persistent step size: each neuron maintains its own integration step size across time

The RKF45 Butcher tableau coefficients follow the standard formulation with stages \(k_1\) through \(k_6\).

5. Update Sequence

Each simulation step processes state updates in this order:

  1. Integration loop: Integrate ODEs from \(t\) to \(t + dt\) using RKF45 substeps, checking for spikes and applying resets within the loop

  2. Refractory decrement: After integration, decrement refractory counter once

  3. Synaptic input delivery: Add spike weights to \(dI_{ex}\) and \(dI_{in}\)

  4. External current update: Store current input \(x\) into one-step-delayed buffer \(I_{stim}\) (to be used in next step)

Parameters:

  • in_size (int, tuple of int) – Shape of the neuron population. Can be an integer (1D) or tuple (multi-dimensional).

  • V_peak (ArrayLike, optional) – Spike detection threshold voltage. Default: 0.0 * u.mV. Used for threshold detection when \(\Delta_T > 0\).

  • V_reset (ArrayLike, optional) – Reset potential after spike. Default: -60.0 * u.mV.

  • t_ref (ArrayLike, optional) – Absolute refractory period duration. Default: 0.0 * u.ms. During refractory period, \(V\) is clamped to \(V_{reset}\) and \(dV/dt = 0\).

  • g_L (ArrayLike, optional) – Leak conductance. Default: 30.0 * u.nS.

  • C_m (ArrayLike, optional) – Membrane capacitance. Default: 281.0 * u.pF. Determines membrane time constant \(\tau_m = C_m / g_L\).

  • E_L (ArrayLike, optional) – Leak reversal potential. Default: -70.6 * u.mV.

  • Delta_T (ArrayLike, optional) – Exponential slope factor. Default: 2.0 * u.mV. Controls spike sharpness; set to 0 for LIF-like behavior.

  • tau_w (ArrayLike, optional) – Adaptation time constant. Default: 144.0 * u.ms.

  • a (ArrayLike, optional) – Subthreshold adaptation coupling. Default: 4.0 * u.nS. Couples adaptation current to membrane potential deviation from \(E_L\).

  • b (ArrayLike, optional) – Spike-triggered adaptation increment. Default: 80.5 * u.pA. Added to \(w\) at each spike.

  • V_th (ArrayLike, optional) – Spike initiation threshold. Default: -50.4 * u.mV. Appears in exponential term and as fallback spike threshold when \(\Delta_T = 0\).

  • tau_syn_ex (ArrayLike, optional) – Excitatory synaptic alpha time constant. Default: 0.2 * u.ms.

  • tau_syn_in (ArrayLike, optional) – Inhibitory synaptic alpha time constant. Default: 2.0 * u.ms.

  • I_e (ArrayLike, optional) – Constant external current. Default: 0.0 * u.pA.

  • gsl_error_tol (ArrayLike, optional) – RKF45 local error tolerance. Default: 1e-6. Smaller values increase accuracy but require smaller integration steps.

  • V_initializer (Callable, optional) – Membrane potential initializer. Default: Constant(-70.6 * u.mV).

  • I_ex_initializer (Callable, optional) – Excitatory current initializer. Default: Constant(0.0 * u.pA).

  • I_in_initializer (Callable, optional) – Inhibitory current initializer. Default: Constant(0.0 * u.pA).

  • w_initializer (Callable, optional) – Adaptation current initializer. Default: Constant(0.0 * u.pA).

  • spk_fun (Callable, optional) – Surrogate gradient function for differentiable spike generation. Default: ReluGrad().

  • spk_reset (str, optional) – Spike reset mode: 'soft' (subtract threshold) or 'hard' (stop gradient). Default: 'hard' (matches NEST behavior).

  • ref_var (bool, optional) – If True, expose a boolean refractory state variable indicating refractory status. Default: False.

  • name (str, optional) – Name of the neuron population.

Parameter Mapping

Parameter

Default

Math equivalent

Description

in_size

(required)

—

Population shape

V_peak

0 mV

\(V_\mathrm{peak}\)

Spike detection threshold (if \(\Delta_T > 0\))

V_reset

-60 mV

\(V_\mathrm{reset}\)

Reset potential

t_ref

0 ms

\(t_\mathrm{ref}\)

Absolute refractory duration

g_L

30 nS

\(g_\mathrm{L}\)

Leak conductance

C_m

281 pF

\(C_\mathrm{m}\)

Membrane capacitance

E_L

-70.6 mV

\(E_\mathrm{L}\)

Leak reversal potential

Delta_T

2 mV

\(\Delta_T\)

Exponential slope factor

tau_w

144 ms

\(\tau_w\)

Adaptation time constant

a

4 nS

\(a\)

Subthreshold adaptation coupling

b

80.5 pA

\(b\)

Spike-triggered adaptation increment

V_th

-50.4 mV

\(V_\mathrm{th}\)

Spike initiation threshold

tau_syn_ex

0.2 ms

\(\tau_{\mathrm{syn,ex}}\)

Excitatory alpha time constant

tau_syn_in

2.0 ms

\(\tau_{\mathrm{syn,in}}\)

Inhibitory alpha time constant

I_e

0 pA

\(I_\mathrm{e}\)

Constant external current

gsl_error_tol

1e-6

—

RKF45 local error tolerance

V_initializer

Constant(-70.6 mV)

—

Membrane initializer

I_ex_initializer

Constant(0 pA)

—

Excitatory current initializer

I_in_initializer

Constant(0 pA)

—

Inhibitory current initializer

w_initializer

Constant(0 pA)

—

Adaptation current initializer

spk_fun

ReluGrad()

—

Surrogate spike function

spk_reset

'hard'

—

Reset mode (hard matches NEST)

ref_var

False

—

Expose boolean refractory indicator
V#

Membrane potential, shape (*in_size,) with unit mV.

Type:

brainstate.HiddenState

dI_ex#

Excitatory alpha auxiliary state (derivative of \(I_{ex}\)), unit pA/ms.

Type:

brainstate.ShortTermState

I_ex#

Excitatory synaptic current, unit pA.

Type:

brainstate.HiddenState

dI_in#

Inhibitory alpha auxiliary state (derivative of \(I_{in}\)), unit pA/ms.

Type:

brainstate.ShortTermState

I_in#

Inhibitory synaptic current, unit pA.

Type:

brainstate.HiddenState

w#

Adaptation current, unit pA.

Type:

brainstate.HiddenState

refractory_step_count#

Remaining refractory time steps (int32 array).

Type:

brainstate.ShortTermState

integration_step#

Current RKF45 integration step size, unit ms. Persists across simulation steps.

Type:

brainstate.ShortTermState

I_stim#

One-step-delayed external current buffer, unit pA.

Type:

brainstate.ShortTermState

last_spike_time#

Time of last spike emission, unit ms. Updated to \(t + dt\) on spike.

Type:

brainstate.ShortTermState

refractory#

Boolean refractory indicator (only if ref_var=True).

Type:

brainstate.ShortTermState, optional

Raises:

ValueError –

  • If \(V_{reset} \geq V_{peak}\) - If \(\Delta_T < 0\) - If \(V_{peak} < V_{th}\) - If \(C_m \leq 0\) - If \(t_{ref} < 0\) - If any time constant \(\leq 0\) - If gsl_error_tol \(\leq 0\) - If \((V_{peak} - V_{th})/\Delta_T\) would cause exponential overflow - If numerical instability detected during integration (\(V < -1000\) mV or \(|w| > 10^6\) pA)

Notes

NEST Compatibility

  • Replicates NEST aeif_psc_alpha dynamics including RKF45 integration and in-loop spike handling

  • Default parameter values match NEST defaults (converted to SI units)

  • Refractory clamping follows NEST semantics: \(V_{eff} = V_{reset}\) during refractory, with \(dV/dt = 0\)

Numerical Considerations

  • Maximum iteration limit: 100,000 substeps per time step (prevents infinite loops)

  • Minimum step size: \(h_{min} = 10^{-8}\) ms

  • Voltage capping during integration: \(V_{eff} = \min(V, V_{peak})\) outside refractory period to prevent exponential overflow

  • Overflow protection: validates that \(\exp((V_{peak} - V_{th})/\Delta_T)\) remains within floating-point range

Multiple Spikes Per Step

With \(t_{ref} = 0\) (default), a neuron can spike multiple times within a single simulation step. The internal adaptation \(w\) accumulates all spike-triggered increments \(b\), but the returned spike tensor is binary (0 or 1) per step.

Surrogate Gradients

The spk_fun parameter controls backpropagation through spikes for gradient-based learning. The surrogate function approximates the derivative of the Heaviside step function during backward passes.

Examples

Basic usage with default parameters:

>>> import brainpy.state as bst
>>> import saiunit as u
>>> import brainstate as bs
>>>
>>> # Create a population of 100 AdEx neurons
>>> neuron = bst.aeif_psc_alpha(100)
>>>
>>> # Initialize states
>>> with bs.environ.context(dt=0.1 * u.ms):
...     neuron.init_all_states()
>>>
>>> # Simulate one step with external current
>>> with bs.environ.context(dt=0.1 * u.ms):
...     spikes = neuron.update(x=500 * u.pA)
>>> spikes.shape
(100,)

Custom parameters for fast-spiking interneuron:

>>> # Fast-spiking configuration
>>> fs_neuron = bst.aeif_psc_alpha(
...     in_size=50,
...     C_m=150 * u.pF,
...     g_L=20 * u.nS,
...     tau_w=30 * u.ms,
...     a=0 * u.nS,  # Minimal subthreshold adaptation
...     b=20 * u.pA,  # Small spike-triggered adaptation
...     V_th=-52 * u.mV,
...     Delta_T=1.5 * u.mV,
...     tau_syn_ex=0.5 * u.ms,
...     tau_syn_in=1.0 * u.ms,
... )

Connecting to upstream spike sources:

>>> import brainevent as be
>>>
>>> # Create presynaptic spike generator
>>> spike_gen = bst.PoissonSpike(100, freq=10 * u.Hz)
>>>
>>> # Create postsynaptic AdEx neurons
>>> neurons = bst.aeif_psc_alpha(50)
>>>
>>> # Create projection with alpha synapses
>>> proj = be.nn.FixedProb(
...     pre=spike_gen,
...     post=neurons,
...     prob=0.2,
...     weight=50.0,  # pA per spike
... )

See also

aeif_cond_alpha

Conductance-based AdEx with alpha synapses

aeif_psc_exp

AdEx with exponential postsynaptic currents

aeif_psc_delta

AdEx with delta-function synaptic currents

iaf_psc_alpha

Leaky integrate-and-fire with alpha currents (set Delta_T=0)

References

__init__(in_size, V_peak=Quantity(0., "mV"), V_reset=Quantity(-60., "mV"), t_ref=Quantity(0., "ms"), g_L=Quantity(30., "nS"), C_m=Quantity(281., "pF"), E_L=Quantity(-70.6, "mV"), Delta_T=Quantity(2., "mV"), tau_w=Quantity(144., "ms"), a=Quantity(4., "nS"), b=Quantity(80.5, "pA"), V_th=Quantity(-50.4, "mV"), tau_syn_ex=Quantity(0.2, "ms"), tau_syn_in=Quantity(2., "ms"), I_e=Quantity(0., "pA"), gsl_error_tol=1e-06, V_initializer=Constant(value=-70.6 mV), I_ex_initializer=Constant(value=0. pA), I_in_initializer=Constant(value=0. pA), w_initializer=Constant(value=0. pA), spk_fun=ReluGrad(alpha=0.3, width=1.0), spk_reset='hard', ref_var=False, name=None)[source]#

Initialize the aeif_psc_alpha neuron model.

See class docstring for detailed parameter descriptions.

get_spike(V=None)[source]#

Compute differentiable spike output using surrogate gradient.

Applies the surrogate spike function to the scaled membrane potential for gradient-based learning. This method is used for backpropagation and does not affect the internal spike detection logic (which uses hard threshold crossing during integration).

Parameters:

V (ArrayLike, optional) – Membrane potential array with unit mV. If None, uses the current self.V.value. Shape: (*in_size,).

Returns:

spike – Differentiable spike output with shape matching V. Values are continuous in the forward pass (soft spikes) but use surrogate gradients in the backward pass. Typically in range [0, 1] depending on surrogate function.

Return type:

Array

Notes

The voltage is scaled before applying the surrogate function:

\[\begin{split}v_{scaled} = \\frac{V - V_{th}}{V_{th} - V_{reset}}\end{split}\]

This normalization ensures the surrogate function operates in a consistent range regardless of the specific voltage parameters.

init_state(**kwargs)[source]#

Initialize all state variables.

Creates and initializes the membrane potential, synaptic currents, adaptation current, refractory counters, integration step size, and stimulus buffer.

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'))[source]#

Advance the neuron state by one simulation time step.

Performs adaptive RKF45 integration of membrane, synaptic, and adaptation dynamics over the interval \([t, t+dt]\), with in-loop spike detection, reset, and refractory handling matching NEST semantics.

Parameters:

x (ArrayLike, optional) –

External current input at the current time step, with unit pA. Shape must be broadcastable to (*in_size,). Default: 0.0 * u.pA.

This input is stored in the one-step-delayed buffer I_stim and will be used in the next time step’s dynamics (matching NEST input handling).

Returns:

spike – Binary spike indicator with shape (*in_size,), dtype float. Value is 1.0 where at least one spike occurred during the integration interval, 0.0 otherwise.

Note: With t_ref=0, neurons may spike multiple times within the step, but the returned tensor is binary per neuron per step. Internal adaptation dynamics accumulate all spike-triggered increments.

Return type:

Array

Notes

Integration Process

  1. Adaptive RKF45 loop: Starting from current state at time \(t\), integrate ODEs using RKF45 with adaptive step sizing until reaching \(t + dt\).

    • Each substep computes 6 stages (\(k_1\) through \(k_6\))

    • Error estimate: \(err = \max|y_5 - y_4|\)

    • Step acceptance: if \(err \leq atol\) or \(h \leq h_{min}\)

    • Step size update: \(h_{new} = h \cdot \min(5, \max(0.2, 0.9(atol/err)^{0.2}))\)

  2. In-loop spike handling: After each accepted substep, check if \(V \geq V_{peak}\) (or \(V \geq V_{th}\) if \(\Delta_T=0\)). If spike detected:

    • Reset: \(V \leftarrow V_{reset}\)

    • Adaptation jump: \(w \leftarrow w + b\)

    • Refractory counter: \(r \leftarrow \lceil t_{ref}/dt \rceil + 1\) (if enabled)

  3. Post-integration cleanup:

    • Decrement refractory counter: \(r \leftarrow r - 1\) (if \(r > 0\))

    • Deliver synaptic inputs: add spike weights to \(dI_{ex}\) and \(dI_{in}\)

    • Store external input: \(I_{stim} \leftarrow x\) (for next step)

    • Update spike time: \(t_{spike} \leftarrow t + dt\) (where spikes occurred)

Refractory Clamping

During refractory period (\(r > 0\)):

  • Effective voltage: \(V_{eff} = V_{reset}\)

  • Voltage derivative: \(dV/dt = 0\)

  • All other state variables evolve normally

Voltage Capping

Outside refractory period, effective voltage is capped to prevent exponential overflow: \(V_{eff} = \min(V, V_{peak})\).

Numerical Stability

  • Raises ValueError if \(V < -1000\) mV (indicates divergence)

  • Raises ValueError if \(|w| > 10^6\) pA (adaptation overflow)

  • Maximum iteration limit: 100,000 substeps per time step