glif_psc

class brainpy.state.glif_psc(in_size, g=Quantity(9.43, 'nS'), E_L=Quantity(-78.85, 'mV'), V_th=Quantity(-51.68, 'mV'), C_m=Quantity(58.72, 'pF'), t_ref=Quantity(3.75, 'ms'), V_reset=Quantity(-78.85, 'mV'), th_spike_add=0.37, th_spike_decay=0.009, voltage_reset_fraction=0.2, voltage_reset_add=18.51, th_voltage_index=0.005, th_voltage_decay=0.09, asc_init=(0.0, 0.0), asc_decay=(0.003, 0.1), asc_amps=(-9.18, -198.94), asc_r=(1.0, 1.0), tau_syn=(2.0,), spike_dependent_threshold=False, after_spike_currents=False, adapting_threshold=False, I_e=Quantity(0., 'pA'), gsl_error_tol=1e-06, V_initializer=None, spk_fun=ReluGrad(alpha=0.3, width=1.0), spk_reset='hard', ref_var=False, name=None)

Current-based generalized leaky integrate-and-fire (GLIF) neuron model.

The glif_psc model implements the five-level GLIF model hierarchy from the Allen Institute [1], featuring alpha-function shaped synaptic currents, after-spike currents (ASC), spike-dependent threshold adaptation, and voltage-dependent threshold modulation. Exact integration via propagator matrices ensures numerical stability and matches NEST’s implementation.

Model Hierarchy

The five GLIF models are:

• GLIF Model 1 (LIF) — Traditional leaky integrate-and-fire

• GLIF Model 2 (LIF_R) — LIF with biologically defined reset rules

• GLIF Model 3 (LIF_ASC) — LIF with after-spike currents

• GLIF Model 4 (LIF_R_ASC) — LIF with reset rules and after-spike currents

• GLIF Model 5 (LIF_R_ASC_A) — LIF with reset rules, after-spike currents, and a voltage-dependent threshold

Model mechanism selection is based on three boolean parameters:

Model	spike_dependent_threshold	after_spike_currents	adapting_threshold
GLIF1	False	False	False
GLIF2	True	False	False
GLIF3	False	True	False
GLIF4	True	True	False
GLIF5	True	True	True

Mathematical Formulation

1. Membrane Dynamics

The membrane potential \(U\) (stored relative to \(E_L\)) evolves according to exact integration (linear dynamics):

\[U(t+dt) = U(t) \cdot P_{33} + (I_e + I_\mathrm{ASC,sum}) \cdot P_{30} + \sum_k \left( P_{31,k} \cdot y_{1,k} + P_{32,k} \cdot y_{2,k} \right)\]

where the propagator matrix elements are:

\[P_{33} = \exp\left(-\frac{dt}{\tau_m}\right), \quad P_{30} = \frac{\tau_m}{C_m} \left(1 - P_{33}\right), \quad \tau_m = \frac{C_m}{g}\]

and \(P_{31,k}\), \(P_{32,k}\) are computed via the IAFPropagatorAlpha algorithm that handles the singularity when \(\tau_m \approx \tau_{\mathrm{syn},k}\).

2. Synaptic Currents (Alpha Function)

Each receptor port has a current modeled by an alpha function with two state variables \(y_{1,k}\) and \(y_{2,k}\):

\[y_{2,k}(t+dt) = P_{21,k} \cdot y_{1,k}(t) + P_{22,k} \cdot y_{2,k}(t)\]

\[y_{1,k}(t+dt) = P_{11,k} \cdot y_{1,k}(t)\]

where:

\[P_{11,k} = P_{22,k} = \exp(-dt / \tau_{\mathrm{syn},k}), \quad P_{21,k} = dt \cdot P_{11,k}\]

On a presynaptic spike of weight \(w\):

\[y_{1,k} \leftarrow y_{1,k} + w \cdot \frac{e}{\tau_{\mathrm{syn},k}}\]

The alpha function is normalized such that an event of weight 1.0 results in a peak current of 1 pA at \(t = \tau_\mathrm{syn}\).

3. After-Spike Currents (GLIF3/4/5)

After-spike currents (ASC) are modeled as exponentially decaying currents with exact integration. Each ASC component \(I_j\) decays with rate \(k_j\):

\[I_j(t+dt) = I_j(t) \cdot \exp(-k_j \cdot dt)\]

The time-averaged ASC over a step uses the stable coefficient:

\[\bar{I}_j = \frac{1 - \exp(-k_j \cdot dt)}{k_j \cdot dt} \cdot I_j(t)\]

On spike, ASC values are reset:

\[I_j \leftarrow \Delta I_j + I_j \cdot r_j \cdot \exp(-k_j \cdot t_\mathrm{ref})\]

4. Spike-Dependent Threshold (GLIF2/4/5)

The spike component of the threshold decays exponentially:

\[\theta_s(t+dt) = \theta_s(t) \cdot \exp(-b_s \cdot dt)\]

On spike, after refractory decay:

\[\theta_s \leftarrow \theta_s \cdot \exp(-b_s \cdot t_\mathrm{ref}) + \Delta\theta_s\]

Voltage reset (with spike-dependent threshold):

\[U \leftarrow f_v \cdot U_\mathrm{old} + V_\mathrm{add}\]

5. Voltage-Dependent Threshold (GLIF5)

The voltage component of the threshold evolves according to:

\[\theta_v(t+dt) = \phi \cdot (U_\mathrm{old} - \beta) \cdot P_\mathrm{decay} + \frac{1}{P_{\theta,v}} \cdot \left(\theta_v(t) - \phi \cdot (U_\mathrm{old} - \beta) - \frac{a_v}{b_v} \cdot \beta \right) + \frac{a_v}{b_v} \cdot \beta\]

where \(\phi = a_v / (b_v - g/C_m)\), \(P_\mathrm{decay} = \exp(-g \cdot dt / C_m)\), \(P_{\theta,v} = \exp(b_v \cdot dt)\), and \(\beta = (I_e + I_\mathrm{ASC,sum}) / g\).

Overall threshold:

\[\theta = \theta_\infty + \theta_s + \theta_v\]

Spike condition (checked after voltage update):

\[U > \theta\]

6. Numerical Integration and Update Order

NEST uses exact integration for the linear subthreshold dynamics (via propagator matrices). The discrete-time update order per simulation step is:

1. Record \(U_\mathrm{old}\) (relative to \(E_L\)).

2. If not refractory:

   1. Decay spike threshold component.

   2. Compute time-averaged ASC and decay ASC values.

   3. Update membrane potential: \(U = U_\mathrm{old} \cdot P_{33} + (I + ASC_\mathrm{sum}) \cdot P_{30} + \sum P_{31} y_1 + P_{32} y_2\).

   4. Compute voltage-dependent threshold component (using \(U_\mathrm{old}\)).

   5. Update total threshold.

   6. If \(U > \theta\): emit spike, apply reset rules.

3. If refractory: decrement counter, hold U at \(U_\mathrm{old}\).

4. Update synaptic current state variables: \(y_2 = P_{21} y_1 + P_{22} y_2\), then \(y_1 = P_{11} y_1\).

5. Add incoming spike current jumps (scaled by \(e / \tau_\mathrm{syn}\)).

6. Update external current input \(I\).

7. Record and save \(U_\mathrm{old}\) for next step.

Parameters:

• in_size (Size) – Shape of the neuron population. Can be tuple of ints or single int.

• g (ArrayLike, optional) – Membrane (leak) conductance. Default: 9.43 nS.

• E_L (ArrayLike, optional) – Resting membrane potential. Default: -78.85 mV.

• V_th (ArrayLike, optional) – Instantaneous threshold voltage (absolute). Default: -51.68 mV.

• C_m (ArrayLike, optional) – Membrane capacitance. Default: 58.72 pF.

• t_ref (ArrayLike, optional) – Absolute refractory period. Default: 3.75 ms.

• V_reset (ArrayLike, optional) – Reset potential (absolute; used for GLIF1/3). Default: -78.85 mV.

• th_spike_add (float, optional) – Threshold additive constant after spike (mV). Default: 0.37.

• th_spike_decay (float, optional) – Spike threshold decay rate (/ms). Default: 0.009.

• voltage_reset_fraction (float, optional) – Voltage fraction coefficient after spike. Default: 0.20.

• voltage_reset_add (float, optional) – Voltage additive constant after spike (mV). Default: 18.51.

• th_voltage_index (float, optional) – Voltage-dependent threshold leak rate (/ms). Default: 0.005.

• th_voltage_decay (float, optional) – Voltage-dependent threshold decay rate (/ms). Default: 0.09.

• asc_init (Sequence[float], optional) – Initial values of after-spike currents (pA). Default: (0.0, 0.0).

• asc_decay (Sequence[float], optional) – ASC decay rates (/ms). Default: (0.003, 0.1).

• asc_amps (Sequence[float], optional) – ASC amplitudes added on spike (pA). Default: (-9.18, -198.94).

• asc_r (Sequence[float], optional) – ASC fraction coefficients (dimensionless). Default: (1.0, 1.0).

• tau_syn (Sequence[float], optional) – Synaptic alpha-function time constants (ms), one per receptor port. Default: (2.0,).

• spike_dependent_threshold (bool, optional) – Enable biologically defined reset rules (GLIF2/4/5). Default: False.

• after_spike_currents (bool, optional) – Enable after-spike currents (GLIF3/4/5). Default: False.

• adapting_threshold (bool, optional) – Enable voltage-dependent threshold (GLIF5). Default: False.

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

• gsl_error_tol (ArrayLike, optional) – Unitless local RKF45 error tolerance, broadcastable and strictly positive. Default: 1e-6.

• V_initializer (Callable, optional) – Membrane potential initializer. Default: Constant(E_L).

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

• spk_reset (str, optional) – Spike reset mode: ‘hard’ or ‘soft’. Default: ‘hard’.

• ref_var (bool, optional) – If True, allocate and expose self.refractory state.

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

Parameter Mapping

Parameter	Default	Math equivalent	Description
in_size	(required)		Population shape
g	9.43 nS	\(g\)	Membrane (leak) conductance
E_L	-78.85 mV	\(E_L\)	Resting membrane potential
V_th	-51.68 mV	\(V_\mathrm{th}\)	Instantaneous threshold (absolute)
C_m	58.72 pF	\(C_\mathrm{m}\)	Membrane capacitance
t_ref	3.75 ms	\(t_\mathrm{ref}\)	Absolute refractory period
V_reset	-78.85 mV	\(V_\mathrm{reset}\)	Reset potential (absolute; GLIF1/3)
th_spike_add	0.37 mV	\(\Delta\theta_s\)	Threshold additive constant after spike
th_spike_decay	0.009 /ms	\(b_s\)	Spike threshold decay rate
voltage_reset_fraction	0.20	\(f_v\)	Voltage fraction after spike
voltage_reset_add	18.51 mV	\(V_\mathrm{add}\)	Voltage additive after spike
th_voltage_index	0.005 /ms	\(a_v\)	Voltage-dependent threshold leak
th_voltage_decay	0.09 /ms	\(b_v\)	Voltage-dependent threshold decay rate
asc_init	(0.0, 0.0) pA		Initial values of ASC
asc_decay	(0.003, 0.1) /ms	\(k_j\)	ASC time constants (decay rates)
asc_amps	(-9.18, -198.94) pA	\(\Delta I_j\)	ASC amplitudes on spike
asc_r	(1.0, 1.0)	\(r_j\)	ASC fraction coefficient
tau_syn	(2.0,) ms	\(\tau_{\mathrm{syn},k}\)	Synaptic alpha-function time constants
spike_dependent_threshold	False		Enable biologically defined reset (GLIF2/4/5)
after_spike_currents	False		Enable after-spike currents (GLIF3/4/5)
adapting_threshold	False		Enable voltage-dependent threshold (GLIF5)
I_e	0.0 pA	\(I_e\)	Constant external current
gsl_error_tol	1e-6	–	Local absolute tolerance for RKF45 error estimate
V_initializer	Constant(E_L)		Membrane potential initializer
spk_fun	ReluGrad()		Surrogate spike function
spk_reset	'hard'		Reset mode
ref_var	False		If True, expose boolean refractory state

V

Membrane potential \(V_\mathrm{m}\) (absolute, mV).

Type:

HiddenState

y1

Synaptic current derivative states (pA), one per receptor port.

Type:

list of HiddenState

y2

Synaptic current states (pA), one per receptor port.

Type:

list of HiddenState

last_spike_time

Last spike time for each neuron (ms).

Type:

ShortTermState

refractory_step_count

Remaining refractory grid steps (int32).

Type:

ShortTermState

integration_step

Persistent RKF45 substep size estimate (ms).

Type:

ShortTermState

I_stim

Buffered external current for next step (pA).

Type:

ShortTermState

_ASCurrents

After-spike current values (pA). Shape: (n_asc, *varshape).

Type:

numpy.ndarray

_ASCurrents_sum

Sum of after-spike currents (pA). Shape: (*varshape).

Type:

numpy.ndarray

_threshold

Total threshold (relative to E_L, in mV). Shape: (*varshape).

Type:

numpy.ndarray

_threshold_spike

Spike component of threshold (mV). Shape: (*varshape).

Type:

numpy.ndarray

_threshold_voltage

Voltage component of threshold (mV). Shape: (*varshape).

Type:

numpy.ndarray

refractory

Optional boolean refractory indicator, available only when ref_var=True.

Type:

ShortTermState

Raises:

• ValueError – If invalid model mechanism combination is specified.

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

• ValueError – If capacitance, conductance, or time constants are not positive.

• ValueError – If voltage_reset_fraction not in [0, 1].

• ValueError – If asc_r values not in [0, 1].

• ValueError – If ASC parameter arrays have mismatched lengths.

Notes

• Default parameter values are from GLIF Model 5 of Cell 490626718 from the Allen Cell Type Database.

• Parameters V_th and V_reset are specified in absolute mV. Internally, membrane potential is tracked relative to E_L, matching NEST’s convention.

• For models with spike-dependent threshold (GLIF2/4/5), the reset condition should satisfy:

  \[E_L + f_v \cdot (V_{th} - E_L) + V_{add} < V_{th} + \Delta\theta_s\]

  Otherwise the neuron may spike continuously.

• Unlike glif_cond which uses an RKF45 ODE integrator, glif_psc uses exact integration via propagator matrices for the linear subthreshold dynamics, matching NEST’s implementation.

• If tau_m is very close to tau_syn, the model numerically behaves as if they are equal, to avoid numerical instabilities (see NEST IAF_Integration_Singularity notebook).

• Synaptic inputs are delivered to receptor ports starting from port 0. Register inputs with keys like ‘receptor_0’, ‘receptor_1’, etc., via the add_delta_input method. Inputs without a receptor label default to receptor port 0.

Examples

GLIF Model 1 (Basic LIF):

>>> import brainpy.state as bp
>>> import saiunit as u
>>> with u.context(dt=0.1 * u.ms):
...     model = bp.glif_psc(100, spike_dependent_threshold=False,
...                         after_spike_currents=False, adapting_threshold=False)
...     model.init_all_states()
...     output = model(350 * u.pA)

GLIF Model 5 (Full Model with Adaptation):

>>> import brainpy.state as bp
>>> import saiunit as u
>>> with u.context(dt=0.1 * u.ms):
...     model = bp.glif_psc(100, spike_dependent_threshold=True,
...                         after_spike_currents=True, adapting_threshold=True)
...     model.init_all_states()
...     output = model(200 * u.pA)

Multi-Receptor Configuration:

>>> import brainpy.state as bp
>>> import saiunit as u
>>> with u.context(dt=0.1 * u.ms):
...     model = bp.glif_psc(100, tau_syn=(2.0, 5.0, 10.0))
...     model.init_all_states()
...     # Register inputs to different receptor ports
...     model.add_delta_input('exc_receptor_0', lambda: 10 * u.pA)
...     model.add_delta_input('inh_receptor_1', lambda: -5 * u.pA)

References

[1]

Teeter C, Iyer R, Menon V, Gouwens N, Feng D, Berg J, Szafer A, Cain N, Zeng H, Hawrylycz M, Koch C, & Mihalas S (2018). Generalized leaky integrate-and-fire models classify multiple neuron types. Nature Communications 9:709.

[2]

Meffin H, Burkitt AN, Grayden DB (2004). An analytical model for the large, fluctuating synaptic conductance state typical of neocortical neurons in vivo. J. Comput. Neurosci. 16:159-175.

[3]

NEST Simulator glif_psc model documentation and C++ source: models/glif_psc.h and models/glif_psc.cpp.

See also

glif_cond

Conductance-based GLIF model with RKF45 integration.

gif_psc_exp_multisynapse

Generalized IF with exponential synapses.

get_spike(V=None)

Generate spike output via surrogate gradient function.

Applies the surrogate gradient function to a normalized voltage signal. The voltage is linearly scaled such that V_th maps to 1 and V_reset maps to 0, providing a normalized input for the surrogate function.

Parameters:

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

Returns:

spike – Spike output (float32). Shape matches the neuron population. Forward pass produces values in [0, 1]; backward pass uses the surrogate gradient specified by spk_fun.

Return type:

jax.numpy.ndarray

Notes

• This method is called internally by the base Neuron class and is typically not invoked directly by users.

• The surrogate function enables gradient-based learning by providing a differentiable approximation to the Heaviside step function.

init_state(**kwargs)

Initialize persistent and short-term state variables.

All GLIF-specific state variables are stored as JAX HiddenState arrays, and pre-computed decay constants are stored as Python floats.

Parameters:

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

property n_receptors

Number of synaptic receptor ports.

Returns:

Number of independent receptor ports, determined by the length of the tau_syn parameter. Each receptor port has its own synaptic time constant and independent alpha-function dynamics.

Return type:

int

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

Perform a single simulation step using exact propagator matrices.

Implements the NEST glif_psc update using the exact IAFPropagatorAlpha integration scheme. All GLIF-specific discrete updates (threshold decay, ASC, voltage-dependent threshold) are applied as vectorised JAX operations, making this method compatible with brainstate.transform.for_loop.

Parameters:

x (ArrayLike, optional) – External current input (pA), applied with one-step delay. Default: 0.0 pA.

Returns:

spike – Binary spike tensor (float32), shape (*varshape).

Return type:

jax.Array